From 6d9fdfee2efdab8a67831ba354b1cd61cbb83faa Mon Sep 17 00:00:00 2001
From: Eric Kooistra <kooistra@astron.nl>
Date: Wed, 21 Aug 2024 16:27:24 +0200
Subject: [PATCH] Add NYQUIST(L).
---
.../lofar2/model/pfb_os/dsp_study_erko.txt | 217 +++++++++++++++---
1 file changed, 184 insertions(+), 33 deletions(-)
diff --git a/applications/lofar2/model/pfb_os/dsp_study_erko.txt b/applications/lofar2/model/pfb_os/dsp_study_erko.txt
index 4aa7d564a5..4cff475aef 100644
--- a/applications/lofar2/model/pfb_os/dsp_study_erko.txt
+++ b/applications/lofar2/model/pfb_os/dsp_study_erko.txt
@@ -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/
--
GitLab