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

Improve PFS definition

parent dd2d6f3b
Branches
No related tags found
1 merge request!419Resolve RTSD-265
Pipeline #91736 passed
......@@ -71,6 +71,7 @@
# Youtube: Guitars 4RL
# Youtube: MATHEMATICAL METHODS AND TECHNIQUES IN SIGNAL PROCESSING
# https://www.allaboutcircuits.com/technical-articles/pipelined-direct-form-fir-versus-the-transposed-structure/
# https://brianmcfee.net/dstbook-site/content/intro.html
#
# FBMC = filter bank-based multicarrier
......@@ -564,7 +565,7 @@ c) s-plane and z-plane
z^-1 --> b2 --> + + <-- -a2 <-- z^-1
. Direct-Form II: 1 / A(z) * B(z), to share z^-1 delay memories, but needs
more bits to avoid internal overflow.
more bits to avoid internal overflow. For FIR Direct-Form I = II.
x[n] --> + -------------v-----> b0 --> + --> y[n]
+ <-- -a1 <-- z^-1 --> b1 --> +
......@@ -868,7 +869,7 @@ c) s-plane and z-plane
. Linearity: a1 x1[n] + a2 x2[n] <==> a1 X1(w) + a2 X2(w)
. Scaling: x(t / a) <==> |a| X(a w)
. Time shift: x(t - T) <==> X(w) exp(-j w T)
x[n - k] <==> X(w) exp(-j w k), t = n Ts
x[n - m] <==> X(w) exp(-j w m), t = n Ts, T = m Ts
. Frequency shift (complex modulation): x[n] exp(+j v n) <==> X(w - v), is
dual of time shift
. Real modulation: x[n] cos(v n) <==> 1/2 [X(w + v) + X(w - v)]
......@@ -1179,6 +1180,40 @@ c) s-plane and z-plane
X_U(e^jw) = X(exp(jw U)), w in [0:2pi>
X(jwL) traverses unit circle U times
- Resampling:
. U --> LPF --> LPF --> D -->
<==>
U --> LPF --> D [LYONS Fig 10.9, CROCHIERE Fig 3.1]
. Resampling is upsampling with commutator yielding L / D samples per input
sample instead of L samples [LYONS 10.11].
- Time domain Down, Up, Resample U/D:
. down [CROCHIERE Eq 2.55]:
+inf +inf
y[m] = sum h[k] x[mD - k] = sum h[mD - n] x[n], using n = mD - k
k=-inf n=-inf
. up [CROCHIERE Eq 2.78, 2.80]:
+inf
y[m] = sum h[m - k] x[k / U], for k / U is integer, else x[] = 0 for
k=-inf inserted zeros
+inf
= sum h[m - rU] x[r] = v[m], using r = m // U - n yields:
r=-inf
+inf
= sum h[nU + m % U] x[m // U - n]
n=-inf
. resample [CROCHIERE Eq 2.88, 2.90]
+inf
y[m] = v[mD] = sum h[mD - rU] x[r], using r = mD // U - n yields:
r=-inf
+inf
= sum h[nU + mD % U] x[mD // U - n]
n=-inf
- Noble identities [LYONS Fig 10.20], HARRIS 2.2.1, VAIDYANATHAN Fig 4.2.3]
. Downsampling : x[n] --> H(z^Q) --> Q:1 --> y[m], is equivalent to:
x[n] --> Q:1 --> H(z) --> y[m]
......@@ -1207,88 +1242,93 @@ c) s-plane and z-plane
- Polyphase decomposition of H(z) [VAIDYANATHAN 4.3, PROAKIS 10.5.2,
CROCHIERE 3.3]:
. FIR structure:
. Direct-Form FIR is first apply delay z^(-q) then apply coefficient. Fits
downsampling because coefficients are then applied at low rate.
. Transposed Direct-Form FIR is first apply coefficient, then delay z^(-q)
result. Fits up sampling because coefficients are then applied at low rate.
result. Fits upsampling because coefficients are then applied at low
rate. The Direct-Form FIR would require fractional delays, which is not
realizable [CROCHIERE 3.3.1].
. The FIR sections in the polyphase branches can be implemented using any
form. The Transposed Direct-Form FIR can thus also be used for down
sampling, to make effcient use of z-1 delay line memory and multipliers
by loading the branch coefficients when they are needed and by passing on
the accumulated partial sums of the branches [HARRIS 5.5.1].
by loading the branch coefficients when they are needed and by passing
on the accumulated partial sums of the branches [HARRIS 5.5.1].
. Commutator model [VAIDYANATHAN 4.3.4, CROCHIERE 3.3.2, HARRIS Fig 5.4,
down 6.12, 6.13, up 7.11, 7.12]
down 6.12, 6.13, up 7.11, 7.12]:
. x[n] switch connects to n = 0 branch and the switch rotates in direction
of increasing n
. Assume FIR delay lines are drawn from top to bottom for phases q = 0 to
Q-1 for both Direct-Form and Transposed Direct-Form, then:
- Typically q = 0 is the connected branch [CROCHIERE], but [HARRIS]
connects q = Q-1 for the down sampler probably because it is the
connects q = Q-1 for the downsampler, probably because it is the
oldest so first phase in the FIR sum.
- downsampling input commutator rotates counter clockwise and yields
one output sample every rotation, because the summation stage is
combinatorial. The downsampled output depends on Q input samples from
the past. The output appears whenever the switch is back at branch
the past. The output appears whenever the switch is at branch
q = 0. A new FIR sum output starts being aggregated when switch is at
branch q = Q-1.
- up sampling output commutator rotates clockwise and yields U output
samples every rotation. The up sampled output creates Q output samples
in the future for every input sample. The first output appears when the
switch is at the connected branch q = 0 and then also a new input
arrives.
- upsampling output commutator rotates counter clockwise and yields U
output samples every rotation. The upsampled output creates Q output
samples in the future for every input sample. The first output appears
when the switch is at the connected branch q = 0 and then also a new
input arrives.
- The counter clockwise commutator for both down and up implies that a
type I polyphase representation is used, so with delays z^(-1) in the
delay line.
. Type I polyphase representation, based on delays z^(-q), yielding counter
clockwise commutator
(Can also notate Hq = Eq, hq = eq)
clockwise commutator.
H(z) = E0(z^Q) + E1(z^Q) z^-1 + E2(z^Q) z^-2 + ... + E_{Q-1}(z^Q) z^-(Q-1)
H(z) = H0(z^Q) + H1(z^Q) z^-1 + H2(z^Q) z^-2 + ... + H{Q-1}(z^Q) z^-(Q-1)
= E0(z^Q) +
z^-1 E1(z^Q) +
z^-2 E2(z^Q) +
= H0(z^Q) + # q = 0
z^-1 H1(z^Q) + # q = 1
z^-2 H2(z^Q) + # q = 2
... +
z^-(Q-1) E_{Q-1}(z^Q)
z^-(Q-1) H{Q-1}(z^Q) # q = Q-1
Q-1
= sum z^(-q) Eq(z^Q) [VAIDYANATHAN Eq 4.3.7, HARRIS Eq 6.6]
= sum z^(-q) Hq(z^Q) [VAIDYANATHAN Eq 4.3.7, HARRIS Eq 6.6]
q=0
where Eq(z^Q) is the z-transform of eq[n]:
where Hq(z^Q) is the z-transform of hq[n]:
eq[n] = h(Qn + q), +q for counter clockwise with delays z^(-1)
hq[n] = h(nQ + q), +q for counter clockwise with delays z^(-1)
[VAIDYANATHAN Eq 4.3.8]
+inf
Eq(z) = sum eq[n] z^(-n), 0 <= q <= Q - 1
Hq(z) = sum hq[n] z^(-n), 0 <= q <= Q - 1
n=-inf
- Note: For Q = 1 the Eq(z) are the single FIR coefficients in b
- Note: For Q = 1 the Hq(z) are the single FIR coefficients in b
. Type II polyphase, based on advances z^(q), yielding clockwise commutator:
. Type II polyphase, based on advances z^q, yielding clockwise commutator:
H(z) = z^(-(Q-1)) z^(Q-1) H(z), multiply by delay * advance = 1
= z^(-(Q-1)) * [z^(Q-1) E0(z^Q) + # q = 0
z^(Q-1) z^-1 E1(z^Q) + # q = 1
z^(Q-1) z^-2 E2(z^Q) + # q = 2
= z^(-(Q-1)) * [z^(Q-1) H0(z^Q) + # q = 0
z^(Q-1) z^-1 H1(z^Q) + # q = 1
z^(Q-1) z^-2 H2(z^Q) + # q = 2
... +
z^(Q-1) z^-(Q-2) E_{Q-2}(z^Q)] # q = Q-2
z^(Q-1) z^-(Q-1) E_{Q-1}(z^Q)] # q = Q-1
z^(Q-1) z^-(Q-2) H{Q-2}(z^Q) + # q = Q-2
z^(Q-1) z^-(Q-1) H{Q-1}(z^Q)] # q = Q-1
= z^(-(Q-1)) * [z^(Q-1) E0(z^Q) + # q = 0
z^(Q-1 - 1) E1(z^Q) + # q = 1
z^(Q-1 - 2) E2(z^Q) + # q = 2
= z^(-(Q-1)) * [z^(Q-1) H0(z^Q) + # q = 0 --> p = Q-1
z^(Q-2) H1(z^Q) + # q = 1 --> p = Q-2
z^(Q-3) H2(z^Q) + # q = 2 --> p = Q-3
... +
z^( 1) E_{Q-2}(z^Q) + # q = Q-2
E_{Q-1}(z^Q)] # q = Q-1
z^1 H{Q-2}(z^Q) + # q = Q-2 --> p = 1
H{Q-1}(z^Q)] # q = Q-1 --> p = 0
= z^(-(Q-1)) * [ E_{Q-1}(z^Q) + # q = Q-1 --> p = 0
z^( 1) E_{Q-2}(z^Q) + # q = Q-2 --> p = 1
= z^(-(Q-1)) * [z^(Q-1) R{Q-1}(z^Q) + # p = Q-1
z^(Q-2) R{Q-2}(z^Q) + # p = Q-2
z^(Q-3) R{Q-3}(z^Q) + # p = Q-3
... +
z^(Q-1 - 2) E2(z^Q) + # q = 2 --> p = Q-3
z^(Q-1 - 1) E1(z^Q) + # q = 1 --> p = Q-2
z^(Q-1) E0(z^Q) + # q = 0 --> p = Q-1
z^1 R1(z^Q) + # p = 1
R0(z^Q)] # p = 0
Q-1
H(z) = sum z^(-(Q-1-p)) Rp(z^Q) [VAIDYANATHAN Eq 4.3.9]
......@@ -1299,20 +1339,27 @@ c) s-plane and z-plane
p=0
where:
rp[n] = h(Qn - p), -p for clockwise with advances z^(+1)
p = Q - 1 - q
Rp(z) = E_{Q-1-p}(z), so flipud Eq phases, but keep coefficient order
per phase
rp[n] = h(nQ - p), -p for clockwise with advances z^(+1)
??? Is Type II polyphase same as Transposed Direct Form ???, because both still
implement H(z), they only yield different implementation structures. For Type 1
and Direct-Form the delay line goes from top to bottom in drawings. For Type 2
and Transposed we also want to draw the delay line from top to bottom (Fig 4.3.4,
Fig 4.3.11, Fig 4.5.6), but then we need to flipud the Eq phases.
. Phase q of h(n) with Q - 1 zeros, so [VAIDYANATHAN Fig 4.3.1]:
Rp(z) = H{Q-1-p}(z), so flipud Hq phases, but keep coefficient order
per phase
h[n] --> z^(-q) --> Q:1 --> eq[n]
. Type I polyphase does not imply Direct Form FIR and type II polyphase is
not the same as Transposed Direct Form FIR [CROCHIERE 3.3.3]:
* Type I uses delay line z^(-1) and +q in hq[n] = h(nQ + q) and yields
counter clockwise commutator for both downsample and upsample structure,
because branches are drawn such that q increases from top to bottom.
For downsampling q starts at q = Q-1, so at bottom to sum Q filtered
samples from the past, so input commutator rotates counter clockwise.
For upsampling q starts at q = 0 at top to produce Q samples until the
next input sample, so output commutator also rotates counter clockwise
because it is located at the output.
* Type II means using an advances line z^(+1) and -p in rp[n] = h[nQ - p]
and yields clockwise commutator, because p increases from bottom to
top. The advances can be implemented by allowing one sample delay at the
low rate.
- Variable bandwidth PFB [HARRIS, FARZAD, CHEN]
. Analysis --> select channels --> synthesis
......@@ -1361,6 +1408,7 @@ c) s-plane and z-plane
noise power at higher frequences will be filtered by the analoge LPF that
filters the DAC output.
13) Single channel down converter [HARRIS 6]
- Analogue I-Q downconverter, yields baseband signal:
xb[n] = exp(-j w_k n) x[n] = I[n] + j Q[n]
......@@ -1368,22 +1416,29 @@ c) s-plane and z-plane
. Q[n] = sin(-w_k n) x[n]
Performs complex multiply for real input, so I is the real result and Q is
the imaginary result.
- Mixer and LPF:
- Mixer local oscillator (LO) and LPF and downsampler D [HARRIS Fig 6.2]:
y[n, k] = exp(-j w_k n) x[n] * h[n] # * is convolution
for channel k converts positive w_k to baseband (because of -j),
w_k = 2pi k / M
D is downsample factor --> y[nD, k]
. LO(n, k) = exp(-j w_k n)
. M is number of frequency bins
. D <= M is downsample factor --> y[mD, k]
* First choose D = M for maximally downsampled (= critically sampled).
* Downsample the LO, using equivalency theorem.
. Down-conversion followed by a LPF is equivalent to a BPF followed by a
down-conversion:
down-conversion [HARRIS Fig. 6.5]:
exp(-j w_k n) --> H_LPF(z) = H_BPF(z) --> exp(-j w_k n)
exp(-j w_k n) --> H_LPF(z) <==> H_BPF(z) --> exp(-j w_k n)
where:
h_npf[n] = h_lpf[n] exp(j w_k n) # mix LPF up to BPF [HARRIS Fig. 6.6]
+inf
H_BPF(z) = sum h_lpf[n] exp(j w_k n) z^-n # mix LPF up to BPF
H_BPF(z) = sum h_lpf[n] exp(j w_k n) z^-n
n=-inf
+inf
......@@ -1392,43 +1447,103 @@ c) s-plane and z-plane
= H_LPF(z exp(-j w_k))
. The rotation rate of the sampled complex sinusoid is w_k radians per sample
at the input and D w_k radians per sample at the output, of the D:1 down
sampler:
. The rotation rate of the sampled complex sinusoid is w_k radians per
LO sample at the input and D w_k radians per LO sample at the output, of
the D:1 downsampler [HARRIS Fig. 6.7].:
H_BPF(z exp(-j w_k)) exp(-j w_k n) --> D
H_BPF(z exp(-j w_k)) --> D --> exp(-j D w_k n)
. The change in down converter rate is due to aliasing of the down sampling.
Now choose channel center frequencies w_k such that they will alias to DC
(zero frequency) as a result of the down sampling to D w_k. This occurs
when D w_k = k 2pi, so when w_k = k 2pi / D and then exp(-j D w_k n) = 1.
H_LPF(z exp(-j w_k)) exp(-j w_k n) --> D
<==>
H_LPF(z exp(-j w_k)) --> D --> exp(-j D w_k n)
H_BPF(z exp(-j w_k)) --> D, when M = D
. The change in LO down converter rate is due to aliasing of the
downsampling. Now choose channel center frequencies w_k such that they
will alias to DC (zero frequency) as a result of the downsampling to
D w_k. This occurs when D w_k = k 2pi, so when w_k = k 2pi / D and then
the downsampled LO term exp(-j D w_k n) = 1 vanishes when D = M
[HARRIS Fig. 6.8].
Hence the channel center frequencies have to be at integer multiples of the
output sample rate fs / D, so that they alias to baseband by the sample
rate change.
<==>
H_LPF(z exp(-j w_k)) --> D, when M = D
Hence the channel center frequencies have to be at integer multiples of
the output sample rate fs / M, so that they alias to baseband by the
sample rate change.
* Downsample the BPF, using Noble identity
. H(z^M) --> D <==> D --> H(z)
. The type I (for downsampling) polyphase FIR representation has a delay
line of z^-1 and Hp(z^M) in each branch. With this structure the
branch terms of H_BPF(z) = H_LPF(z exp(-j w_k)) reduce to Hp_BPF(z^M)
= Hp_LPF(z^M), because (z exp(-j w_k))^M = z^M, when D = M. With the
Noble identity this yields Hp_LPF(z^M) --> M ==> M --> Hp_LPF(z), so the
coefficients of the prototype filter.
Replacing z by z exp(-j w_k) the input delay line of z^-p becomes
z^-p exp(j w_k p). The constant phasors exp(j w_k p) can be placed
anywhere along each branch, so choose to place them just before the
summation [HARRIS Fig. 6.14].
. What remains is a downsampling polyphase FIR structure with the
coefficients of the prototype H_LPF, and a phasor in each branch. The
sum of the phasors for one k is the same as the DFT result for bin k.
- Type I polyphase representation of H(z), based on delays z^(-q), yielding
counter clockwise commutator
. H(z) = H0(z^M) + H1(z^M) z^-1 + H2(z^M) z^-2 + ... + H_{M-1}(z^M) z^-(M-1)
<==>
z^-p --> M --> Hp_LPF(z) * exp(-j w_k p) --> sum all p --> y[mM, k]
= H0(z^M) +
z^-1 H1(z^M) +
z^-2 H2(z^M) +
... +
z^-(M-1) H_{M-1}(z^M)
and y[mM, k] is [HARRIS Eq. 6.9 = 9.1]:
M-1
= sum z^(-q) Hq(z^M) [VAIDYANATHAN Eq 4.3.7, HARRIS Eq 6.6]
q=0
y[mM, k] = sum yp(mM) exp(j w_k p), with w_k = 2pi k / M
p=0
* Implementation can be shift in a sample and calculate one branch, or shift
in a block of D samples and calculate all branches. Hence the block size
is the downsample factor D.
* Now choose D < M for non-maximally downsampling (= oversampled), but keep
the polyphase FIR structure for M, in order to preserve using the
prototype H_LPF FIR coefficients and having the frequency grid for M bins.
. Shift in D samples and output sum of all branches after every D input
samples.
. A time shift of D samples causes a (bin) frequency dependent phase shift
of:
theta(w) = (D Ts) * (w_k / Ts)
= D * w_k, with w_k = 2pi k / M
For D = M this term vanishes, but for D < M it appears in y at the
output, similar as with the DFT shift property:
Time shift: DFT(x[n - D]) <==> DFT(x) exp(-j w_k D).
. Apply Noble identity:
. Therefore for every new output sample needs to be counter rotated by:
exp(-j w_k D m) = exp(-j 2pi k D / M m)
. This can also be interpreted as to calculate bin k for downsampled output
samples m = 0, 1, ..., calculate the output y for bin k + (mD) % M.
Instead of first caluclating the output y for bin k and then multiplying
y by exp(-j 2pi k D / M m).
exp(j w_k) exp(j w_k D m) = exp(j w_k(1 + mD))
14) Polyphase filterbank (PFB) [HARRIS Fig 6.21, 9.21]
. The PFB implements M single channel down converters to output all k bins.
The output rate per branch p in the single channel down converter is a
factor M less, so all k = 0:M-1 bins of y[mM, k] can be calculated by using
the IDFT (because +j in phasor) to sum the p = 0:M-1 branches for each k
= 0:M-1 bins, because DFT is:
N-1
X(k) = sum x(n) exp(-j w_k n)
n=0
Hence the PFB can be interpreted as M independent single channel down
converters in parallel, that share the same downsampling polyphase FIR
structure, with the coefficients of the H_LPF prototype filter, and use the
IDFT to calculate all M bins. The PFB output has a commutator that starts
at bin 0, so bin k = M-1 is the the newest and output last.
14) Polyphase filterbank [HARRIS 6]
15) Quadrature Mirror Filter (QMF) [CROCHIERE 7.7, PROAKIS 10.9.6]
......@@ -1553,9 +1668,10 @@ Appendix B: Teaser talk: Oversampled polyphase filterbank (OPFB)
- Goal of multirate processing is to process at the rate just sufficient for
the signal bandwidth, not higher
4 Comparison with sampling theorem of analogue signal
- The sampling theorem [CROCHIERE Fig 2.2, 2.3]
- LPF selects the baseband signal, BPF can select any band at multiple of
subband rate fs / Q
- Equivalence criterium mixer + LPF + D = BPF + D
subband rate fdown = fs / Q
- Equivalence criterium LO + LPF <==> BPF + LO [HARRIS 6]
5 Single channel mixer
6 Filterbank for D = U = Q channels --> DFT section
7 Oversampled filterbank
......
%% Cell type:markdown id:c69a2eb8 tags:
# Multirate mixer
Author: Eric Kooistra, May 2024
Purpose:
* Practise DSP [1].
* Use multirate processing to implement a mixer
References:
1. dsp_study_erko, summary of DSP books
2. chapter 6 in [HARRIS]
%% Cell type:code id:02689e50 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from fractions import Fraction
```
%% Cell type:code id:65235f50 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.firfilter import filterbank_frequency_response
from rtdsp.fourier import dtft
from rtdsp.multirate import down, up, maximal_downsample_bpf
from rtdsp.plotting import plot_power_spectrum, plot_magnitude_spectrum
```
%% Output
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
%% Cell type:code id:c49515de tags:
``` python
# Filterbank
Ntaps = 8 # number of taps per poly phase FIR filter
Ndft = 16 # DFT size
Ncoefs = Ndft * Ntaps
```
%% Cell type:code id:c5c90a6b tags:
``` python
# Samples
fs = 1.0 # sample rate
Ts = 1 / fs # sample period
```
%% Cell type:code id:6d3a14bc tags:
``` python
# Subbands
Nsub = Ndft // 2 # number of subbands in fs / 2
fsub = fs / Ndft # subband frequency
Tsub = 1 / fsub # subband period
```
%% Cell type:markdown id:15bf6804 tags:
# 1 Prototype FIR low pass filter
%% Cell type:code id:9aa3a1ae tags:
``` python
# 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
# firwin hanning filter
hpFactor = 1.1
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)
```
%% Cell type:code id:0a69b385 tags:
``` python
# Plot impulse response
plt.figure(1)
plt.plot(hPrototype, 'g')
plt.title('Impulse response')
plt.grid(True)
```
%% Output
%% Cell type:markdown id:e3542672 tags:
%% Cell type:code id:0b979a1f tags:
``` python
# Filterbank aggregate frequency response
h, f, HFprototype = dtft(hPrototype)
HFbank = filterbank_frequency_response(HFprototype, 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)
# Plot transfer function (use frequency in fsub units)
fLim = None
fLim = (0, 2)
dbLim = None
dbLim = (-100, 5)
voltLim = None
plt.figure(1)
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)
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
```
%% Output
%% Cell type:markdown id:e130f714 tags:
# 2 Generate input data
%% Cell type:code id:d76e42f5 tags:
``` python
# Time
Nsim = 10 # number of subband periods to simulate
Nsamples = Nsim * Ndft
# . 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]
```
%% Cell type:code id:48d4a3b3 tags:
``` python
# Carriers
# . freq = center subband yields constant baseband I, Q signal
# . phase = 0 yields Q = 0
subbands = np.array([1.0]) # in range(Nsub)
freqs = subbands * fsub # in Hz
phases = [0.0] # in degrees
ampl = 1
xData = np.zeros(Nsamples) # = x[n]
for freq, phase in zip(freqs, phases):
xData += ampl * np.cos(2 * np.pi * freq * n_s + np.radians(phase))
```
%% Cell type:code id:0ae28649 tags:
``` python
# Additive Gaussian White Noise (AGWN)
SNR_dB = 100 # signal to noise ratio between one carrier and the noise
seed = None
seed = 1
rng = np.random.default_rng(seed)
mu = 0.0
sigma = ampl * np.sqrt(0.5) / 10**(SNR_dB / 20)
noise = rng.normal(mu, sigma, Nsamples)
xData += noise
# Check SNR, each extra carrie adds 3 dB
powCarriers = np.sum(xData**2) / Nsamples
powNoise = np.sum(noise**2) / Nsamples
snr_db = 10 * np.log10(powCarriers / powNoise)
print('powCarriers = %f' % powCarriers)
print('powNoise = %f' % powNoise)
print('SNR = %f dB' % snr_db)
```
%% Output
powCarriers = 0.500000
powNoise = 0.000000
SNR = 100.959873 dB
%% Cell type:markdown id:06e69ebc tags:
# 3 Single channel downconverter
%% Cell type:code id:adc33e70 tags:
``` python
# Downsample rate
Ndown = Ndft
#Ndown = Ndft // 2
fdown = fs / Ndown # downsampled data rate
Tdown = 1 / fdown # downsampled data period
# Time
Msamples = Nsamples // Ndown
# . 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]
```
%% Cell type:markdown id:127b23a8 tags:
# 3.1 Full rate: LO --> LPF --> D
Down convert bin kLo to baseband, then LPF still at sample rate and then downsample [HARRIS Fig 6.2].
%% Cell type:code id:9d3fb1c8 tags:
``` python
# Mixer local oscillator (LO) for channel k
kLo = np.fix(np.round(subbands[0]))
w_k = 2 * np.pi * kLo / Ndft
LO = np.exp(-1j * w_k * n_s)
```
%% Cell type:code id:50334d52 tags:
``` python
# x[n] --> LO --> LPF --> D --> y[mD, k] [HARRIS Fig 6.2]
xLoData = xData * LO
yData = signal.lfilter(hPrototype, [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 0x7f202f761760>
<matplotlib.legend.Legend at 0x7f48da4b6b80>
%% Cell type:markdown id:06436339 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 by 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:bd3cf3b3 tags:
``` python
x = np.array([0,1,2,3,4,5,6,7,8,9,10,11])
x.reshape((4,3))
```
%% Output
array([[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8],
[ 9, 10, 11]])
%% Cell type:code id:25317d4e 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
if np.all(np.isclose(LOdown, loD)):
print('PASSED')
else:
print('FAILED')
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--')
```
%% Output
Ndft = 16
Ndown = 16
w_k = 0.39269908169872414
D_w_k = 6.283185307179586
PASSED
%% Cell type:code id:6f50284f tags:
``` python
# Verify that LOdown == 1 when Ndown == Ndft
if Ndown == Ndft:
if np.all(np.isclose(LOdown, 1.0)):
print('PASSED')
else:
print('FAILED')
else:
plt.plot(m_sub, LOdown.real, 'r-')
plt.plot(m_sub, LOdown.imag, 'r--')
```
%% Output
PASSED
%% Cell type:code id:3ae19d32 tags:
``` python
# x[n] --> BPF --> D --> LOdown --> y[m D, k] [HARRIS Fig 6.7]
hBpf = hPrototype * np.exp(1j * w_k * np.arange(Ncoefs))
yBpfData = signal.lfilter(hBpf, [1.0], xData)
yBpfDown = down(yBpfData, Ndown)
yBpfDownLo = yBpfDown * LOdown # = y[m D, k]
if np.all(np.isclose(yDown, yBpfDownLo)):
# True for any Ndft, Ndown, because LOdown is in equation of yBpfDownLo
print('PASSED')
else:
print('FAILED')
```
%% Output
PASSED
%% Cell type:markdown id:f0642906 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:fe3e5fb8 tags:
### 3.3.1 Maximally downsampled (= critically sampled)
%% Cell type:code id:7a8f3e37 tags:
``` python
print('Ndft =', Ndft)
if Ndown == Ndft:
yDownBpf = maximal_downsample_bpf(xData, Ndown, kLo, hPrototype)
yDownBpfLo = yDownBpf # = yDownBpf * LOdown, because LOdown = 1 when Ndown == Ndft
if np.all(np.isclose(yDown, yDownBpfLo)):
print('PASSED')
else:
print('FAILED')
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--')
```
%% Output
Ndft = 16
> downsample_bpf():
> Log downsample_bpf():
. len(x) = 160
. Ndown = 16
. Nx = 145
. Nxp = 10
. len(y) = 10
. k = 1.0
PASSED
%% Cell type:markdown id:73a63d82 tags:
### 3.3.2 Fractionally downsampled (= oversampled)
### 3.3.2 Non-maximally downsampled (= oversampled)
%% Cell type:code id:6dfb2975 tags:
``` python
```
......
Source diff could not be displayed: it is too large. Options to address this: view the blob.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment