Merge branch 'RTSD-118' into 'master'

Resolve RTSD-118

Closes RTSD-118

See merge request !366

applications/lofar2/model/pfb_os/README.txt

+Author: Eric Kooistra, nov 2023
+
+* Practise DSP [1].
+
+* Try to reproduce LOFAR subband filter FIR coefficients using scipy instead
+  of MATLAB.
+  The pfs_coeff_final.m from the Filter Task Force (FTF) in 2005 use fircls1
+  with r_pass and r_stop to define the ripple. In addition it post applies a
+  Kaiser window with beta = 1 to make the filter attenuation a bit more deep
+  near the transition. The pfir_coeff.m from Apertif also uses fircls1.
+  Both use fircls1 with N = 1024 FIR coefficients and then Fourier
+  interpolation to achieve Ncoefs = 1024 * 16 FIR coefficients. Both scripts
+  can not exactly reproduce the actual LOFAR1 coefficients, therefore these
+  are loaded from a file Coeffs16384Kaiser-quant.dat
+
+* Try low pass filter design methods using windowed sync, firls, remez [3]
+  The windowed sync method, firls leased squares method and remez method all
+  yield comparable results, but firls and remez perform slightly better near
+  the transition band. The firls and remez functions from scipy.signal use
+  transition bandwidth and weights between pass and stop band to influence
+  the transition region and ripple. For remez the ripple is constant in the
+  pass band and stop band, for firls the ripple is largest near the band
+  transition.
+
+* It is possible to design a good FIR filter using Python scipy. Possibly with
+  some extra help of a filter design and analysis (FDA) tool like pyfda [2].
+
+[1] dsp_study_erko.txt, summary of DSP books
+[2] pyfda, dsp, at https://github.com/chipmuenk
+[3] Try FIR filter design methods
+    * dsp.py import for Python jupyter notebooks
+    * filter_design_firls.ipynb
+    * filter_design_remez.ipynb
+    * filter_design_windowed_sync.ipynb

applications/lofar2/model/pfb_os/dsp.py

+#! /usr/bin/env python3
+###############################################################################
+#
+# Copyright 2022
+# ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+###############################################################################
+
+# Author: Eric Kooistra
+# Purpose: Utilities and functions for DSP
+# Description:
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+c_interpolate = 10
+c_atol = 1e-15
+c_rtol = 1e-8   # 1/2**32 = 2.3e-10
+
+
+###############################################################################
+# Utilities
+###############################################################################
+
+def ceil_div(num, den):
+    """ Return integer ceil value of num / den """
+    return int(np.ceil(num / den))
+
+
+def ceil_log2(num):
+    """ Return integer ceil value of log2(num) """
+    return int(np.ceil(np.log2(num)))
+
+
+def ceil_pow2(num):
+    """ Return power of 2 value that is equal or greater than num """
+    return 2**ceil_log2(num)
+
+
+def pow_db(volts):
+    """Voltage to power in dB"""
+    return 20 * np.log10(np.abs(volts) + c_atol)
+
+
+def is_even(n):
+    """Return True if n is even, else False when odd."""
+    return n % 2 == 0
+
+
+def is_symmetrical(x, anti=False):
+    """Return True when x[n] = +-x[N-1 - n], within tolerances, else False."""
+    rtol = c_rtol
+    atol = np.min(np.abs(x[np.nonzero(x)])) * rtol
+    n = len(x)
+    h = n // 2
+    if is_even(n):
+        if anti:
+            return np.allclose(x[0:h], np.flip(-x[h:]), rtol=rtol, atol=atol)
+        else:
+            return np.allclose(x[0:h], np.flip(x[h:]), rtol=rtol, atol=atol)
+    else:
+        if anti:
+            return np.allclose(x[0:h], np.flip(-x[h + 1:]), rtol=rtol, atol=atol) and np.abs(x[h]) < atol
+        else:
+            return np.allclose(x[0:h], np.flip(x[h + 1:]), rtol=rtol, atol=atol)
+
+
+def read_coefficients_file(filepathname):
+    coefs = []
+    with open(filepathname, 'r') as fp:
+        for line in fp:
+            if line.strip():  # skip empty line
+                s = int(line)   # one coef per line
+                coefs.append(s)
+    return coefs
+
+
+###############################################################################
+# Filter design
+###############################################################################
+
+def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0):
+    """Derive FIR coefficients for prototype low pass filter using ifft of
+    magnitude frequency response.
+
+    Input:
+    - Npoints: Number of points of the DFT in the filterbank
+    - 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
+    - HF: frequency transfer function of h
+    """
+    # Magnitude frequency reponse
+    HF = np.zeros([Npoints])
+    HF[0] = bandEdgeGain
+    HF[1 : Npass - 1] = 1
+    HF[Npass - 1] = bandEdgeGain
+    # Zero center HF to make it even
+    HF = np.roll(HF, -(Npass // 2))
+    f = np.arange(0, 1, 1 / Npoints)
+    # Filter impulse response
+    h = np.fft.ifft(HF).real  # imag is 0 for even HF
+    h = np.roll(h, Npoints // 2)
+    return h, f, HF
+
+
+def fourier_interpolate(HFfilter, Ncoefs):
+    """Use Fourier interpolation to create final FIR filter coefs.
+
+    HF contains filter transfer function for N points, in order 0 to fs. The
+    interpolation inserts Ncoefs - N zeros and then performs IFFT to get the
+    interpolated impulse response.
+
+    Use phase correction depenent on interpolation factor Q to have fractional
+    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.
+    """
+    N = len(HFfilter)
+    K = N // 2
+
+    # Interpolation factor Q can be fractional, for example:
+    # - Ncoefs = 1024 * 16
+    #   . firls: N = 1024 + 1  --> Q = 15.98439
+    #   . remez: N = 1024 --> Q = 16
+    Q = Ncoefs / N
+
+    # Apply phase correction (= time shift) to make the impulse response
+    # exactly symmetric
+    f = np.arange(N) / Ncoefs
+    p = np.exp(-1j * np.pi * (Q - 1) * f)
+    HF = HFfilter * p
+
+    HFextended = np.zeros(Ncoefs, dtype=np.complex_)
+    if is_even(N):
+        # Copy DC and K - 1 positive frequencies in lower part (K values)
+        HFextended[0:K] = HF[0:K]  # starts with DC at [0]
+        # Create the upper part of the spectrum from the K - 1 positive
+        # frequencies and fs/2 in lower part (K values)
+        HFflip = np.conjugate(np.flip(HF[0:K + 1]))  # contains DC at [0]
+        HFextended[-K:] = HFflip[0:K]  # omit DC at [K]
+        # K + K = N values, because N is even and K = N // 2
+    else:
+        # Copy DC and K positive frequencies in lower part (K + 1 values)
+        HFextended[0:K + 1] = HF[0:K + 1]  # starts with DC at [0]
+        # Create the upper part of the spectrum from the K positive
+        # frequencies in the lower part (K values)
+        HFflip = np.conjugate(np.flip(HF[0:K + 1]))  # contains DC at [0]
+        HFextended[-K:] = HFflip[0:K]  # omit DC at [K]
+        # K + 1 + K = N values, because N is odd and K = N // 2
+    hInterpolated = np.fft.ifft(HFextended)
+    if np.allclose(hInterpolated.imag, np.zeros(Ncoefs), rtol=c_rtol, atol=c_atol):
+        print('hInterpolated.imag ~= 0')
+    else:
+        print('WARNING: hInterpolated.imag != 0 (max(abs) = %f)' % np.max(np.abs(hInterpolated.imag)))
+    return hInterpolated.real
+
+
+###############################################################################
+# DFT
+###############################################################################
+
+def dtft(coefs, Ndtft=None, zeroCenter=True, fftShift=True):
+    """Calculate DTFT of filter impulse response or window.
+
+    Use DFT with Ndtft points, to have frequency resolution of 2 pi / Ndtft.
+    Define h by extending coefs with Ndtft - M zeros, where M = len(coefs).
+    This DFT approaches the DTFT using bandlimited interpolation, similar to
+    INTERP() which interpolates by factor L = Ndtft / M [JOS1].
+
+    Input:
+    . coefs: filter impulse response or window coefficients
+    . Ndtft: number of points in DFT to calculate DTFT
+    . zeroCenter: when True zero center h to have even function that aligns
+      with cos() in DTFT, for zero-phase argument (+-1 --> 0, +-pi). Else
+      apply h as causal function.
+    . fftShift: when True fft shift to have -0.5 to +0.5 frequency axis,
+      else use 0 to 1.0.
+    Return:
+    . h: zero padded coefs
+    . f: normalized frequency axis for HF, fs = 1
+    . HF: dtft(h), the frequency transfer function of h
+    """
+    M = len(coefs)
+    if Ndtft is None:
+        Ndtft = ceil_pow2(M * c_interpolate)
+    # Time series, causal with coefs at left in h
+    h = np.concatenate((coefs, np.zeros([Ndtft - M])))
+    if zeroCenter:
+        # Zero center h to try to make it even
+        h = np.roll(h, -(M // 2))
+    # DFT
+    HF = np.fft.fft(h)
+    # Normalized frequency axis, fs = 1, ws = 2 pi
+    f = np.arange(0, 1, 1 / Ndtft)  # f = 0,pos, neg
+    if fftShift:
+        # FFT shift to center HF, f = neg, 0,pos
+        f = f - 0.5
+        HF = np.roll(HF, Ndtft // 2)
+    return h, f, HF
+
+
+###############################################################################
+# Plotting
+###############################################################################
+
+def plot_time_response(h, markers=False):
+    """Plot time response (= impulse response, window, FIR filter coefficients).
+
+    Input:
+    . h: time response
+    . markers: when True plot time sample markers in curve
+    """
+    if markers:
+        plt.plot(h, '-', h, 'o')
+    else:
+        plt.plot(h, '-')
+    plt.title('Time response')
+    plt.ylabel('Voltage')
+    plt.xlabel('Sample')
+    plt.grid(True)
+
+
+def plot_spectra(f, HF, fs=1.0, fLim=None, dbLim=None):
+    """Plot spectra for power, magnitude, phase, real, imag
+
+    Input:
+    . f: normalized frequency axis for HF (fs = 1)
+    . HF: spectrum, e.g. frequency transfer function HF = DTFT(h)
+    . fs: sample frequency in Hz, scale f by fs, fs >= 1
+    """
+    Hmag = np.abs(HF)
+    Hphs = np.angle(HF)
+    Hpow_dB = pow_db(HF)  # power response
+    fn = f * fs
+    if fs > 1:
+        flabel = 'Frequency [fs / %d]' % fs
+    else:
+        flabel = 'Frequency [fs]'
+
+    plt.figure(1)
+    plt.plot(fn, Hpow_dB)
+    plt.title('Power spectrum')
+    plt.ylabel('Power [dB]')
+    plt.xlabel(flabel)
+    if fLim:
+        plt.xlim(fLim)
+    if dbLim:
+        plt.ylim(dbLim)
+    plt.grid(True)
+
+    plt.figure(2)
+    plt.plot(fn, HF.real, 'r')
+    plt.plot(fn, HF.imag, 'g')
+    plt.title('Complex voltage spectrum')
+    plt.ylabel('Voltage')
+    plt.xlabel(flabel)
+    plt.legend(['real', 'imag'])
+    if fLim:
+        plt.xlim(fLim)
+    plt.grid(True)
+
+    plt.figure(3)
+    plt.plot(fn, Hmag)
+    plt.title('Magnitude spectrum')  # = amplitude
+    plt.ylabel('Voltage')
+    plt.xlabel(flabel)
+    if fLim:
+        plt.xlim(fLim)
+    plt.grid(True)
+
+    plt.figure(4)
+    plt.plot(fn, Hphs)
+    plt.title('Phase spectrum (note -1: pi = -pi)')
+    plt.ylabel('Phase [rad]')
+    plt.xlabel(flabel)
+    if fLim:
+        plt.xlim(fLim)
+    plt.grid(True)
+
+
+def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None, showRoll=False):
+    """Plot power spectrum
+
+    Input:
+    . f: normalized frequency axis for HF (fs = 1)
+    . HF: spectrum, e.g. frequency transfer function HF = DTFT(h)
+    . fmt: curve format string
+    . fs: sample frequency in Hz, scale f by fs, fs >= 1
+    """
+    if fs > 1:
+        flabel = 'Frequency [fs / %d]' % fs
+    else:
+        flabel = 'Frequency [fs]'
+
+    plt.plot(f * fs, pow_db(HF), fmt)
+    plt.title('Power spectrum')
+    plt.ylabel('Power [dB]')
+    plt.xlabel(flabel)
+    if fLim:
+        plt.xlim(fLim)
+    if dbLim:
+        plt.ylim(dbLim)
+    plt.grid(True)
+
+
+def plot_two_power_spectra(f1, HF1, name1, f2, HF2, name2, fs=1.0, fLim=None, dbLim=None, showRoll=False):
+    """Plot two power spectra in same plot for comparison
+
+    Input:
+    . f1,f2: normalized frequency axis for HF1, HF2 (fs = 1)
+    . HF1, HF2: spectrum, e.g. frequency transfer function HF = DTFT(h)
+    . fs: sample frequency in Hz, scale f by fs, fs >= 1
+    """
+    if fs > 1:
+        flabel = 'Frequency [fs / %d]' % fs
+    else:
+        flabel = 'Frequency [fs]'
+
+    if showRoll:
+        plt.plot(f1 * fs, pow_db(HF1), 'r',
+                 f1 * fs, np.roll(pow_db(HF1), len(f1) // 2), 'r--',
+                 f2 * fs, pow_db(HF2), 'b',
+                 f2 * fs, np.roll(pow_db(HF2), len(f2) // 2), 'b--')
+        plt.legend([name1, '', name2, ''])
+    else:
+        plt.plot(f1 * fs, pow_db(HF1), 'r',
+                 f2 * fs, pow_db(HF2), 'b')
+        plt.legend([name1, name2])
+    plt.title('Power spectrum')
+    plt.ylabel('Power [dB]')
+    plt.xlabel(flabel)
+    if fLim:
+        plt.xlim(fLim)
+    if dbLim:
+        plt.ylim(dbLim)
+    plt.grid(True)

applications/lofar2/model/pfb_os/dsp_study_erko.txt

+###############################################################################
+# Copyright 2023
+# ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+###############################################################################
+#
+# Author: Eric Kooistra
+# Purpose: DSP theory summary
+#
+# References:
+#
+# * [LYONS] Understanding Digital Signal Processing, 3rd edition
+# * [PROAKIS] Digital Signal Processing, 3rd edition
+# * [HARRIS] Multirate Signal Processing for Communication Systems
+# * [CROCHIERE] Multirate Signal Processing
+# * [JOS1] Mathematics of the Discrete Fourier Transform
+# * [JOS2] Introduction to Digital Filters
+# * [JOS3] Physical Audio Signal Processing
+# * [JOS4] Spectral Audio Signal Processing
+#
+
+1) Linear Time Invariant (LTI) system [LYONS 1.6]
+
+- Time shift in input causes equal time shift in output.
+- Order of sequential LTI operations can be rearranged with no change in final
+  output
+- LTI output is linear combination (weighted sum) of delayed copies of the
+  input signal [JOS1 4.3.11]
+- Transforms use discrete or continuous range of complex phasors [JOS1 4.3.12]
+  * For causal signals x(n) = 0 for n < 0, yields unilateral (one sided)
+    transforms. Used to solve difference equations with initial conditions
+    [PROAKIS 3.5]. Only unique for causal signals, because these are 0 for
+    n < 0.
+  * DFT: Every signal x(n) can be expressed as a linear combination of complex
+    sinusoids W_N^kn = exp(j w_k t_n). The coefficients of projecting x(n) on
+    W_N^kn for n = 0,1,...,N-1 yield the DFT of x is X(k) for k = 0,1,...,N-1.
+  * DTFT: For N --> inf, linear combination of exp(j w t_n) = exp(j w T)^n
+  * z-transform: linear combination of z^n. Generalization (analytic
+    continuation) of the DTFT on the unit circle in the complex plane, to the
+    entire complex z-plane.
+  * FT: integral 0 --> inf, linear combination of exp(j w t)
+  * Laplace transform: integral 0 --> inf, linear combination of exp(s t),
+    s = o + jw, so analytic continuation of FT (on jw axis) to the s-plane.
+- Analogue (Laplace) and digital (z) complex planes [JOS1 4.3.13]
+    . Transform of growing functions, important for transient behaviour and
+      stability analysis
+    . Poles and zeros, frequency reponse z = exp(jw) [PROAKIS 4.4.6]
+    . Laplace: Every point in s-plane corresponds to a complex sinusoid
+      A exp(s t), t >= 0. Frequency axis s = j w
+    . Z: Every point in z-plane corresponds to sampled z = exp(s T), z^n
+
+- FT, DTFT, DFT [MATLAB sinusoids-and-fft-frequency-bins.m]
+  . Fourier transform (FT): continuous time <-> continuous frequency
+  . Discrete-time Fourier transform (DTFT): discrete time <-> continuous
+    frequency. If there is no aliasing, then the DTFT is the same as the
+    Fourier transform up to half the sampling frequency.
+  . Discrete Fourier transform (DFT): discrete time <-> discrete frequency.
+    For a signal that is nonzero only in the interval 0 <= n < N, an Ndtft
+    point DFT exactly samples the DTFT when N_dtft >= N.
+
+- Transforms are used because the time-domain mathematical models of systems
+  are generally complex differential equations. Transforming these complex
+  differential equations into simpler algebraic expressions makes them much
+  easier to solve. Once the solution to the algebraic expression is found,
+  the inverse transform will give you the time-domain repsponse.
+  . Fourier is a subset of Laplace. Laplace is a more generalized transform.
+    Fourier is used primarily for steady state signal analysis, while Laplace
+    is used for transient signal analysis. Laplace is good at looking for the
+    response to pulses, step functions, delta functions, while Fourier is
+    good for continuous signals.
+    Many of the explanations just mention that the relationship between the
+    Laplace and Fourier transforms is that s = o + jw, so the Fourier
+    transform becomes a special case of the laplace transform. Better
+    explanations deals that Laplace is used for stability studies and Fourier
+    is used for sinusoidal responses of systems. Systems are stable if the
+    real part of s is negative, that is to say there is a transient that will
+    vanish in time, in those cases, it is enough to use Fourier. Of course
+    you will lose the insight of the transient part. Laplace should be able to
+    determine the full response of a system, be it stable or unstable,
+    including transient parts.
+
+2) Windows [JOS4 3]
+- Tabel [PROAKIS 8.2.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
+- Properties of rectangular window with M points from [JOS4 3.1.2]:
+  . Zero crossings at integer multiples of 2 pi / M = Ndtft / M [LYONS Eq.
+    3.45]
+  . Main lobe width is 4 pi / M
+  . 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.
+  . First side lobe only 13 dB down from the main-lobe peak.
+  . Side lobes roll off at approximately 6dB per octave.
+  . A phase term arises when we shift the window to make it causal, while
+    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
+  . Pass band ripple (allowed gain deviation):
+      r_pass = 10**(r_pass_dB / 20) - 1
+  . Stop band edge frequency: w_stop
+  . Stop band ripple (allowed leakage level):
+      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
+- 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)
+    - 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
+  . LS and Remez can do bandpass (= flat), differentiator, hilbert
+  . Linear phase filter types (filter order is Ntaps - 1, fNyquist = f2/2):
+      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
+
+2) Finite Impulse Response (FIR) filters
+- FIR filters perform time domain Convolution by summing products of shifted
+  input samples and a sequence of filter coefficients [LYONS 5.2].
+- Convolution equation [LYONS Eq. 5.6]:
+
+         N-1
+  y(n) = sum h(k)(x(n-k) = h(k) * x(n)
+         k=0
+
+- Impulse response h(k) are the FIR filter coefficients:
+
+  x(n) --> x(n-1) --> ... --> x(N-1) --\
+        |          |                   |
+       h(0)       h(1)                h(N-1)
+        \----------\-- ... ------------\--> + --> 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)
+
+- 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 ~= 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
+
+- Linear phase FIR filter
+  . Even or odd symmetrical h(n) = +-h(M - 1 - n), n = 0,1,...,N-1
+    [PROAKIS 8.2.1]. Reason for using FIR [LYONS 5.10.3]
+  . Group delay (= envelope delay[LYONS 5.8]) of symmetrical FIR filter is:
+      G = (N_taps - 1) / 2 [Ts]  [LYONS 5.10.3]
+  . 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
+  . When Ntaps is odd, then half of the coefs are 0.
+
+
+3) 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:
+
+  W_N = exp(-j 2pi / N) is primitive Nth root of unity
+  W_N^k = exp(-j 2pi k / N)
+  W_N^kn = exp(-j 2pi k / N * n) = exp(-j w_k * t_n)
+  . w_k = k 2pi fs / N
+  . t_n = n Ts
+  . fs / N is the frequency sampling interval in Hz
+  . Ts = 1 / fs is time sampling interval in seconds
+  . ws = 2pi fs is the frequency sampling interval in rad/s
+
+- Normalized frequency axis [LYONS 3.5, 3.13.4]:
+  . fs = 1, ws = 2pi fs = 2pi
+  .  -fs/2, 0, fs/2 [Hz]
+       -pi, 0, pi  [rad]
+      -0.5, 0, 0.5 [fs]
+  . N even, e.g. N = 4:
+
+         <------ N = 4 ------->
+         0           fs/2           ( fs  )
+         |             |            ( |   )
+     n = 0    1        2    3       ( |   )
+         0/4  1/4      2/4  3/4     ( 4/4 )
+         DC   positive |    negative
+                       |
+                       \--> fftshift([0, 1, 2, 3]) = [2, 3, 0, 1]
+
+  . N odd, e.g. N = 5:
+         <------- N = 5 -------->
+         0           fs/2          ( fs  )
+         |             |           ( |   )
+     n = 0    1    2   | 3    4    ( |   )
+         0/5  1/5  2/5 | 3/5  4/5  ( 5/5 )
+         DC   positive | negative
+                       |
+                       \-->fftshift([0, 1, 2, 3, 5]) = [3, 4, 0, 1, 2]
+
+  . With K = N // 2:
+    . N even : DC, K - 1 positive, fs/2, K - 1 negative frequencies
+    . N odd  : DC, K positive, K negative frequencies
+
+- The DFT project a length N signal x on a set of N sampled complex sinusoids
+  that are generated by the Nth roots of unity [JOS4 6.1]. These sinusoids form
+  an orthogonal basis and are the only frequencies that have a whole number of
+  periods in N samples [JOS1 6.2].
+
+- Discrete Fourier Transform of x(n), n, k = 0,1,...,N-1 [JOS1 5.1, 6.6, 7.1],
+  [PROAKIS 5.1.2, 5.1.3]:
+
+                  N-1
+  X(w_k) = X(k) = sum x(n) W_N^kn
+                  n=0      exp(-j w_k t_n)
+                           exp(-j 2pi k n / N)
+  Inverse DFT:
+
+             N-1
+  x(n) = 1/N sum X(k) exp(+j w_k t_n)
+             k=0      s_k(n)
+    with:
+    . s_k(n) = exp(+j w_k t_n)
+
+- Matrix formulation, DFT as linear transformation [JOS1, PROAKIS 5.1.3]:
+
+  DFT:
+  XN       = WN                                            xN
+
+  |X(0)  |   |1  1     1              ... 1             | |x(0)  |
+  |X(1)  |   |1  W_N   W_N^2          ... W_N^(N-1)     | |x(1)  |
+  |X(2)  | = |1  W_N^2 W_N^4          ... W_N^2(N-1)    | |x(2)  |
+  |...   |   |...                     ...               | |...   |
+  |X(N-1)|   |1  W_N^(N-1) W_N^2(N-1) ... W_N^(N-1)(N-1)| |x(N-1)|
+
+  IDFT:
+  xN = WN^-1 XN = 1/N conj(WN) XN, so
+
+  WN conj(WN) = N IN, where IN i N x N identity matrix
+
+- Real input: X(k) = conj(X(N - k))
+
+- Spectral leakage or cross-talk occurs for frequencies other then w_k, because
+  this cause truncation distortian in the periodic extension of x(n) =
+   x(n + mN) [JOS1 7.1.2].
+
+- DTFT Ndtft >= N [LYONS 3.11, JOS1 7.2.9]
+  . Zero padding N to Ndtft increases the interpolation density, but does not
+    increase the frequency resolution in the ability to resolve, to distinguish
+    between closely space features in the spectrum. The frequency resolution
+    can only be increased by increasing N
+- N point DFT of rectangular function yields aliased sinc (= Dirchlet kernel):
+
+  x(n) = 1 for K samples
+  X(m) = c sin(pi * m * K / N) / sin(pi * m / N)
+
+  . c = 1 for symmetric rect -(K-1)/2 to + (K-1)/2
+  . m = 0: X(0) = K
+  . m first zero crossing = N / K --> main lobe width is 2N / K
+  . 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
+
+4) 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,
+  so there are D or U phases [LYONS 10.7]. The LPF with all phases is called
+  the protype 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]
+
+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.
+
+Upsampling + LPF = interpolation:
+- 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]
+
+
+5) 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
+    extension
+- FLIP_n(x) = x(-n) reversal
+  . x(0) = x(-0)
+  . x(-n) = x(N - n) for modulo N
+- SHIFT_L,n(x) = x(n - L)
+
+- ZEROPAD_M,m(x) = x(n) for |m| < N / 2
+                 = 0 else
+  . zero centered, zero-phase, periodic
+  . ZEROPAD_10([1,2,3,4]) = [1,2,0,0,0,0,0,0,3,4]
+  . ZEROPAD_10([1,2,3,4,5]) = [1,2,3,0,0,0,0,0,4,5]
+- CAUSALZEROPAD_M,m(x) = x(n) for m < N
+                       = 0 for N < m < M
+  . CAUSALZEROPAD_10([1,2,3,4]) = [1,2,3,4,0,0,0,0,0,0]
+- INTERP_L,k'(X) = X(w_k'), for X(w_k) with k = 0,1,...,N-1, where:
+    w_k' = 2pi k' / M, k' = 0,1,...,M-1
+    M = L * N
+  . Interpolates a signal by factor L using bandlimited interpolation
+
+- STRETCH_L,m(x) = x(n) for n = m / L is an integer
+                 = 0 else
+  . Upsampling by inserting L-1 zeros after every x sample
+- REPEAT_L,m(x) = x(m), where:
+    m = 0,1,...,M-1 and M = L * N and n = m modulo N
+  . REPEAT_2([1,2,3,4]) = [1,2,3,4,1,2,3,4]
+- DOWNSAMPLE_L,m = x(mL), where N = L * M and m = 0,1,...,M-1
+  . DOWNSAMPLE_L(STRETCH_L(x) = x
+  . DOWNSAMPLE_2([1,2,3,4,5,6]) = [1,3,5]
+                 L-1
+- ALIAS_L,m(x) = sum x(m + lM), with m = 0,1,...,M-1 and N = L*M
+                 l=0
+  . ALIAS_3([1,2,3,4,5,6] = [1,2] + [3,4] + [5,6] = [9,12]
+
+- DFT_N,k(x) = X(k) with N points
+- DFT relations:
+  . STRETCH_L <==> REPEAT_L
+  . DOWNSAMPLE_L(x) <==> 1/L * ALIAS_L(X)  [JOS1 7.4.11]
+  . ZEROPAD_LN(x) <==> INTERP_L(X)  [JOS1 7.4.12]
+  . conj(x) <==> FLIP(conj(X))
+  . FLIP(conj(x)) <==> conj(X)
+  . FLIP(x) <==> FLIP(X)
+  . x even <==> X even
+  . SHIFT_L(x) <==> exp(-j w_k L) X(k)
+
+
+
+https://learning.anaconda.cloud/
+https://realpython.com/python-scipy-fft/
+
+https://mne.tools/0.24/auto_tutorials/preprocessing/25_background_filtering.html