Updates.

applications/lofar2/model/pfb_os/dsp_study_erko.txt

@@ -91,19 +91,25 @@
     determine the full response of a system, be it stable or unstable,
     including transient parts.
 
+
 2) Windows [JOS4 3]
-- Tabel [PROAKIS 8.2.2]
+- Window types table [HARRIS 3.2, PROAKIS 8.2.2]
+- Kaiser beta and side lobe levels figure [HARRIS 3.2]
 - Rectangular window with length M [LYONS 3.13]
   . Dirichlet kernel or aliased sinc:
       HF(m) = c sin(pi * m * M / Ndtft) / sin(pi * m / Ndtft)
   . Ndtft = M yields all-ones form that defines the DFT frequency response
     to an input sinusoidal and it is also the of a single DFT bin:
       HF(m) = sin(pi * m) / sin(pi * m / M)
-           ~= Ndtft * sinc(pi * m) for Ndtft = M >~ 10
+           ~= Ndtft * sinc(m) for Ndtft = M >~ 10
 - 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.
     3.45]
-  . Main lobe width is 4 pi / M
+  . Main lobe width is 4 pi / M between zero crossings.
+  . Main lobe cutoff frequency is about fpass ~= 0.6 / M, where fpass is
+    positive frequencies -6 dB power BW [EK]
   . As M increases, the main lobe narrows (better frequency resolution).
   . M has no effect on the height of the side lobes (same as the Gibbs
     phenomenon for truncated Fourier series expansions.
@@ -113,6 +119,7 @@
     the window transform is real in the zero-phase case (i.e., centered
     about time 0).
 
+
 3) Low pass filter (LPF)
 - Design parameters [JOS4 4.2]
   . Pass band edge frequency: w_pass
@@ -123,28 +130,51 @@
       r_stop = 10**(r_stop_dB / 20)
 - Ideal LPF [JOS4 4.1]
   . w_cutoff = w_pass = w_stop, r_pass = r_stop = 0
-  . sinc(t) = sin(pi t) / (pi t) [JOS4 3.1, numpy)
-  . h_ideal(n) = 2 f_c sinc(2 f_c n), n in Z
-    - f_c is normalized cutoff frequency with fs = 1
+  . sinc(t) = sin(pi t) / (pi t) [JOS4 3.1, numpy, MATLAB]
+  . 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
+  . 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.
 - 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)
+      with default Hamming).
     - Least squares designs (scipy.signal.firls(), MATLAB firls, fircls1)
       . firls = least squares
       . fircls, fircls1 = constrained ls with pass, stop ripple
     - The Remez or Parks-McClellan algorithm (scipy.signal.remez(), MATLAB
       firpm)
-  . MATLAB filters yield n + 1 coefs
+  . MATLAB filters specify order n, so yield n + 1 coefs
   . LS and Remez can do bandpass (= flat), differentiator, hilbert
-  . Linear phase filter types (filter order is Ntaps - 1, fNyquist = f2/2):
+  . LS and Remez for large Ntaps (> 1000) can fail, but can be achieved
+    using Fourier interpolation a filter with less coefs.
+    by stuffing zeros in the frequency domain [LYONS 13.27, 13.28]
+  . Linear phase filter types (filter order is Ntaps - 1, fNyquist = fs/2)
+    [LYONS 9.4.1]:
       Type  Ntaps  Symmetry  H(0) H(fs/2)
       I     Odd    Even      any  any  --> LPF, HPF
       II    Even   Even      any  0    --> LPF
       III   Odd    Odd       0    0    --> differentiator, hilbert
       IV    Even   Odd       0    any  --> differentiator, hilbert
+  . 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(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(k pi), for fc = fs / 2,
+              series [1, -1]
+
 
 2) Finite Impulse Response (FIR) filters
 - FIR filters perform time domain Convolution by summing products of shifted
@@ -163,20 +193,30 @@
         \----------\-- ... ------------\--> + --> y(n)
 
 - Convolution in time domain is equivalent to multiplication in frequency
-  domain
-  y(n) = h(k) * x(n) ==> DFT ==> Y(m) = H(m) X(m)
+  domain:
+    y(n) = h(k) * x(n) ==> DFT ==> Y(m) = H(m) X(m)
 
 - Number of FIR coefficients (Ntaps)
   . Trade window main-lobe width for window side-lobe levels and in turn filter
     transition bandwidth and side-lobe levels
   . Transition bandwidth: df = fstop - fpass
-  . Window based design [HARRIS 3.2, LYONS 5.10.5]:
+  . Ntaps must be odd if a passband includes fNyquist = fs / 2 [scipy firwin].
+  . Kaiser window based design [HARRIS 3.2, Fig. 3.8 for beta]:
     - Ntaps ~= fs / df * (Atten(dB) - 8) / 14
   . Remez = Parks-McClellan [HARRIS 3.3, LYONS 5.6]:
     - yield a Chebychev-type filter
     - Steeper transition than window based, but constant stopband peak levels
-    - Ntaps = f(fs, fpass, fstop, passband ripple +-d1, stopband ripple +-d2)
-           ~= fs / df * Atten(dB) / 22
+    - Ntaps = func(fs, fpass, fstop, pb ripple, sb ripple) [HARRIS 3.3,
+      LYONS 5.10.5]:
+           ~= fs / df * Atten(dB) / 22, df = abs(fstop - fpass)
+    - Choose transition region specification in order 4 pi / Ntaps, and not
+      too wide, because then the transition band is not smooth [JOS4 4.5.2].
+  . Equiripple vs 1/f ripple [HARRIS 3.3.1], rate of decay depends on order of
+    the discontinuity in time domain:
+    . 0-order 1/f decay implies discontinuous amplitude like with rect
+      window, this is advantagous in case of multi rate where multiple
+      regions fold (alias) into the pass band.
+    . -1 order no decay implies impulses at end-points of impulse response.
 
 - Linear phase FIR filter
   . Even or odd symmetrical h(n) = +-h(M - 1 - n), n = 0,1,...,N-1
@@ -186,12 +226,74 @@
   . Design using windowed sinc, because windows are symmetrical [PROAKIS 8.2.2,
     LYONS 5.3.2, DSPGUIDE 16]
 
-- Half band FIR filter [LYONS 5.7]
-  . Symmetrical frequency response about fs / 2, so fpass + fstop = fs / 2
+- Half band FIR filter [LYONS 5.7, HARRIS 8]
   . When Ntaps is odd, then half of the coefs are 0.
-
-
-3) Discrete Fourier Transform (DFT)
+    - Ntaps + 1 is multiple of 4
+    - Gain at fs / 4 is 0.5
+  . Symmetrical frequency response about fc = fs / 4, so fpass + fstop = fs / 2
+  . Design:
+    - window: h_lp[n] * w[n] [HARRIS 8.3]
+    - remez:
+      . same weigths for stop and pass band [HARRIS 8.4]
+      . construct 2N + 1 half band from N (non zero) half band design, trick
+        to use different weights [HARRIS 8.5]
+  . Low pass, fc = fs / 4, sinc(t) = sin(pi t) / (pi t) [numpy]
+    . h_lp[n] = h_ideal(n) = 2 f_cutoff sinc(2 f_cutoff n) = 1/2 sinc(n / 2), n in Z
+    . h_hp[n] = h_lp[n] exp(j pi n) = h_lp[n] * cos(n pi)
+  . Low pass and high pass half band filters are complementary [HARRIS 8.2],
+    so h_lp[n] + h_hp[n] = d(n - N), 0 <= n <= 2N (causal)
+
+
+4) Hilbert transform (HT) and analytic signal [LYONS 9]
+- The Hilbert transform (HT) creates a 90 degrees phase-shifted version of a
+  real signal. The HT can be defined as a filter with impulse response ht,
+  . x_ht[t] = xi[t] = ht[t] * xr[t].
+    - Analytical signal: xa = xr + j xi
+    - E.g.: xr = cos(wt) --> xi = sin(wt) --> xa = exp(jwt)
+  . HT(w) = j, 0, -j for f < 0, 0, > 0
+    - Analytic signal is zero for negative frequencies, and has only the
+      positive frequencies of xr, but with 2x magnitude. Therefore Xa(w) is
+      also called a one sided spectrum (of Xr(w)).
+  . ht[t] = IFT(HT(w)), yields:
+    ht[n] = 0, when n is even
+          = 2 / (pi n Ts), when n is odd
+- Purpose analytic signal:
+  . measuring instantaneous characteristics (magnitude sqrt(xr^2 + xi^2),
+    phase atan(xi / xr), frequency d/dt(phase)) of time signals, e.g. for AGC,
+    envelop, demodulation
+  . frequency translation.
+- Anti-symmetric response:
+  . Linear phase, with discontinuity of pi at f = 0, as is the purpose of HT
+  . Type IV, Ntaps even --> HT(0) = 0, HT(fs/2) = any
+  . Type III, Ntaps odd --> HT(0) = 0, HT(fs/2) = 0 --> so BPF, advantages:
+    - odd coefs are zero
+    - group delay is (Ntaps - 1) / 2 is integer
+- Design methods:
+  . HT window design: ht[n] * w[n], best use Ntaps is odd
+  . Analytic signal generation using complex BPF [LYONS 9.5, HARRIS 8.5]:
+      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_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
+        feasible
+      . h_cos is symmetrical, h_sin is anti-symmetrical
+      . if BPF fcenter = fs / 4 then mixer = exp(pi / 2 (k - D)), so then half
+        of h_cos and h_sin coefs are 0, alternately
+      . if LPF BW is half band, so f_cutoff = fcenter, then half of h_sin coefs
+        are 0 and all but one of h_cos coefs are 0.
+  . Analytic signal generation using complex down conversion [LYONS 8]
+  . Using the DFT [LYONS 9.4.2], like with signal.hilbert. The signal.hilbert
+    generate the analytic signal of a real signal. It uses the DFT of the real
+    signal to determine the analytic signal xa = IDFT(DFT(xr) * 2U) =
+    xr + j ht. With xr = d(n - Ntaps // 2) this yields ht = imag(xa).
+  . Half band [13.1, 13.37]
+
+- IIR Hilbert transform [HARRIS 10.7]
+
+5) Discrete Fourier Transform (DFT)
 - The N roots of unity [JOS1 3.12, 5.1, PROAKIS 5.1.3, LYONS 4.3]. Note JOS
   uses +j in W_N because inproduct is with conj(W_N), others use -j because
   then W_N can be used directly in equation and matrix:
@@ -293,9 +395,9 @@
   . K = N yields all-ones form that defines the DFT frequency response
     to an input sinusoidal and it is also the of a single DFT bin:
       X(m) = sin(pi * m) / sin(pi * m / K)
-          ~= K * sinc(pi * m) for K = N >~ 10
+          ~= K * sinc(m) for K = N >~ 10
 
-4) Multirate processing:
+6) Multirate processing:
 - Linear Time Variant (LTV) process, because it depends on when the
   downsampling and upsampling start.
 - Polyphase filtering ensures that only the values that remain are calculated,
@@ -318,7 +420,7 @@ Upsampling + LPF = interpolation:
   pass band of sin(x)/x [LYONS 10.5.1]
 
 
-5) Signal operators [JOS1 7.2]
+Appendix 7) Signal operators [JOS1 7.2]
 
 - Operator(x) is element of C^N for all x element of C^N
   . assume modulo N indexing for n in x(n), so x(n) = x(n + mN) or periodic

applications/lofar2/model/pfb_os/filter_design_firls.ipynb

@@ -144,6 +144,7 @@
     "#   and the weight. Choose 0 ~< w_tb ~< 1.0 fpass, to ensure the FIR filter design converges\n",
     "#   and improve the passband ripple and stopband attenuation. A to large transition band \n",
     "#   also gives the design too much freedom and causes artifacts in the transition.\n",
+    "# . scipy firls() defines fpass relative to fs, so can use fpass as cutoff frequency.\n",
     "Q = Ntaps\n",
     "N = Ncoefs // Q + 1  # + 1, because firls only supports odd number of FIR coefficients\n",
     "f_pb = fpass * Q # pass band cut off frequency\n",