Merge branch 'RTSD-162' into 'master'

Resolve RTSD-162

Closes RTSD-162

See merge request !368

applications/lofar2/model/pfb_os/dsp.py

@@ -58,10 +58,21 @@ def pow_db(volts):
 
 
 def is_even(n):
-    """Return True if n is even, else False when odd."""
+    """Return True if n is even, else False when odd.
+
+    For all n in Z.
+    """
     return n % 2 == 0
 
 
+def is_odd(n):
+    """Return True if n is odd, else False when even.
+
+    For all n in Z, because result of n % +2 is 0 or +1.
+    """
+    return n % 2 == 1
+
+
 def is_symmetrical(x, anti=False):
     """Return True when x[n] = +-x[N-1 - n], within tolerances, else False."""
     rtol = c_rtol
@@ -90,21 +101,53 @@ def read_coefficients_file(filepathname):
     return coefs
 
 
+def one_bit_quantizer(x):
+    """Quantize 0 and positive x to +1 and negative x to -1."""
+    return np.signbit(x) * -1 + 2
+
+
+def impulse_at_zero_crossing(x):
+    """Create signed impulse at zero crossings of x."""
+    diff = np.diff(one_bit_quantizer(x))
+    return np.concatenate((np.array([0]), diff))
+
+
 ###############################################################################
 # Filter design
 ###############################################################################
 
+def nof_taps_kaiser_window(fs, fpass, fstop, atten_db):
+    """Number of FIR LPF taps using Kaiser window based design
+
+    Reference: [HARRIS 3.2, Fig. 3.8 for beta]
+    """
+    df = fstop - fpass
+    return int((fs / df) * (atten_db - 8) / 14)
+
+
+def nof_taps_remez(fs, fpass, fstop, atten_db):
+    """Number of FIR LPF taps using remez = Parks-McClellan based design.
+
+    Reference: [HARRIS 3.3, LYONS 5.6]
+    """
+    df = fstop - fpass
+    return int((fs / df) * (atten_db / 22))
+
+
 def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0):
     """Derive FIR coefficients for prototype low pass filter using ifft of
     magnitude frequency response.
 
+    The Npoints defines the double sided spectrum, so Npass = 2 * fpass / fs
+    * Npoints, where fpass is the positive cutoff frequency of the LPF.
+
     Input:
     - Npoints: Number of points of the DFT in the filterbank
-    - Npass: Number of points with gain > 0 in pass band
+    - Npass: Number of points with gain > 0 in pass band.
     - bandEdgeGain : Gain at band edge
     Return:
     - h: FIR coefficients from impulse response
-    . f: normalized frequency axis for HF, fs = 1
+    - f: normalized frequency axis for HF, fs = 1
     - HF: frequency transfer function of h
     """
     # Magnitude frequency reponse
@@ -132,6 +175,8 @@ def fourier_interpolate(HFfilter, Ncoefs):
     time shift of hInterpolated, to make it symmetrical. Similar as done in
     pfs_coeff_final.m and pfir_coeff.m. Use upper = conj(lower), because that
     is easier than using upper from HFfilter.
+
+    Reference: LYONS 13.27, 13.28
     """
     N = len(HFfilter)
     K = N // 2
@@ -173,6 +218,58 @@ def fourier_interpolate(HFfilter, Ncoefs):
     return hInterpolated.real
 
 
+###############################################################################
+# Hilbert transform filter
+###############################################################################
+
+def hilbert_response(Ntaps):
+    """Calculate impulse response of Hilbert filter truncated to Ntaps
+    coefficients.
+
+    h(t) = 1 / (pi t) ( 1 - cos(ws t / 2), l'Hoptial's rule: h(0) = 0
+
+    For n = ...-2, -1, 0, 1, 2, ..., and Ts = 1, and Ntaps is
+    . odd: t = n Ts
+        ht[n] = 1 / (pi n)) ( 1 - cos(pi n))
+              = 0, when n is even
+              = 2 / (pi n), when n is odd
+    . even: t = (n + 0.5) Ts
+        ht(n) = 1 / (pi m)) ( 1 - cos(pi m)), with m = n + 0.5
+    """
+    Npos = Ntaps // 2
+    if is_even(Ntaps):
+        ht = np.zeros(Npos)
+        for n in range(0, Npos):
+            m = n + 0.5
+            ht[n] = 1 / (np.pi * m) * (1 - np.cos(np.pi * m))
+        ht = np.concatenate((np.flip(-ht), ht))
+    else:
+        Npos += 1
+        ht = np.zeros(Npos)
+        for n in range(1, Npos, 2):
+            ht[n] = 2 / (np.pi * n)
+        ht = np.concatenate((np.flip(-ht[1:]), ht))
+    return ht
+
+
+def hilbert_delay(Ntaps):
+    """Delay impulse by (Ntaps - 1) / 2 to align with hilbert_response(Ntaps).
+
+    Analytic signal htComplex = htReal + 1j * htImag where:
+    . htReal = hilbert_delay(Ntaps)
+    . htImag = hilbert_response(Ntaps)
+
+    Only support integer delay D = (Ntaps - 1) / 2, so Ntaps is odd, then return
+    htReal, else return None.
+    """
+    if is_even(Ntaps):
+        return None
+    D = (Ntaps - 1) // 2
+    htReal = np.zeros(Ntaps)
+    htReal[D] = 1.0
+    return htReal
+
+
 ###############################################################################
 # DFT
 ###############################################################################
@@ -221,7 +318,7 @@ def dtft(coefs, Ndtft=None, zeroCenter=True, fftShift=True):
 # Plotting
 ###############################################################################
 
-def plot_time_response(h, markers=False):
+def plot_time_response(h, name='', markers=False):
     """Plot time response (= impulse response, window, FIR filter coefficients).
 
     Input:
@@ -232,7 +329,7 @@ def plot_time_response(h, markers=False):
         plt.plot(h, '-', h, 'o')
     else:
         plt.plot(h, '-')
-    plt.title('Time response')
+    plt.title('Time response %s' % name)
     plt.ylabel('Voltage')
     plt.xlabel('Sample')
     plt.grid(True)
@@ -247,7 +344,7 @@ def plot_spectra(f, HF, fs=1.0, fLim=None, dbLim=None):
     . fs: sample frequency in Hz, scale f by fs, fs >= 1
     """
     Hmag = np.abs(HF)
-    Hphs = np.angle(HF)
+    Hphs = np.unwrap(np.angle(HF))
     Hpow_dB = pow_db(HF)  # power response
     fn = f * fs
     if fs > 1:

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",