Skip to content
Snippets Groups Projects
Commit 11d64388 authored by Eric Kooistra's avatar Eric Kooistra
Browse files

Addd more description.

parent 24315517
Branches
No related tags found
1 merge request!420Resolve RTSD-264
%% Cell type:markdown id:93474f68 tags:
 
# Multirate mixer
 
Author: Eric Kooistra, May 2024
Author: Eric Kooistra, May - Dec 2024
 
Purpose:
* Practise DSP [1].
* Use multirate processing to implement a single channel mixer
* Extend single channel mixer to DFT filterbank
 
Description:
 
For a real input the mixer can either downconvert the positive or the negative frequeny band. Convention is to downconvert the positive frequency band [2, 3]. This leads to using IDFT in the analysis DFT filterbank and IDFT in the synthesis filterbank [2, 3]. By using IDFT(x) = DFT(fold(x)) it is possible to change to a DFT in the analysis filterbank, this corresponds to using a clockwise commutator [3] and type 3 polyphase filter structure [4].
 
Sampling and anlogue signal yields the low pass spectrum and its replicas around multiples of fs. This endles replication is why the discrete samples can represent and anlogue signal. Downsampling of the digitized signal again causes aliasing of all replicas to baseband. Normally the interest is in the baseband replica. However using a BPF instead of an LPF any of the replicas can be selected for baseband. Similar for upsampling (synthesis) the baseband signal can be upconverted to the wanted replica by using a BPF instead of an LPF. The first step in upsampling is inserting zeros to increase the sample rate, without changing the digital signal. However this increased sample rate does make the replicas available in the digital domain.
With a single channel mixer the downsampling by M implies that at the downsampled rate there is processing power to process M - 1 more channels. This fits a critically sampled DFT filterbank with K = M. Whereby the DFT should typically be an FFT.
The single channel mixer and DFT filterbank can share a prototype LPF for all channels, because the downsampling and upsampling performs the demodulation and modulation, to and from baseband at 0 Hz. When Ros = K / M > 1 then there remains and offset frequency that is compensated by a phasor term W_K^(knM). With the DFT filterbank this phasor that can be implemented as a circular shift (modulo K) at the input of the DFT.
For perfect reconstruction by an analysis-synthesis filterbank the aliasing must cancel and distortion must be zero. Considering only adjancent channels is called pseudo QMF (two-channel), for the other channels the stop band attenuation suppresses the crosstalk. For zero distortion the transfer function of the channels must be power complementary. With an DFT filterbank (= complex modulated filterbank) the aliasing does not cancel. For critically sampled filterbanks this leads to the Modified DFT (MDFT) or Cosine modulated fiterbank [4]. By using oversampling the aliasing can be avoided. Hence using some oversampling and power complementary channel filters the oversampled DFT filterbank can achieve almost perfect reconstruction.
Remarks:
* With oversampling the subband bandwidth can be increased to have a flat reponse up to fsub / 2. The transition band is then from fsub / 2 to Ros * fsub / 2 The subband filter is then not power complementary (around fsub / 2), so the subbands are not suitable for synthesis. It may be feasible achieve almost perfect reconstruction by first removing the transition band from the subbands, by separating the subbands into fine channels using a DFT, then zero the fine channels of the transistion band, and then synthesize the subband using IDFT.
Features:
* Use full rate model of single channel down converter and up converter as expected exact golden reference result for the efficient polyphase implementation and WOLA implementation and refBunton.
* Use SNR of input and difference with time aligned reconstructed output, to show that perfect reconstruction depends on the pass band gain being one and the stopband gain approaching zero. Hence the SNR for center bin sine wave inputs improves, towards perfect reconstruction, when Ntaps is increased.
* Use analysis DFT filterbank and synthesis IDFT filterbank and verify that it yields the same result as with the single channel pipeline.
* Compare model with PFB reconstruction by reconstruct.m of John Bunton for Ntaps = 12 and Ndft = 192, using refBunton.
* Support Ncoefs <= Ndelays = Ntaps * Ndft, to verify correct implementation of coefs and delay line
* Support WG with offset bin center frequency or with noise, to verify that result is correct for any frequency, not only the center frequency
* Support Ncoefs <= Ndelays = Ntaps * Ndft, to verify correct implementation of coefs and delay line.
* Support WG with offset bin center frequency (e.g. wgSub = 1.5) or with noise (via SNR_WG_dB < 100), to verify that result is correct for any frequency, not only the center frequency.
* Support refBunton with asymmetrical hPrototype, to verify correct order implementation of coefs in analysis and synthesis. The refBunton has slightly asymmetrical hPrototype due to that it uses interpolation to increase the number of coefficient after designing the FIR filter with fircls1().
 
References:
1. dsp_study_erko, summary of DSP books
2. chapter 7 in [CROCHIERE]
3. chapter 6 downconverter, 7 upconverter, 9 filterbank in [HARRIS]
4. chapter 1.1.2, 8.6 in [FLIEGE]
4. chapter 1.1.2 (commutator), 7.3 (pseudo QMF), 7.3.3 (MDFT), 7.4 (Cos filterbank) 8.6 in [FLIEGE]
 
%% Cell type:code id:8043fa7b tags:
 
``` python
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from fractions import Fraction
```
 
%% Cell type:code id:1c1ba454 tags:
 
``` python
from sys import maxsize
np.set_printoptions(threshold=maxsize)
```
 
%% Cell type:code id:fc530dbc tags:
 
``` python
# Auto reload module when it is changed
%load_ext autoreload
%autoreload 2
 
# Add rtdsp module path to $PYTHONPATH
import os
import sys
module_path = os.path.abspath(os.path.join('../'))
if module_path not in sys.path:
sys.path.insert(0, module_path)
 
# Import rtdsp
from rtdsp.utilities import ceil_div, verify_result, is_integer_value, is_symmetrical, pow_db, snr_db, \
read_coefficients_file
from rtdsp.firfilter import filterbank_frequency_response
from rtdsp.fourier import dtft
from rtdsp.multirate import down, up
from rtdsp.singlechannel import unit_circle_loops_phasor_arr, \
maximal_downsample_bpf, non_maximal_downsample_bpf, \
maximal_upsample_bpf, non_maximal_upsample_bpf
from rtdsp.dftfilterbank import analysis_dft_filterbank, synthesis_dft_filterbank
from rtdsp.plotting import plot_power_spectrum, plot_magnitude_spectrum, plot_phase_spectrum
```
 
%% Cell type:code id:f527cbb6 tags:
 
``` python
# Verification control
enExit = True # Exit or continue when a verification fails
vb = 0 # Verbosity
```
 
%% Cell type:markdown id:fbf4bf7b tags:
 
# 0 Time, bins and rates
 
%% Cell type:code id:5b37a1dc tags:
 
``` python
### Filterbank
refBunton = True
#refBunton = False
if refBunton:
Ntaps = 12
Ndft = 192
# To compare SNR results with reconstruct.m
normalizeDcGain = False
flipAPrototype = True
# To compare SNR performance with generic hPrototype
#normalizeDcGain = True
else:
Ntaps = 8 # number of taps per polyphase FIR filter
Ndft = 16 # DFT size
#Ntaps = 4
#Ndft = 8
Ntaps = 12
Ndft = 192
#Ntaps = 12
#Ndft = 192
Ntaps = 32
# Generic hPrototype is symmetrical, so then coefs flip makes no difference
flipAPrototype = False
Ndelays = Ndft * Ntaps
Ncoefs = Ndelays
#Ncoefs = Ndelays - 2 # try shorter length to verify impact of hPairGroupDelay
#Ncoefs = Ndelays - 1 # try odd length
print('Ntaps =', Ntaps)
print('Ndft =', Ndft)
print('Ndelays =', Ndelays)
print('Ncoefs =', Ncoefs)
 
# Waveform generator
wgSub = 1.5 # in range(Nsub)
wgPhase = 0.0 # in degrees
wgModulation = 0 # for amplitude modulation (AM) frequency fsub / wgModulation, use 0 for no AM
# use >> 1 for fraction of fsub
print()
print('wgSub =', wgSub)
print('wgPhase =', wgPhase)
print('wgModulation =', wgModulation)
 
# Additive Gaussian White Noise (AGWN), use SNR_WG_dB >= 100 to have only WG carrier,
# so no noise and SNR = +inf
SNR_WG_dB = 100
 
# Mixer LO bin
kLo = int(np.round(wgSub))
print('kLo =', kLo)
 
# Downsample rate for analysis
if refBunton:
Ndown = 144
else:
Ndown = Ndft * 3 // 4 # oversampled PFB
#Ndown = Ndft * 7 // 8
#Ndown = Ndft * 15 // 16
#Ndown = Ndft // 4
#Ndown = Ndft # Critically sampled PFB
Ros = Fraction(Ndft, Ndown)
print()
print('Ros =', Ros)
print('. Ndown =', Ndown)
print('. Ndft =', Ndft)
 
# Upsample rate (= downsample rate) for synthesis
Nup = Ndown
```
 
%% Output
 
Ntaps = 12
Ndft = 192
Ndelays = 2304
Ncoefs = 2304
wgSub = 1.5
wgPhase = 0.0
wgModulation = 0
kLo = 2
Ros = 4/3
. Ndown = 144
. Ndft = 192
 
%% Cell type:code id:e5680c7b tags:
 
``` python
# Samples
fs = 1.0 # sample rate
Ts = 1 / fs # sample period
```
 
%% Cell type:code id:74ca764f tags:
 
``` python
# Subbands
Nsub = Ndft // 2 # number of subbands in fs / 2
fsub = fs / Ndft # subband center frequency
Tsub = 1 / fsub # subband period
 
fdown = fs / Ndown # downsampled data rate
Tdown = 1 / fdown # downsampled data period
 
# Two single side bands, one keeping only downconverted positive subband frequencies
nofSsb = 1 if kLo == 0 else 2
 
print('Tsub =', Tsub)
print('Tdown =', Tdown)
```
 
%% Output
 
Tsub = 192.0
Tdown = 144.0
 
%% Cell type:code id:786af296 tags:
 
``` python
# Time in number of subband periods to simulate
# . use small Nsim for more detail in point plots,
# . use large Nsim for more accurracy in SNR estimates or when wgSub < 0.5,
# . force Nsim as multiple of Ndown to avoid mismatch in array lengths.
Nsim = ceil_div(3 * Ntaps, Ndown) * Ndown
Nsim *= 4
print('Nsim = %6d [Tsub]' % Nsim)
 
# Time in number input samples to simulate
Nsamples = Nsim * Ndft
# . input time index n for up rate
n_i = np.arange(Nsamples) # sample index, time in sample period units [Ts]
n_s = n_i * Ts # time in seconds
n_sub = n_s / Ndft # time in subband period units [Tsub]
 
# Time in number of downsampled samples
Msamples = Nsamples // Ndown
# . downsampled time index m for down rate, n = m D, so m = n // D
m_i = np.arange(Msamples) # downsampled sample index
m_s = down(n_s, Ndown) # = m_i * Tdown, time in seconds
m_sub = m_s / Ndft # time in subband period units [Tsub]
 
print('Nsamples = %6d [Ts]' % Nsamples)
print('Msamples = %6d [Tdown]' % Msamples)
```
 
%% Output
 
Nsim = 576 [Tsub]
Nsamples = 110592 [Ts]
Msamples = 768 [Tdown]
 
%% Cell type:markdown id:8cd835d6 tags:
 
# 1 Prototype FIR low pass filter
 
%% Cell type:code id:1bb76ada tags:
 
``` python
if refBunton:
# Load LPF prototype FIR filter coefficients from GenPfilter12.m. These coefficients
# are not symmetrical (so not exactly linear phase), which makes it possible to verify
# whether the coefficients are applied in the correct order.
hPrototype = read_coefficients_file('../pfb_bunton_annotated/cl.txt')
if Ncoefs < Ndelays:
dDiv = (Ndelays - Ncoefs) // 2
dRem = (Ndelays - Ncoefs) % 2
hPrototype = hPrototype[dDiv : -dDiv - dRem]
if normalizeDcGain:
hPrototype /= np.sum(hPrototype)
print('refBunton: is_symmetrical(hPrototype) =', is_symmetrical(hPrototype))
else:
# Use windowed sync (= firwin) prototype FIR filter
# . For sinc() the ideal bandwidth is 2pi / Ndft = fs / Ndft = fsub,
# . Use half power bandwidth factor hpFactor to tune half power cutoff frequency of LPF.
# . Default hpFactor = 1.0 yields flat filterbank aggregate frequency response for
# HFbank firwin hanning filter, but for perfect reconstruction HFpowerbank needs to
# be flat.
# Optimum flat HFpowerbank
if Ntaps == 8 and Ndft == 16:
if Ndft == 16 and Ntaps == 8:
hpFactor = 1.108
elif Ntaps == 12 and Ndft == 192:
elif Ndft == 16 and Ntaps == 32:
hpFactor = 1.027
elif Ndft == 192 and Ntaps == 12:
hpFactor = 1.065
else:
# Default for optimum flat HFbank
hpFactor = 1.0
BWbin = fs / Ndft # bandwidth of one bin
BWpass = hpFactor * BWbin
fpass = BWpass / 2 # bin at DC: -fpass to +fpass
fcutoff = fpass
hPrototype = signal.firwin(Ncoefs, fcutoff, window='hann', fs=fs)
verify_result(is_symmetrical(hPrototype), ': is_symmetrical(hPrototype)', enExit)
print('DC gain = %.10f' % np.sum(hPrototype))
```
 
%% Output
 
refBunton: is_symmetrical(hPrototype) = False
DC gain = 1.0051947660
 
%% Cell type:code id:6ebc94aa tags:
 
``` python
# Plot impulse response
plt.figure(1)
plt.plot(hPrototype, 'g')
plt.title('Impulse response')
plt.grid(True)
```
 
%% Output
 
 
%% Cell type:markdown id:4d9b4324 tags:
 
 
%% Cell type:code id:03da0b30 tags:
 
``` python
# Single FIR filter group delay
hGroupDelay = (Ncoefs - 1) / 2
 
# Group delay for LPF down and LPF up in series, is integer for Ncoefs is even and odd
hPairGroupDelay = Ncoefs - 1 # = 2 * hGroupDelay
 
print('hGroupDelay = %.1f samples' % hGroupDelay)
print('hPairGroupDelay = %d samples' % hPairGroupDelay)
```
 
%% Output
 
hGroupDelay = 1151.5 samples
hPairGroupDelay = 2303 samples
 
%% Cell type:code id:3abeee86 tags:
 
``` python
# Filterbank aggregate frequency response
h, f, HFprototype = dtft(hPrototype)
print('DC gain = %.10f' % np.abs(HFprototype[len(HFprototype) // 2]))
HFbank = filterbank_frequency_response(HFprototype, Ndft)
 
# Magnitude squared response, to have correlator (power) response.
HFpowerbank = filterbank_frequency_response(HFprototype**2, Ndft)
 
# Filterbank bin 1, 2 frequency responses, HFprototype is for bin 0
Lprototype = len(HFprototype)
Lbin = Lprototype // Ndft
HF1 = np.roll(HFprototype, 1*Lbin)
HF2 = np.roll(HFprototype, 2*Lbin)
```
 
%% Output
 
DC gain = 1.0051947660
 
%% Cell type:code id:8008d6b7 tags:
 
``` python
# Plot transfer function (use frequency in fsub units)
fLim = None
dbLim = None
voltLim = None
 
plt.figure(1)
fLim = (0, 2)
if Ntaps <= 32:
dbLim = (-120, 5)
plot_power_spectrum(f, HFprototype, 'r--', fs / fsub, fLim, dbLim) # bin 0
plot_power_spectrum(f, HF1, 'b--', fs / fsub, fLim, dbLim) # bin 1
plot_power_spectrum(f, HF2, 'm--', fs / fsub, fLim, dbLim) # bin 2
plot_power_spectrum(f, HFbank, 'g', fs / fsub, fLim, dbLim) # all bins
 
plt.figure(2)
#fLim = (0, 1)
voltLim = (0.9, 1.1)
plot_magnitude_spectrum(f, HFprototype, 'r--', fs / fsub, fLim, voltLim) # bin 0
plot_magnitude_spectrum(f, HF1, 'b--', fs / fsub, fLim, voltLim) # bin 1
plot_magnitude_spectrum(f, HF2, 'm--', fs / fsub, fLim, voltLim) # bin 2
plot_magnitude_spectrum(f, HFbank, 'g--', fs / fsub, fLim, voltLim) # all bins
plot_magnitude_spectrum(f, HFpowerbank, 'g', fs / fsub, fLim, voltLim) # all bins
```
 
%% Output
 
 
 
%% Cell type:code id:372445f4 tags:
 
``` python
fLim = None
#fLim = [-2.1, 2.1]
aLim = None
plot_phase_spectrum(f, HFprototype, fmt='g', fs=fs / fsub, fLim=fLim, aLim=aLim)
```
 
%% Output
 
 
%% Cell type:code id:2c7b858c tags:
 
``` python
# FIR filter group delay
w, gd = signal.group_delay((hPrototype, [1.0]), w=1024, fs=1)
#w, gd = signal.group_delay((np.convolve(hPrototype, hPrototype), [1.0]), w=1024, fs=1)
#w, gd = signal.group_delay((np.convolve(hPrototype, np.flip(hPrototype)), [1.0]), w=1024, fs=1)
#plt.title('FIR filter group delay')
plt.plot(w, gd)
plt.ylabel('Group delay [samples]')
plt.xlabel('Frequency [rad/sample]')
if refBunton:
plt.ylim([0, 5000])
print('FIR filter group delay = %f' % hGroupDelay)
```
 
%% Output
 
FIR filter group delay = 1151.500000
 
/tmp/ipykernel_45544/4263550011.py:2: UserWarning: The filter's denominator is extremely small at frequencies [0.270, 0.295, 0.319, 0.344, 0.368, 0.393, 0.417, 0.442, 0.466, 0.491, 0.515, 0.540, 0.565, 0.589, 0.614, 0.638, 0.663, 0.687, 0.712, 0.736, 0.761, 0.785, 0.810, 0.834, 0.859, 0.884, 0.908, 0.933, 0.957, 0.982, 1.006, 1.031, 1.055, 1.080, 1.104, 1.129, 1.154, 1.178, 1.203, 1.227, 1.252, 1.276, 1.301, 1.325, 1.350, 1.374, 1.399, 1.424, 1.448, 1.473, 1.497, 1.522, 1.546, 1.571, 1.595, 1.620, 1.644, 1.669, 1.694, 1.718, 1.743, 1.767, 1.792, 1.816, 1.841, 1.865, 1.890, 1.914, 1.939, 1.963, 1.988, 2.013, 2.037, 2.062, 2.086, 2.111, 2.135, 2.160, 2.184, 2.209, 2.233, 2.258, 2.283, 2.307, 2.332, 2.356, 2.381, 2.405, 2.430, 2.454, 2.479, 2.503, 2.528, 2.553, 2.577, 2.602, 2.626, 2.651, 2.675, 2.700, 2.724, 2.749, 2.773, 2.798, 2.823, 2.847, 2.872, 2.896, 2.921, 2.945, 2.970, 2.994, 3.019, 3.043, 3.068, 3.093, 3.117], around which a singularity may be present
/tmp/ipykernel_27035/4263550011.py:2: UserWarning: The filter's denominator is extremely small at frequencies [0.270, 0.295, 0.319, 0.344, 0.368, 0.393, 0.417, 0.442, 0.466, 0.491, 0.515, 0.540, 0.565, 0.589, 0.614, 0.638, 0.663, 0.687, 0.712, 0.736, 0.761, 0.785, 0.810, 0.834, 0.859, 0.884, 0.908, 0.933, 0.957, 0.982, 1.006, 1.031, 1.055, 1.080, 1.104, 1.129, 1.154, 1.178, 1.203, 1.227, 1.252, 1.276, 1.301, 1.325, 1.350, 1.374, 1.399, 1.424, 1.448, 1.473, 1.497, 1.522, 1.546, 1.571, 1.595, 1.620, 1.644, 1.669, 1.694, 1.718, 1.743, 1.767, 1.792, 1.816, 1.841, 1.865, 1.890, 1.914, 1.939, 1.963, 1.988, 2.013, 2.037, 2.062, 2.086, 2.111, 2.135, 2.160, 2.184, 2.209, 2.233, 2.258, 2.283, 2.307, 2.332, 2.356, 2.381, 2.405, 2.430, 2.454, 2.479, 2.503, 2.528, 2.553, 2.577, 2.602, 2.626, 2.651, 2.675, 2.700, 2.724, 2.749, 2.773, 2.798, 2.823, 2.847, 2.872, 2.896, 2.921, 2.945, 2.970, 2.994, 3.019, 3.043, 3.068, 3.093, 3.117], around which a singularity may be present
w, gd = signal.group_delay((hPrototype, [1.0]), w=1024, fs=1)
 
 
%% Cell type:code id:fa4b54d6 tags:
 
``` python
# Select impulse response for downconverter / analysys and for upconverter / synthesis
# . Use flipped aPrototype to make SNR results for refBunton with CW equal to reconstruct.m.
# The SNR improves slightly for flipped aPrototype, because refBunton hPrototype is
# asymmetrical. The combination of convolute(np.flip(hPrototype), hPrototype) yields an
# overall linear phase, so constant group delay for the reconstructed signal.
if flipAPrototype:
aPrototype = np.flip(hPrototype)
else:
aPrototype = hPrototype
sPrototype = hPrototype
#hPrototype
```
 
%% Cell type:markdown id:ce7e672e tags:
 
# 2 Generate input data
 
%% Cell type:code id:e74340e4 tags:
 
``` python
# x[n] = carrier waveform
# . freq = center subband yields constant baseband I, Q signal
# . phase = 0 yields Q = 0
wgSubIsInt = is_integer_value(wgSub)
wgFreq = wgSub * fsub # in Hz
wgAmpl = 1
xData = wgAmpl * np.cos(2 * np.pi * wgFreq * n_s + np.radians(wgPhase))
#xData = n_i
# Apply some modulation
if wgModulation:
xData *= 1 + 0.5 * np.cos(2 * np.pi * fsub / wgModulation * n_s)
```
 
%% Cell type:code id:bd2add56 tags:
 
``` python
# Add AGWN when SNR_WG_dB < 100
seed = None
seed = 1
rng = np.random.default_rng(seed)
mu = 0.0
sigma = wgAmpl * np.sqrt(0.5) / 10**(SNR_WG_dB / 20)
noise = rng.normal(mu, sigma, Nsamples)
if SNR_WG_dB < 100:
xData += noise
# Check SNR
powCarrier = np.sum(xData**2) / Nsamples
powNoise = np.sum(noise**2) / Nsamples
snrWg_db = 10 * np.log10(powCarrier / powNoise)
print('powCarrier = %f' % powCarrier)
print('powNoise = %f' % powNoise)
print('snrWg_db = %f dB' % snrWg_db)
```
 
%% Cell type:code id:7106ad3f tags:
 
``` python
plt.plot(n_sub, xData, 'b.-')
plt.xlim([0,10])
```
 
%% Output
 
(0.0, 10.0)
 
 
%% Cell type:markdown id:c6f01b91 tags:
 
# 3 Single channel downconverter
 
%% Cell type:markdown id:0ad89dcf tags:
 
# 3.1 Full rate: LO --> LPF --> D
 
Downconvert bin kLo to baseband, then LPF still at sample rate and then downsample [HARRIS Fig 6.2].
 
%% Cell type:code id:da53b25e tags:
 
``` python
# Down mixer local oscillator (LO) for channel k
w_k = 2 * np.pi * kLo / Ndft
LO = np.exp(-1j * w_k * n_s)
```
 
%% Cell type:code id:9acf0ec2 tags:
 
``` python
# x[n] --> LO --> LPF --> D --> y[mD, k] [HARRIS Fig 6.2]
xLoData = xData * LO
yData = signal.lfilter(aPrototype, [1.0], xLoData) # = y[n, k], Eq. 6.1
yDown = down(yData, Ndown) # = y[mD, k]
 
plt.plot(n_sub, yData.real, 'g-')
plt.plot(n_sub, yData.imag, 'g--')
plt.plot(m_sub, yDown.real, 'b.')
plt.plot(m_sub, yDown.imag, 'b.')
plt.title('Baseand signal for channel %d' % kLo)
plt.xlabel('Time [Tsub]')
plt.ylabel('yData')
#plt.ylim([-0.6, 0.6])
plt.legend(['y real', 'y imag'])
```
 
%% Output
 
<matplotlib.legend.Legend at 0x7f1042f35fd0>
<matplotlib.legend.Legend at 0x7fb4e035dee0>
 
 
%% Cell type:markdown id:6433c3f4 tags:
 
## 3.2 LO at downsampled rate: BPF --> D --> LOdown
 
Use BPF centered at kLo (is LPF shifted by +kLo) still at sample rate, then downsample and do down conversion from kLo to baseband at downsampled rate [HARRIS Fig 6.7].
 
If Ndown = Ndft, then D * w_k = D * 2pi * k / Ndft is multiple of 2pi, so then LOdown = 1.
 
%% Cell type:code id:b8036250 tags:
 
``` python
# LO --> D --> LOdown == loD
#
# LOdown = exp(-j * D w_k * m) = loD
# = exp(-j * 2 pi D / Ndft * k * m)
LOdown = down(LO, Ndown)
 
D_w_k = Ndown * w_k
loD = np.exp(-1j * D_w_k * m_i)
 
print('Ndft =', Ndft)
print('Ndown =', Ndown)
print('w_k =', w_k)
print('D_w_k =', D_w_k)
print('')
 
# Verify that LO data rotates with w_k and LO down with D * w_k rad/s
result = np.all(np.isclose(LOdown, loD))
if not result:
plt.plot(m_sub, LOdown.real, 'r-')
plt.plot(m_sub, LOdown.imag, 'r--')
plt.plot(m_sub, loD.real, 'g-')
plt.plot(m_sub, loD.imag, 'g--')
verify_result(result, ': LOdown == LoD', enExit)
```
 
%% Output
 
Ndft = 192
Ndown = 144
w_k = 0.06544984694978735
D_w_k = 9.42477796076938
PASSED: LOdown == LoD
 
%% Cell type:code id:0663df66 tags:
 
``` python
# Verify that LOdown == 1 when Ndown == Ndft
if Ndown == Ndft:
result = np.all(np.isclose(LOdown, 1.0))
if not result:
plt.plot(m_sub, LOdown.real, 'r-')
plt.plot(m_sub, LOdown.imag, 'r--')
verify_result(result, ': LOdown == 1, for Ros = 1', enExit)
```
 
%% Cell type:code id:3a039428 tags:
 
``` python
# x[n] --> BPF --> D --> LOdown --> y[m D, k] [HARRIS Fig 6.7]
aBpf = aPrototype * np.exp(1j * w_k * np.arange(Ncoefs))
yBpfData = signal.lfilter(aBpf, [1.0], xData)
yBpfDown = down(yBpfData, Ndown)
yBpfDownLo = yBpfDown * LOdown # = y[m D, k]
 
# result is True for any Ndft, Ndown, because LOdown is in equation of yBpfDownLo
result = np.all(np.isclose(yDown, yBpfDownLo))
verify_result(result, ': yDown == yBpfDownLo, for Ros >= 1', enExit)
```
 
%% Output
 
PASSED: yDown == yBpfDownLo, for Ros >= 1
 
%% Cell type:markdown id:589b5e1e tags:
 
## 3.3 BPF and LO at downsampled rate: D --> poly BPF --> LOdown
 
Partition the BPF FIR filter H(z) in Ndown polyphases to have Hp(z^Ndown) per polyphase branch p, so that the down sampling can be done before the BPF by using the Noble identity.
 
%% Cell type:markdown id:6d43919d tags:
 
### 3.3.1 Maximally downsampled (= critically sampled)
 
%% Cell type:code id:327236c2 tags:
 
``` python
print('Ros =', Ros)
if Ndown == Ndft:
yMaxDownBpf = maximal_downsample_bpf(xData, Ndown, kLo, aPrototype)
yMaxDownBpfLo = yMaxDownBpf # = yMaxDownBpf * LOdown, because LOdown = 1 when Ndown == Ndft
 
result = np.all(np.isclose(yDown, yMaxDownBpfLo))
if not result:
plt.plot(m_sub, yDown.real, 'g.-')
plt.plot(m_sub, yDown.imag, 'g.--')
plt.plot(m_sub, yMaxDownBpfLo.real, 'r-')
plt.plot(m_sub, yMaxDownBpfLo.imag, 'r--')
verify_result(result, ': yDown == yMaxDownBpfLo, for Ros = 1', enExit)
```
 
%% Output
 
Ros = 4/3
 
%% Cell type:markdown id:c154ebe1 tags:
 
### 3.3.2 Non-maximally downsampled (= oversampled)
 
%% Cell type:code id:aefa8615 tags:
 
``` python
yDownBpf = non_maximal_downsample_bpf(xData, Ndown, kLo, Ndft, aPrototype)
yDownBpfLo = yDownBpf # = yDownBpf * LOdown * LOshift, because LOdown = 1 when Ndown == Ndft
# and LOshift phase shifts compensate for time shift due to Ndown < Ndft
 
result = np.all(np.isclose(yDown, yDownBpfLo))
if not result:
plt.plot(m_sub, yDown.real, 'g.-')
plt.plot(m_sub, yDown.imag, 'g.--')
plt.plot(m_sub, yDownBpfLo.real, 'r-')
plt.plot(m_sub, yDownBpfLo.imag, 'r--')
verify_result(result, ': yDown == yDownBpfLo for any Ros', enExit)
```
 
%% Output
 
> non_maximal_downsample_bpf():
. len(x) = 110592
. Nx = 110449
. Nblocks = 768
. len(yc) = 768
. Ros = 1.3333333333333333
. Ndown = 144
. Ndft = 192
. k = 2
PASSED: yDown == yDownBpfLo for any Ros
 
%% Cell type:code id:860a1f90 tags:
 
``` python
# Only check bin phase and bin ampl for CW at center of bin, and with averaging over
# least Ntaps number of subband periods
verifyBin = wgSubIsInt and not wgModulation and SNR_WG_dB >= 100 and Nsim >= 2 * Ntaps
```
 
%% Cell type:markdown id:06e66708 tags:
 
 
%% Cell type:code id:f814c9b9 tags:
 
``` python
# Verify subband bin phase
if verifyBin:
print('wgSub = %.1f, Nsim = %d [Tsub]' % (wgSub, Nsim))
# The phaseMargin >> c_atol, because it depends on the stop band attenuation of the
# aPrototype LPF. This is because the LO downconverts the positive frequency band
# of the WG cos() wave, so the negative frequency band will appear in the stop band.
# Use fixed approriate phase margin for carrier only, so with no AWGN.
phaseMargin = 1 # 0.01
binPhase = np.mean(np.angle(yDownBpfLo[Ntaps:], deg=True))
diffPhase = np.abs(wgPhase - binPhase)
snrPhase_dB = pow_db(1 / diffPhase)
print('binPhase = %.3f, diffPhase = %.3f, phaseMargin = %.3f [deg], snrPhase_dB = %.1f dB' %
(binPhase, diffPhase, phaseMargin, snrPhase_dB))
result = np.isclose(wgPhase, binPhase, atol=phaseMargin)
verify_result(result, ': wgPhase ~= binPhase', enExit)
```
 
%% Cell type:code id:dd7a9503 tags:
 
``` python
# Verify subband bin amplitude
if verifyBin:
print('wgSub = %.1f, Nsim = %d [Tsub]' % (wgSub, Nsim))
amplMargin = wgAmpl * 0.01
binAmpl = np.mean(np.abs(yDownBpfLo[Ntaps:]))
print('wgAmpl = %.3f, binAmpl = %f, amplMargin = %.3f' % (wgAmpl, binAmpl, amplMargin))
result = np.isclose(wgAmpl / nofSsb, binAmpl, atol=amplMargin)
verify_result(result, ': wgAmpl / nofSsb ~= binAmpl', enExit)
```
 
%% Cell type:markdown id:830631a1 tags:
 
# 4 Single channel upconverter
 
%% Cell type:markdown id:ec976ece tags:
 
# 4.1 Full rate: U --> LPF --> LOp
 
Upsample baseband, then LPF at baseband [HARRIS Fig 7.9] and upconvert to bin kLo. Make real bandpass signal by adding bin -kLo.
 
%% Cell type:code id:3cf6aa74 tags:
 
``` python
# Up mixer local oscillator (LO) for channel k
# . LOp = 1 / LO with p to indicate positive sign in LOp = np.exp(+1j * w_k * n_s)
# . Adjust LOp by hPairGroupDelay to be able to align reconstructed output data with input data
LOadjust = np.exp(-1j * w_k * hPairGroupDelay * Ts)
LOp = np.exp(+1j * w_k * n_s)
```
 
%% Cell type:code id:373ab5bb tags:
 
``` python
# Back to back, down converter up converter, reference output yrUpLpfLo
# With (inefficient) processing at full rate, because down sampling is done after LO and LPF and
# upsampling is done before LPF and LOp.
# yDown ycUp ycUpLpf
# y[mD, k] --> U ------> LPF --------> LOp --> ycUpLpfLo
ycUp = up(yDown, Nup) # insert Nup - 1 zeros
ycUpLpf = Nup * signal.lfilter(sPrototype, [1.0], ycUp) # interpolate by Nup
ycUpLpfLo = ycUpLpf * LOp * LOadjust # upconvert to positive bin kLo
yrUpLpfLo = ycUpLpfLo.real * nofSsb # = ycUpLpfLo + np.conj(ycUpLpfLo), add negative bin -kLo
```
 
%% Cell type:code id:a5c60be5 tags:
 
``` python
# Plot original real xData recovered yrUpLpfLo
xDelayed = xData[0:len(yrUpLpfLo) - hPairGroupDelay]
xyDiff = yrUpLpfLo[hPairGroupDelay:] - xDelayed
if 1:
plt.plot(xDelayed, 'r.-')
plt.plot(yrUpLpfLo[hPairGroupDelay:], 'b.-')
plt.plot(xyDiff, 'g')
else:
plt.plot(xyDiff, 'g')
plt.xlim([0, 1000])
 
if not wgModulation:
offset = Ncoefs
xLen = len(xDelayed[offset:])
if xLen < Ncoefs:
print('Too low Nsim for proper yrSNR calculation %d < %d' % (xLen, Ncoefs))
verify_result(False, '', enExit)
yrSNR = snr_db(xDelayed[offset:], xyDiff[offset:])
yrAmpl = np.sqrt(np.mean(np.abs(yrUpLpfLo[Ncoefs:]**2)) * nofSsb)
print('wgSub = %f' % wgSub)
print('kLo = %d' % kLo)
print('xLen = %d input samples' % xLen)
print('yrSNR = %.5f [dB], single bin' % yrSNR)
print('yrAmpl = %f' % yrAmpl)
```
 
%% Output
 
wgSub = 1.500000
kLo = 2
xLen = 105985 input samples
yrSNR = 5.37011 [dB], single bin
yrAmpl = 0.461804
 
 
%% Cell type:markdown id:18e8463e tags:
 
# 4.2 LOp at downsampled rate: LOpDown --> U --> BPF
 
Use LOp at downsampled rate. Do upconversion by U from baseband to kLo at sampled rate. Use BPF centered at kLo (is LPF shifted by +kLo) at sample rate.
 
If Nup = Ndft, then U w_k = U 2pi * k / Ndft is multiple of 2pi, so then LOpDown = 1.
 
%% Cell type:code id:7049249e tags:
 
``` python
# LOpDown --> U --> LOp
#
# LOpD = exp(-j * U w_k * m)
# = exp(-j * 2 pi U / Ndft * k * m)
LOpDown = down(LOp, Nup)
 
U_w_k = Nup * w_k
LOpD = np.exp(1j * U_w_k * m_i)
 
# Verify that LOp rotates with w_k and LOpDown with U * w_k rad/s
result = np.all(np.isclose(LOpDown, LOpD))
if not result:
plt.plot(m_sub, LOpDown.real, 'r-')
plt.plot(m_sub, LOpDown.imag, 'r--')
plt.plot(m_sub, LOpD.real, 'g-')
plt.plot(m_sub, LOpD.imag, 'g--')
verify_result(result, ': LOpDown == LOpD', enExit)
```
 
%% Output
 
PASSED: LOpDown == LOpD
 
%% Cell type:code id:64cc34f3 tags:
 
``` python
# Verify that LOpDown == 1 when Nup == Ndft
if Nup == Ndft:
result = np.all(np.isclose(LOpDown, 1.0))
if not result:
plt.plot(m_sub, LOpDown.real, 'r-')
plt.plot(m_sub, LOpDown.imag, 'r--')
verify_result(result, ': LOpDown == 1, for Ros = 1', enExit)
```
 
%% Cell type:code id:71a91beb tags:
 
``` python
# ycDownLo ycLoUp
# y[mD, k] --> LOpDown ---------> U -------> BPF --> ycLoUpBpf
sBpf = sPrototype * np.exp(1j * w_k * np.arange(Ncoefs))
ycDownLo = yDown * LOpDown * LOadjust # upconvert to positive bin kLo
ycLoUp = Nup * up(ycDownLo, Nup) # insert Nup - 1 zeros
ycLoUpBpf = signal.lfilter(sBpf, [1.0], ycLoUp) # interpolate by Nup with BPF at kLo
yrLoUpBpf = ycLoUpBpf.real * nofSsb # = ycLoUpBpf + np.conj(ycLoUpBpf), add negative bin -kLo
 
# result is True for any Ndft, Ndown, because LOdown is in equation of yBpfDownLo
result = np.all(np.isclose(yrUpLpfLo, yrLoUpBpf))
if not result:
plt.plot(yrUpLpfLo, 'r')
plt.plot(yrLoUpBpf, 'b')
#plt.xlim([0, 850])
verify_result(result, ': yrUpLpfLo == yrLoUpBpf, for Ros >= 1', enExit)
```
 
%% Output
 
PASSED: yrUpLpfLo == yrLoUpBpf, for Ros >= 1
 
%% Cell type:markdown id:2830cdc7 tags:
 
# 4.3 BPF and LOp at downsampled rate: LOpDown --> poly BPF --> U
 
Partition the BPF FIR filter H(z) in Nup polyphases to have Hp(z^Nup) per polyphase branch p, so that the upsampling can be done after the BPF by using the Noble identity.
 
%% Cell type:markdown id:28fced9c tags:
 
# 4.3.1 Maximally upsampled (= critically sampled)
 
%% Cell type:code id:a3aae48a tags:
 
``` python
print('Ndft =', Ndft)
if Nup == Ndft:
# ycDownLo = yDown * LOdown = yDown, because LOdown = 1 when Ndown == Ndft
ycLoBpfMaxUp = maximal_upsample_bpf(yDown, Nup, kLo, sPrototype)
yrLoBpfMaxUp = ycLoBpfMaxUp.real * nofSsb # add negative bin -kLo to make real
 
result = np.all(np.isclose(yrUpLpfLo, yrLoBpfMaxUp))
if not result:
plt.plot(n_sub, yrUpLpfLo, 'g.-')
plt.plot(n_sub, yrLoBpfMaxUp, 'r-')
plt.xlim([20, 30])
verify_result(result, ': yrUpLpfLo ==yrLoBpfMaxUp, for Ros = 1', enExit)
```
 
%% Output
 
Ndft = 192
 
%% Cell type:markdown id:825df86e tags:
 
# 4.3.2 Nonmaximally upsampled (= oversampled)
 
%% Cell type:code id:3c2c8ec5 tags:
 
``` python
if 0:
ycLoBpfUp = non_maximal_upsample_bpf(yDown, Nup, kLo, Ndft, sPrototype)
yrLoBpfUp = ycLoBpfUp.real * nofSsb # add negative bin -kLo to make real
 
result = np.all(np.isclose(yrUpLpfLo, yrLoBpfUp))
if not result:
plt.plot(n_sub, yrUpLpfLo, 'g.-')
plt.plot(n_sub, yrLoBpfUp, 'r.-')
plt.xlim([12,24])
verify_result(result, ': yrUpLpfLo == yrLoBpfUp, for Ros >= 1', enExit)
```
 
%% Cell type:markdown id:ee9daf9f tags:
 
# 5 Compare with DFT filterbank
 
Can use 'cw' or 'ccw' independently for analysis and synthesis PFB, because with fold() the IDFT can be expressed as a DFT and vice versa. However to have back to back DFT - IDFT in pipeline, use analysis 'cw' and synthesis 'ccw' [CROCHIERE 7.2.3].
 
%% Cell type:code id:4b55557c tags:
 
``` python
# Select analysis filterbank
aStructure = 'wola'
aStructure = 'pfs'
aCommutator = 'cw'
 
# Select synthesis filterbank
sStructure = 'wola'
sCommutator = 'ccw'
```
 
%% Cell type:code id:abb548f6 tags:
 
``` python
# Parameters
print('wgSub = %f' % wgSub)
print('kLo = %d' % kLo)
print('Ndft = %d' % Ndft)
print('Ndown = %d' % Ndown)
```
 
%% Output
 
wgSub = 1.500000
kLo = 2
Ndft = 192
Ndown = 144
 
%% Cell type:markdown id:d8bb436d tags:
 
# 5.1 Analysis
 
%% Cell type:code id:a6de2be1 tags:
 
``` python
# Verify that different analysis_dft_filterbank() structures yield exactly the same
print('Verify all bins:')
Ac0 = analysis_dft_filterbank(xData, Ndown, Ndft, aPrototype, 'pfs', commutator='cw', verbosity=vb)
Ac1 = analysis_dft_filterbank(xData, Ndown, Ndft, aPrototype, 'pfs', commutator='ccw', verbosity=vb)
Ac2 = analysis_dft_filterbank(xData, Ndown, Ndft, aPrototype, 'wola', verbosity=vb)
verify_result(np.all(np.isclose(Ac0, Ac1)), ': Ac0 == Ac1', enExit)
verify_result(np.all(np.isclose(Ac0, Ac2)), ': Ac0 == Ac2', enExit)
```
 
%% Output
 
Verify all bins:
PASSED: Ac0 == Ac1
PASSED: Ac0 == Ac2
 
%% Cell type:code id:a9ca3483 tags:
 
``` python
# Verify that analysis_dft_filterbank() bin kLo yields same as single channel reference
# . Only need to verify one Ac[kLo], when Ac0 = Ac1 = Ac2 = Ac3.
Ac = analysis_dft_filterbank(xData, Ndown, Ndft, aPrototype, aStructure, aCommutator)
yDownBin = Ac[kLo]
#yDownBin = yDownBin[:len(yDown)]
 
result = np.all(np.isclose(yDown, yDownBin))
if not result:
if 1:
plt.plot(m_sub, yDown.real, 'g.-')
plt.plot(m_sub, yDown.imag, 'g.--')
plt.plot(m_sub, yDownBin.real, 'r.-')
plt.plot(m_sub, yDownBin.imag, 'r.--')
else:
plt.plot(m_sub, yDown.real - yDownBin.real, 'b.-')
plt.plot(m_sub, yDown.imag - yDownBin.imag, 'b.--')
plt.xlim([0, 20])
verify_result(result, ': yDown == yDownBin', enExit)
```
 
%% Output
 
> analysis_dft_filterbank():
. len(x) = 110592
. Nblocks = 768
. Ros = 1.3333333333333333
. Ndown = 144
. Ndft = 192
. structure = pfs
. commutator = cw
PASSED: yDown == yDownBin
 
%% Cell type:markdown id:d3f49b8e tags:
 
# 5.2 Synthesis
 
%% Cell type:code id:caf25198 tags:
 
``` python
# Verify that different synthesis_dft_filterbank() structures yield exactly the same
Sr0 = synthesis_dft_filterbank(Ac, Nup, Ndft, sPrototype, 'wola', 'cw', verbosity=vb)
Sr1 = synthesis_dft_filterbank(Ac, Nup, Ndft, sPrototype, 'wola', 'ccw', verbosity=vb)
verify_result(np.all(np.isclose(Sr0, Sr1)), ': Sr0 == Sr1', enExit)
```
 
%% Output
 
PASSED: Sr0 == Sr1
 
%% Cell type:code id:47d5bf5b tags:
 
``` python
# Keep bin kLo and force all other bins to zero, to have exact comparison with full rate
# single channel reference.
Yc = Ac * 0
if kLo == 0:
# DC bin
Yc[kLo] = yDownBin
else:
Yc[kLo] = yDownBin
Yc[Ndft - kLo] = np.conjugate(yDownBin)
 
# Single bin synthesis from PFB
yr = synthesis_dft_filterbank(Yc, Nup, Ndft, sPrototype, sStructure, sCommutator)
yr = yr[0 : len(yrUpLpfLo)]
 
result = np.all(np.isclose(yrUpLpfLo, yr))
if not result:
plt.plot(n_sub, yrUpLpfLo, 'g.-')
plt.plot(n_sub, yr, 'r.-')
plt.xlim((20, 21))
verify_result(result, ': yrUpLpfLo == yr', enExit)
```
 
%% Output
 
> synthesis_dft_filterbank():
. Nblocks = 768
. Ros = 1.3333333333333333
. Nup = 144
. Ndft = 192
. structure = wola
. commutator = ccw
PASSED: yrUpLpfLo == yr
 
%% Cell type:code id:4818ad5a tags:
 
``` python
# Compare output with input
# . For Ros >= 1 the SNR for center bin wgSub approaches perfect reconstruction when Ntaps
# is increased, because then the accuracy of the interpolated samples improves.
# In the time domain this is understood because then there are more weighted values per
# polyphase. In the frequency domain this is understood by the improved bandstop
# suppression, because inserting zeros causes replicas that all need to be suppressed
# by the interpolating LPF.
# . For off-center wgSub the SNR reduces more due to that the gain of both LPF is then < 1.
 
# All bins synthesis from PFB
# . Use yr for comparison between single channel pipeline and using PFB for single bin
# . Use sr to see impact of aliasing in the other bins, due to limited stopband attenuation
# of both LPF.
sr = synthesis_dft_filterbank(Ac, Nup, Ndft, sPrototype, sStructure, sCommutator)
sr = sr[0 : len(yrUpLpfLo)]
 
# Output time aligned with input
x_yr = yr[hPairGroupDelay:]
x_sr = sr[hPairGroupDelay:]
 
# SNR
# . For single bin yr the SNR is sligthly better than for sr, when wgSub is at bin center, because
# then the LPF has gain 1 and the aliasing from the other bins is not in.
# . For all bin sr the SNR is much better than for yr, when wgSub is off bin center or with input
# noise is added, provided that the gain at the edge of the bin is sqrt(0.5).
# . The SNR for CW input calculated for refBunton and with construct.m are about equal within
# +-0.0005 dB. The difference is probably due to that the selected ranges for calculating the
# SNR are not exactly the same. Varying the offset shows a similar spread in SNR.
offset = Ncoefs
#offset = 10000
endset = 70000
xLen = len(xDelayed[offset:endset])
if xLen < Ncoefs:
print('Too low Nsim for proper SNR_yr calculation %d < %d' % (xLen, Ncoefs))
verify_result(False, '', enExit)
xDelayed_snapshot = xDelayed[offset:endset]
x_yr_snapshot = x_yr[offset:endset]
x_sr_snapshot = x_sr[offset:endset]
x_yr_diff_snapshot = xDelayed_snapshot - x_yr_snapshot
x_sr_diff_snapshot = xDelayed_snapshot - x_sr_snapshot
SNR_yr = snr_db(xDelayed_snapshot, x_yr_diff_snapshot)
SNR_sr = snr_db(xDelayed_snapshot, x_sr_diff_snapshot)
 
print('wgSub = %f' % wgSub)
print('kLo = %d' % kLo)
print('xLen = %d input samples' % xLen)
print('yrSNR = %.5f [dB], single bin' % yrSNR)
print('SNR_yr = %.5f [dB], single bin' % SNR_yr)
print('SNR_sr = %.5f [dB], all bins' % SNR_sr)
```
 
%% Output
 
> synthesis_dft_filterbank():
. Nblocks = 768
. Ros = 1.3333333333333333
. Nup = 144
. Ndft = 192
. structure = wola
. commutator = ccw
wgSub = 1.500000
kLo = 2
xLen = 67696 input samples
yrSNR = 5.37011 [dB], single bin
SNR_yr = 5.37011 [dB], single bin
SNR_sr = 22.18392 [dB], all bins
 
%% Cell type:code id:2c1ec80d tags:
 
``` python
# Plot original real xData recovered with sr for all bins
if 1:
plt.plot(xDelayed_snapshot, 'r.-')
plt.plot(x_sr_snapshot, 'g.-')
plt.plot(x_sr_diff_snapshot, 'y.-')
else:
plt.plot(x_sr_diff_snapshot, 'y.-')
plt.xlim([0, 250])
#plt.ylim([-0.01, 0.01])
```
 
%% Output
 
(0.0, 250.0)
 
 
%% Cell type:code id:3a647a98 tags:
 
``` python
```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment