Improve multirate section.

applications/lofar2/model/pfb_os/dsp_study_erko.txt

@@ -31,6 +31,7 @@
 # * [JOS2] Introduction to Digital Filters, 2007
 # * [JOS3] Physical Audio Signal Processing, 2010
 # * [JOS4] Spectral Audio Signal Processing, 2011
+# * [SP4COMM] Signal Processing for Communications, 2008, Paolo Prandoni and Martin Vetterli
 #
 # * [WIKI] https://en.wikipedia.org/wiki/Bilinear_transform
 # * [WIKI] https://en.wikipedia.org/wiki/Discrete_cosine_transform
@@ -61,6 +62,8 @@
 #     . "Realisering van Digitale Signaalbewerkende Systemen, Toepassingen",
 #        5N290, TUE, P.C.M. Sommen, --> DAC slide 19, 20,
 # * [PM-REMEZ] https://pm-remez.readthedocs.io/en/latest/
+# * [SELESNICK] Ivan Selesnick
+#   . https://eeweb.engineering.nyu.edu/iselesni/EL713/zoom/mrate.pdf
 #
 # https://ocw.mit.edu/courses/6-341-discrete-time-signal-processing-fall-2005/
 # Youtube: Guitars 4RL
@@ -1032,42 +1035,92 @@ c) s-plane and z-plane
   . H_i(z^Q) is upsampled-by-Q version of H(z), so with Q-1 zero coefficients
     in the H(z) power series, starting at phase i.
 
-- LPF + downsampling = decimation [LYONS 10.9, PROAKIS 10, CROCHIERE Fig 3.2]:
-  . Downsampling: Y(z) = X(z^(1/Q))
+- LPF + downsampling = decimation [LYONS 10.9, PROAKIS 10, CROCHIERE Fig 3.2,
+  VAIDYANATHAN Fig 4.1.4, SP4COMM 11.1.2]:
+  . Downsampling:
 
-      y(n) = x(n Q), because y removes Q-1 values from x
+      x_D[n] = x[nD], because x_D removes D-1 values from x
 
-             +inf             +inf              +inf
-      Y(z) =  sum y(n) z^-n =  sum x(Qn) z^-n =  sum x(k) z^-k/Q = X(z^(1/Q))
-             n=-inf           n=-inf            k=-inf
+    Define x' at x time grid, but with x' = 0 for samples in x that will be
+    discarded:
 
-  . Spectrum, evaluate Y(z) on unit circle [PROAKIS Eq 10.2.9]:
+      x'[n] = d[n] x[n],   with d[n] = 1, for n % D = 0
+                                     = 0, otherwise
 
-                 Q-1
-      Y(w) = 1/Q sum H((w - 2 pi k) / Q) X((w - 2 pi k) / Q), w = w_y
+    Use x' to first express z-transform X_D(z) in X'(z) and then z-transform
+    of X'(z) in X(z). The z-transform X_D(z) = X'(z^(1/D)), because x' is 0
+    when k is not a multiple of D:
+
+               +inf
+      X_D(z) = sum x_D[n] z^(-n)
+              n=-inf
+
+               +inf                +inf
+             = sum x'[nD] z^(-n) = sum x'[k] z^(-k/D) = X'(z^(1/D))
+              n=-inf              k=-inf
+
+    Make use of the identity that the selector d[n] is equal to the normalized
+    sum of the D roots of unity:
+
+                 D-1
+      d[n] = 1/D sum W_D^(-kn),  with W_D = exp(-j 2pi/D)
+                 k=0
+
+    so:
+
+              +inf                       D-1  +inf
+      X'(z) = sum d[n] x[n] z^(-n) = 1/D sum  sum x[n] (z W_D^k)^(-n)
+             n=-inf                      k=0 n=-inf
+
+                  D-1
+            = 1/D sum X(z W_D^k)
+                  k=0
+
+    Combining the results yields the z-transform of downsampled signal in
+    terms of z-transform of the original signal:
+
+                   D-1
+      X_D(z) = 1/D sum X(z^(1/D) W_D^k)
                    k=0
 
-  . Discarding samples folds the spectrum around multiples of pi / Q =
-    fxNyquist / Q = fyNyquist. First the LPF has to remove all folds.
+    The Fourier transform of the downsampled signal is obtained by evaluating
+    X_D(z) on the unit circle with z = exp(jw) [PROAKIS Eq 10.2.9]:
+
+                      D-1
+      X_D(e^jw) = 1/D sum X(exp(j (w - k 2pi) / D))
+                      k=0
+
+    The resulting spectrum is the scaled sum of D superimposed copies of the
+    original spectrum X(e^jω), and each copy is shifted in frequency by a
+    multiple of 2pi/D and the result is stretched by a factor of D. This is
+    similar to sampling of an analogue signal that creates a periodization of
+    the analogue spectrum, but now the spectra are already inherently
+    2pi periodic, and downsampling creates D − 1 additional interleaved copies.
+    For the spectral copies not to overlap, the maximum (positive) frequency
+    the original spectrum must be less than pi / D. This is the non-aliasing
+    condition for the downsampling operator.
+
+  . Discarding samples folds the spectrum around multiples of pi / D =
+    fxNyquist / D. First the LPF has to remove all folds.
   . Do not calculate samples that will be thrown away.
 
 - Upsampling + LPF = interpolation:
-  . Upsampling: Y(z) = X(z^Q), because y inserts Q-1 zeros in x
+  . Upsampling:
 
-      y[m] = x[m / Q], when m is multiple of Q, else 0, so equivalently
-      y[m] = x[n], when m = Q n else 0
+      x_U[n] = x[n / U], for n % U = 0
+             = 0, otherwise, because x_U inserts U-1 zeros in x
 
-            +inf            +inf              +inf
-      Y(z) = sum y(m) z^-m = sum y(Qn) z^-Qn = sum x(n) z^-Qn = X(z^Q)
-            m=-inf          n=-inf            n=-inf
+              +inf              +inf
+      X_U(z) = sum x_U[n] z^-n = sum x[k] z^-kU = X(z^U), with n = kU
+              n=-inf            k=-inf
 
-  . Spectrum, evaluate Y(z) on unit circle [PROAKIS Eq 10.3.3]:
+  . Spectrum, evaluate X_U(z) on unit circle [PROAKIS Eq 10.3.3]:
 
-      Y(w) = Q X(w Q), w = w_y
+      X_U(e^jw) = X(exp(jw U))
 
-  . Inserting zeros replicates the spectrum around multiples of pi / Q =
-    fyNyquist / Q = fxNyquist. Then the LPF has to remove all replicas and
-    by that it interpolates to fill in the zeros.
+  . Inserting zeros replicates the spectrum around multiples of pi / U. Then
+    the LPF has to remove all replicas and by that it interpolates to fill in
+    the zeros.
   . Do not calculate samples that will be inserted as zeros.
   . Using zero order hold would be a naive approach, because then all samples
     need to be calculated and the LPF then needs to compensate for the non-flat