Add NYQUIST(L).

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
@@ -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
@@ -307,11 +317,6 @@ c) s-plane and z-plane
             = 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
@@ -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).
@@ -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))
 
-13) Quadrature Mirror Filter (QMF) [CROCHIERE 7.7, PROAKIS 10.9.6]
+  . 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)
+
+  . 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/