Resolve RTSD-271 "C"

Closes RTSD-271

applications/lofar2/model/pfb_os/dsp_study_erko.txt

@@ -160,7 +160,7 @@ a) DTFT [LYONS 3.52, MATLAB]
 
 - Discrete time to continuous frequency domain:
 
-           +inf
+          +inf
     X(w) = sum x[n] exp(-jw n)
           n=-inf
 
@@ -168,7 +168,7 @@ a) DTFT [LYONS 3.52, MATLAB]
 b) z-transform [LYONS 6.3, MATLAB]
 - Decrete time to z-domain:
 
-           +inf
+          +inf
     X(z) = sum x[n] z^-n,  z = r exp(jw) = r (cos(w) + j sin(w)),
           n=-inf
 
@@ -563,7 +563,7 @@ c) s-plane and z-plane
 - For FIR b = h. For IIR it is not possible to directly derive b, a from h
   [LYONS 6.1]. Therefor use z-transform [LYONS 6.3]:
 
-           +inf             +inf
+          +inf             +inf
     H(z) = sum h[n] z^-n  = sum h[n] r^-n exp(-j w n)
           n=-inf           n=-inf
 
@@ -830,30 +830,47 @@ c) s-plane and z-plane
       X(m) = sin(pi * m) / sin(pi * m / K)
           ~= K * sinc(m) for K = N >~ 10
 
-- Fourier transform theorems [JOS4 B]
+- DTFT properties [JOS4 B, PROAKIS 4.3]
+  . Linearity: a1 x1[n] + a2 x2[n] <==> a1 X1(w) + a2 X2(w)
   . Scaling: x(t / a) <==> |a| X(a w)
-  . Shift: x(t - T) <==> exp(-j w T) X(w)
-  . Modulation: x(t) exp(j v t) <==> X(w - v), is dual of shift
+  . Time shift: x(t - T) <==> X(w) exp(-j w T)
+                x[n - k] <==> X(w) exp(-j w k), t = n 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)]
+  . Conjugation: x*[n] <==> X*(-w)
   . Convolution:
       x * y <==> X Y
       x y <==> 1 / (2 pi) X * Y
-  . flip(x) <==> flip(X)
-  . d(t) <==> 1, dirac pulse with area 1 at t = 0
+  . 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).
+  . d[n] <==> 1, dirac pulse with area 1 at n = 0
+    d[n - k] <==> exp(-j w k), dirac pulse with area 1 at n = k
 
                   +inf
-  . d_train_P(t) = sum d(t - m P), period P
-                  m=-inf
+  . d_train_P[n] = sum d(n - k P), period P
+                  k=-inf
     <==>
-                         +inf
-    d_train_P(f) = 1 / P  sum d(f - m / P)
-                         m=-inf
+                            +inf
+    d_train_P(w) = 2 pi / P  sum d(w - 2 pi k / P)
+                            k=-inf
+
+  . Sampling: x_d(t) = x_a(t) d_train_Ts(t)
+              <==>                             +inf
+              X_d(f) = X_a * d_train_fs(f) = fs sum X_a(f - k fs)
+                                               k=-inf
 
-  . sampling: x_d(t) = x(t) d_train_Ts(t)
-              <==>
-              X_d(f) = X * d_train_fs(f)
-                         +inf
-                     = fs sum X(f - k fs)
-                         k=-inf
+    - The sampling theorem [PROAKIS 4.2.9, CROCHIERE 2.1]:
+        The digital spectrum is a periodic repetition of the scaled analogue
+        spectrum with period fs. If spectrum X_a = 0 for |f| >= B, then for
+        fs >= 2 B there is no overlapping aliasing (= spectral folding) and
+        then it is possible to reconstruct x_a from x_d using an LPF.
+    - The sinc() is the ideal interpolation formula:
+
+              +inf        sin(pi (t - nT) / T)  +inf
+      x_a(t) = sum x_d[n] -------------------- = sum x_d[n] sinc((t - nT) / T)
+              n=-inf        (pi (t - nT) / T)   n=-inf
 
 
 10) Short Term Fourier Transform (STFT) [JOS4 7, 8]
@@ -979,9 +996,15 @@ c) s-plane and z-plane
 12) Multirate processing:
 - Linear Time Variant (LTV) process, because it depends on when the
   downsampling and upsampling start.
+- Sampling and sampling rate conversion can be viewed as a modulation process
+  in which the spectrum of the digital signal contains periodic repetitions of
+  the baseband signal (images) spaced at harmonics of the sampling frequency.
+  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].
 - 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 protype filter.
+  the prototype filter.
 - For large D or U use two stage D = D1 * D2 or U = U1 * U2, where D1 > D2 and
   U1 < U2 [LYONS 10.8.2]
 
@@ -994,48 +1017,59 @@ c) s-plane and z-plane
 
 - Noble identities [LYONS Fig 10.20], [VAIDYANATHAN Fig 4.2.3]
 
-  up sampling   : x[n] --> up Q --> H(z^Q) --> y[m], is equivalent to:
-                           H(z) --> up Q
+  . Up sampling   : x[n] --> up Q --> H(z^Q) --> y[m], is equivalent to:
+                             H(z) --> up Q
 
-  down sampling : x[n] --> H(z^Q) --> down --> y[m], is equivalent to:
-                  x[n] --> down Q --> H(z) --> y[m]
+  . Down sampling : x[n] --> H(z^Q) --> down --> y[m], is equivalent to:
+                    x[n] --> down Q --> H(z) --> y[m]
 
-  . Hi(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
+  . 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:
-  . Do not calculate samples that will be thrown away.
-  . Discarding samples folds the spectrum, first the LPF has to remove all
-    folds.
-  . Sequence w(m) is an upsampled-by-D version of sequence x(n), and sequence
-    x(n) is a downsampled-by-D version of sequence w(m) [LYONS 10.9].
-    . Downsampling: W(z) = X(z^D)
+- LPF + downsampling = decimation [LYONS 10.9, PROAKIS 10, CROCHIERE Fig 3.2]:
+  . Downsampling: Y(z) = X(z^(1/Q))
 
-      w(m) = x(m / D), when m is multiple of D else 0
+      y(n) = x(n Q), because y removes Q-1 values from x
 
-              +inf             +inf               +inf
-      W(z) =  sum w(m) z^-m =  sum w(Dk) z^-Dk =  sum x(k) z^-Dk = X(z^D)
-             m=-inf           k=-inf             k=-inf
+             +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
 
-    . Upsampling: W(z^1/D) = X(z)
+  . Spectrum, evaluate Y(z) on unit circle [PROAKIS Eq 10.2.9]:
 
-      w(u) = x(m), when u is Dm else 0
+                 Q-1
+      Y(w) = 1/Q sum H((w - 2 pi k) / Q) X((w - 2 pi k) / Q), w = w_y
+                 k=0
 
-              +inf             +inf               +inf
-      W(z) =  sum w(u) z^-u =  sum w(Dm) z^-Dm =  sum x(m) z^-Dm = X(z^D)
-             u=-inf           m=-inf             m=-inf
+  . Discarding samples folds the spectrum around multiples of pi / Q =
+    fxNyquist / Q = fyNyquist. 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
+
+      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
+
+            +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
+
+  . Spectrum, evaluate Y(z) on unit circle [PROAKIS Eq 10.3.3]:
+
+      Y(w) = Q X(w Q), w = w_y
+
+  . 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.
   . Do not calculate samples that will be inserted as zeros.
-  . Inserting zeros replicates the spectrum, the LPF remove all replicas and by
-    that it interpolates to fill in the 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
     pass band of sin(x)/x [LYONS 10.5.1]
 
 - Fractional time delay [CROCHIERE 6.3]
-  . Up sampling M --> LPF --> z^(-L) --> down sampling M yields semi
-    allpass filter and delay of L / M samples
+  . Up sampling Q --> LPF --> z^(-d) --> down sampling Q yields semi allpass
+    filter and delay of d / Q samples.
 
 - Oversampling ADC and DAC
   . Every oversampling factor of 4 yields 1 extra bit, because then 1 / 4 of the
@@ -1074,6 +1108,54 @@ c) s-plane and z-plane
     output.
 
 
+13) 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]
+          |-- h1[n] --> Down Q --> x1[m] --> Up Q --> f1[n] --|
+
+           Q = 2
+
+          X^(z) = T(z)X(z) + A(z)X(-z), with:
+                . T(z) = H0(z)F0(z) + H1(z)F1(z), transfer part
+                . A(z) = H0(-z)F0(z) + H1(-z)F1(z), aliasing part
+
+- Choose:
+    h0[n] = h[n]           <==>  H0(w) = H(w), prototype LPF
+    h1[n] = (-1)^n h[n]    <==>  H1(w) = H(w - pi), mirror image HPF
+    f0[n] = Q h[n]         <==>  F0(w) = Q H(w)
+    f1[n] = -Q (-1)^n h[n] <==>  F1(w) = -Q H(w - pi), to eliminate aliasing,
+                                 so A(w) = 0
+
+    then: X^(w) = T(w) X(w), with T(w) = H^2(w) - H^2(w - pi)
+
+- Get HPF from LPF using frequency shift (complex modulation) by pi = fNyquist,
+  so: h1[n] = h[n] exp(+j pi n) <==> H(w - pi)
+            = h[n] cos(j pi n)
+            = h[n] (-1)^n
+- For perfect reconstruction T(w) = 1. This can only be achieved for a two tap
+  FIR filter, because Q = 2, so each phase then becomes a delay.
+- Choose linear phase (= symmetric) FIR filter:
+
+   h[n] = h[N - 1 - n), for n = 0,1,...,N-1
+
+   H(w) = Hr(w) exp(-j w (N - 1) / 2), where Hr is real function, so:
+
+      Hr^2(w) = |H(w)|^2
+
+      H^2(w)      =            |H(w)     |^2 exp(-j w (N - 1))
+      H^2(w - pi) = (-1)^(N-1) |H(w - pi)|^2 exp(-j w (N - 1))
+
+      T(w) = |H(w)|^2 - (-1)^(N-1) |H(w - pi)|^2
+
+   For N is odd, T(pi / 2) = 0, therefore choose N is even:
+
+      T(w) = |H(w)|^2 + |H(w - pi)|^2
+
+   For approximate reconstruction optimize for both maximum attenuation in
+   stop band of H(w) and all pass for T(w).
+
+
 Appendix A) Signal operators [JOS1 7.2]
 
 - Operator(x) is element of C^N for all x element of C^N