Skip to content
Snippets Groups Projects

Resolve RTSD-265

Merged Resolve RTSD-265
RTSD-265 into master
Merged Eric Kooistra requested to merge RTSD-265 into master
Viewing commit 6d9fdfee
Prev Next
Show latest version
1 file
+ 184
− 33
Compare changes
  • Side-by-side
  • Inline
applications/lofar2/model/pfb_os/dsp_study_erko.txt
+ 184
− 33
@@ -190,7 +190,7 @@ b) z-transform [LYONS 6.3, MATLAB]
z-transform on the unit circle.
- Properties [PROAKIS 3.3]:
. shift property: x[n - k] <--> X(z) z^-k
. time reversal: x[-n] <--> X(z^-1]
. time reversal: x[-n] <--> X(z^-1)
. convolution: x[n] * y[n] <--> X(z) Y(z)
. conjugation: x*[n] <--> X*(z*)
. real part: Re{x[n]} <--> 1/2 [X(z) + X*(z*)]
@@ -207,7 +207,7 @@ b) z-transform [LYONS 6.3, MATLAB]
See section on IIR for inverse z-transform that shows that the b, a
coefficients of H(z) in the z-domain directly yield the same coefficients
for the difference equation in the time domain).
for the difference equation in the time domain.
c) s-plane and z-plane
- The z-plane is not only sampling of s-domain, but also maping from
@@ -236,7 +236,7 @@ c) s-plane and z-plane
- Properties of rectangular window with M points from [JOS4 3.1.2]:
. M determines transition band width, window shape determines stop band
side lobe level [HARRIS 8.2]
. Zero crossings at integer multiples of 2 pi / M = Ndtft / M [LYONS Eq.
. Zero crossings at integer multiples of 2pi / M = Ndtft / M [LYONS Eq.
3.45]
. Main lobe width is 4 pi / M between zero crossings.
. Main lobe cutoff frequency is about fpass ~= 0.6 / M, where fpass is
@@ -265,17 +265,27 @@ c) s-plane and z-plane
. h_ideal(n) = 2 f_cutoff sinc(2 f_cutoff n), n in Z
- f_cutoff is normalized cutoff frequency with fs = 1
- f_cutoff is positive frequencies -6 dB power BW
- sinc(n / N) is NYQUIST(N) and has bandwidth 2pi / N [JOS4 9.5, 9.7]
. cutoff frequency remarks:
- [Scipy firwin] defines f_c relative to fNyquist = fs / 2 instead of fs,
so need to specify f_c = f_cutoff / 2.
- [Scipy firls, remez] define fpass relative to fs, so specify fc =
f_cutoff.
- Nyquist filter = NYQUIST(L)
. A Nyquist filter has 0.5 gain, so -3 dB gain, cutoff frequency at
pi / L, or fNyquist / L, and are also known as L-th band filters.
. The impulse response is zero at every rL-th sample for abs(r) = 1, 2,...
. With L samples per symbol this yields zero intersymbol interference (ISI)
. For example: windowed sinc() (= Portnoff window), Scipy firwin() uses
windowed sinc().
. For sinc() the ideal bandwidth is 2pi / L
- LPF FIR filter design [LYONS 5.3]
. Methods based on desired response characteristics [MNE]:
- Frequency-domain design (construct filter in Fourier domain and use an
IFFT to invert it, MATLAB fir2)
- Windowed FIR design (scipy.signal.firwin(), firwin2(), and MATLAB fir1
with default Hamming).
IFFT to invert it, MATLAB fir2, scipy.signal.firwin2() uses window after
IFFT)
- Windowed FIR design (scipy.signal.firwin(), and MATLAB fir1 with default
Hamming).
- Least squares designs (scipy.signal.firls(), MATLAB firls, fircls1)
. firls = least squares
. fircls, fircls1 = constrained ls with pass, stop ripple
@@ -297,21 +307,16 @@ c) s-plane and z-plane
. Note coefs = flip(h) is important for anti-symmetric
- Band Pass Filter (BPF):
. h_bp[k] = h_lp[k] * s_shift[k]
= h_lp[k] * cos(2 pi fc k)
= h_lp[k] * cos(2pi fc k)
= h_lp[k] * cos(k pi / 2), for half band fc = fs / 4,
series [1, 0, -1, 0]
- High Pass Filter (HPF):
. h_hp[k] = h_lp[k] * s_shift[k]
= h_lp[k] * cos(2 pi fc k)
= h_lp[k] * cos(2pi fc k)
= h_lp[k] * cos(k pi), for fc = fs / 2,
series [1, -1]
- Nyquist filter
. A Nyquist filter has 0.5 gain, so -3 dB gain, cutoff frequency at
pi / L, or fNyquist / L, and are also known as L-th band filters.
The impulse response is zero every L-th sample.
5) Finite Impulse Response (FIR) filters
- FIR filters perform time domain Convolution by summing products of shifted
@@ -500,7 +505,7 @@ c) s-plane and z-plane
xr --> complex BPF --> xa = xI + j xQ, where xQ = HT(xI)
- creates pair of in-phase (I) and quadrature (Q) channels
- h_lp[k] is Ntaps LPF with two sided BW equal to BPF BW
- h_bp[k] = h_lp[k] * exp(2 pi fcenter / fs (k - D)), mixer
- h_bp[k] = h_lp[k] * exp(2pi fcenter / fs (k - D)), mixer
= h_cos[k] + j h_sin[k]
. k = 0 ... Ntaps - 1
. D = (Ntaps - 1) / 2, but same for both I and Q, so even Ntaps is as
@@ -653,7 +658,7 @@ c) s-plane and z-plane
1 + z^-1
w_d = 2 arctan(w_a Ts / 2), with -pi <= w_d <= +pi,
f_d = w_d / (2 pi) fs
f_d = w_d / (2pi) fs
. scipy.signal.bilinear(b, a, fs)
scipy.signal.bilinear(z, p, k, fs)
@@ -712,7 +717,7 @@ c) s-plane and z-plane
The phase shift is now -pi at the fb break or center frequency and controlled
by parameter d:
d = -cos(2 pi fb / fs)
d = -cos(2pi fb / fs)
The parameter c determines the bandwidth (BW) or steepness of the phase slope
through fb:
@@ -870,7 +875,7 @@ c) s-plane and z-plane
. Conjugation: x*[n] <==> X*(-w)
. Convolution:
x * y <==> X Y
x y <==> 1 / (2 pi) X * Y
x y <==> 1 / (2pi) X * Y
. flip(x) <==> flip(X), so when signal is folded (time reversed) about the
origin in time, then its magnitude spectrum remains unchanged, and the
phase spectrum changes sign (phase reversal).
@@ -882,7 +887,7 @@ c) s-plane and z-plane
k=-inf
<==>
+inf
d_train_P(w) = 2 pi / P sum d(w - 2 pi k / P)
d_train_P(w) = 2pi / P sum d(w - 2pi k / P)
k=-inf
. Sampling: x_d(t) = x_a(t) d_train_Ts(t)
@@ -925,6 +930,8 @@ c) s-plane and z-plane
resolution in time, R <= W. R = 1 is 'sliding FFT'. R > 1 is downsampling
in time. N > W and R < W cause that the output sample rate increases.
- NYQUIST(N) [JOS4 9.5]
A window w is NYQUIST(N) when w(rN) = 0 for abs(r) = 1, 2, ...
- OLA = Overlap Add and uses rectangular window w of length Lx on x and step
with Lx through x. The purpose of OLA is to use FFT and implement
convolution with impulse response h as multiplication of X(f) H(f) in
@@ -941,6 +948,8 @@ c) s-plane and z-plane
- Any window is COLA(1).
- Retangular window is COLA(W), but also COLA(W / k).
- Hamming window is COLA(W / 2) and COLA(W / 4)
- Portnoff window: w[n] = w0[n] sinc(n / N) [JOS4 9.5, 9.7]
. is COLA(2pi / N) and is NYQUIST(N)
. With COLA the sum X_m(w) = X(w) = DTFT_w(x)
. With filtering, then STFT is used like OLA to perform spectral
modifications via H(f).
@@ -948,9 +957,9 @@ c) s-plane and z-plane
OLA) and one after IFFT. Typically both windows are equal and the sqrt(w),
so that together they are COLA [JOS4 8.6].
. COLA contraint in frequency domain [JOS4 8.3.2]. Window transform is zero
for all harmonics of the frame rate 2 pi / R:
for all harmonics of the frame rate 2pi / R:
W(w_k) = 0, with w_k = 2 pi k / R and for |k| = 1,2,...,R-1
W(w_k) = 0, with w_k = 2pi k / R and for |k| = 1,2,...,R-1
This is the weak COLA contraint. The strong COLA constraint, that is
useful if the signal is modified in the frequency domain, requires that
@@ -1032,6 +1041,9 @@ c) s-plane and z-plane
This property can be used to advantage when dealing with bandpass signals by
associating the bandpass signal with one of these images instead of with the
baseband [CROCHIERE 2.4.2].
??? Sampling werkt dus al als demodulatie (mixen naar baseband). Igv Ros = 1
past D, U = Nfft, zie Harris. Igv Ros > 1 ontstaat er een extra menging, die
gecompenseerd moet worden. Hoe past het allemaal bij elkaar ???
- Polyphase filtering ensures that only the values that remain are calculated,
so there are D or U phases [LYONS 10.7]. The LPF with all phases is called
the prototype filter. Do not calculate samples that will be:
@@ -1042,7 +1054,7 @@ c) s-plane and z-plane
- Sampling, downsampling and upsampling
. Sampling causes the analoge spectrum to alias around k 2pi, similar for
downsamping the the digital spectrum aliases around k 2pi / D, as if the
downsamping the digital spectrum aliases around k 2pi / D, as if the
analogue signal was sampled directly at the downsampled rate [LYONS 10.1].
. Downsampled spectrum [LYONS 10.3.2]
1. Draw original spectrum beyond -2pi to + 2pi, to show 0 and at least one
@@ -1089,7 +1101,9 @@ c) s-plane and z-plane
- Downsampling [LYONS 10.1, PROAKIS 10, CROCHIERE Fig 3.2, VAIDYANATHAN Fig
4.1.4, JOS4 11.1, SP4COMM 11.1.2]:
x_D[n] = x[nD], because x_D removes D-1 values from x
Index n for up rate, index m for down rate.
x_D[m] = x[mD], because x_D removes D-1 values from x
Define x' at x time grid, but with x' = 0 for samples in x that will be
discarded. This operation has no name in literature, probably because it
@@ -1188,7 +1202,7 @@ c) s-plane and z-plane
z^(-k) = z^(-k(d*D - u*U))
= z^(-kdD) * z^(kuU)
which then can be pulled through a D and a U
which then can be pulled through a D and a U.
- Polyphase decomposition of H(z) [VAIDYANATHAN 4.3, PROAKIS 10.5.2,
CROCHIERE 3.3]:
@@ -1197,14 +1211,32 @@ c) s-plane and z-plane
down sampling 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.
. Commutator direction from oldest phase (q = Q-1 at end of delay line) to
current phase (q = 0 direct path). 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:
- down sampling input commutator rotates counter clockwise and yields 1
sample every rotation, because the summation stage is combinatorial
- up sampling output commutator rotates clockwise and yields U samples
every rotation
. 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].
. Commutator model [VAIDYANATHAN 4.3.4, CROCHIERE 3.3.2, HARRIS Fig 5.4,
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
oldest so first phase in the FIR sum.
- down sampling input commutator rotates counter clockwise and yields
one output sample every rotation, because the summation stage is
combinatorial. The down sampled output depends on Q input samples from
the past. The output appears whenever the switch is back 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.
. Type I polyphase representation, based on delays z^(-q), yielding counter
clockwise commutator
@@ -1220,7 +1252,7 @@ c) s-plane and z-plane
z^-(Q-1) E_{Q-1}(z^Q)
Q-1
= sum z^(-q) Eq(z^Q) [VAIDYANATHAN Eq 4.3.7]
= sum z^(-q) Eq(z^Q) [VAIDYANATHAN Eq 4.3.7, HARRIS Eq 6.6]
q=0
where Eq(z^Q) is the z-transform of eq[n]:
@@ -1282,6 +1314,9 @@ c) s-plane and z-plane
h[n] --> z^(-q) --> Q:1 --> eq[n]
- Variable bandwidth PFB [HARRIS, FARZAD, CHEN]
. Analysis --> select channels --> synthesis
. Requires Ros > 1 for no in band distortion, typically Ros = 2
- Fractional time delay [CROCHIERE 6.3]
. Up sampling Q --> LPF --> z^(-d) --> down sampling Q yields semi allpass
@@ -1326,8 +1361,76 @@ 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]
. I[n] = cos(-w_k n) x[n]
. 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:
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]
. Down-conversion followed by a LPF is equivalent to a BPF followed by a
down-conversion:
exp(-j w_k n) --> H_LPF(z) = H_BPF(z) --> exp(-j w_k n)
where:
+inf
H_BPF(z) = sum h_lpf[n] exp(j w_k n) z^-n # mix LPF up to BPF
n=-inf
+inf
= sum h_lpf[n] (z exp(-j w_k)^-n
n=-inf
= 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:
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)
13) Quadrature Mirror Filter (QMF) [CROCHIERE 7.7, PROAKIS 10.9.6]
. 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_BPF(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 / D, so that they alias to baseband by the sample
rate change.
- 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)
= H0(z^M) +
z^-1 H1(z^M) +
z^-2 H2(z^M) +
... +
z^-(M-1) H_{M-1}(z^M)
M-1
= sum z^(-q) Hq(z^M) [VAIDYANATHAN Eq 4.3.7, HARRIS Eq 6.6]
q=0
. Apply Noble identity:
14) Polyphase filterbank [HARRIS 6]
15) Quadrature Mirror Filter (QMF) [CROCHIERE 7.7, PROAKIS 10.9.6]
|-- h0[n] --> Down Q --> x0[m] --> Up Q --> f0[n] --|
x[n] --| +--> x^[n]
@@ -1425,6 +1528,54 @@ Appendix A) Signal operators [JOS1 7.2]
Appendix B: Teaser talk: Oversampled polyphase filterbank (OPFB)
1 Introduction
- Repeatedly applying the DFT in time to get timeseries samples of the
subbands
- Size of DFT yields number of subbands
- For critically sampled PFB the sample rate of the subbands just fits the
Nyquist criterium of fs = BWsub
- Oversampled PFB to:
. cover signals near edges of the subbands
. better reconstruct original signal from subbands
2 Data processing approach to understand the OPFB
- Downsampled STFT filterbank [JOS4 9.8]
- Repeatedly window section of the signal and apply DFT to block of the
signal to get timeseries per subband
- Window larger than DFT block size yields the sharper subbands.
- Increment in steps of less then block size for oversampled subband series
3 Use multirate signal processing approach to understand the OPFB
- Downsampling D: LPF and then dicard D-1 samples
- Upsampling U: insert U-1 zeros and LPF
- Goal in implemention is to not calculate samples that will be discarded or
are zero
- 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
- LPF selects the baseband signal, BPF can select any band at multiple of
subband rate fs / Q
- Equivalence criterium mixer + LPF + D = BPF + D
5 Single channel mixer
6 Filterbank for D = U = Q channels --> DFT section
7 Oversampled filterbank
- extra phase rotation after DFT is shift before DFT, because of DFT shift
theorem (for phasors in DFT a shift in time is same as a phase rotation)
- fractional Ros
8 Reconstruction
- Aliasing 0, response 1
9 Applications
- Variable bandpass filter in telecom, resource usage compared to FIR BPF is
much less, so more complexity in algorithm, but much more efficient in
implementation.
- Reconstruction of beamformed data
10 Conclusion
- I still only understand parts (e.g. I do not understand what they mean
with paraunitary, PFB viewed as WOLA, STFT)
11 References: HARRIS, CROCHIERE
https://learning.anaconda.cloud/
https://realpython.com/python-scipy-fft/