More on IIR.

applications/lofar2/model/pfb_os/dsp_study_erko.txt

@@ -21,14 +21,24 @@
 #
 # 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
+# * [LYONS] Understanding Digital Signal Processing, 3rd edition 2011
+# * [PROAKIS] Digital Signal Processing, 3rd edition 1996
+# * [HARRIS] Multirate Signal Processing for Communication Systems, 2004,
+# * [STREAMLINING] Streamlining DSP, 2007, 2nd edition 2012
+# * [CROCHIERE] Multirate Signal Processing, 1983
+# * [VAIDYANATHAN] Multirate Systems and Filterbanks, 1993
+# * [JOS1] Mathematics of the Discrete Fourier Transform, 2008
+# * [JOS2] Introduction to Digital Filters, 2007
+# * [JOS3] Physical Audio Signal Processing, 2010
+# * [JOS4] Spectral Audio Signal Processing, 2011
+#
+# * [WIKI] https://en.wikipedia.org/wiki/Bilinear_transform
+# * [CHIPMUNK] https://github.com/chipmuenk : Python Frequency Design Analysis and DSP
+# * [WHDLWHIZ] https://vhdlwhiz.com/articles/ : FIR filter design using DSP blocks
+# * [BIQUAD]
+#   . https://www.w3.org/TR/audio-eq-cookbook/
+#   . https://webaudio.github.io/Audio-EQ-Cookbook/Audio-EQ-Cookbook.txt
+#   . Configure the Coefficients for Digital Biquad Filters in TLV320AIC3xxx Family (pdf)
 #
 
 1) Linear Time Invariant (LTI) system [LYONS 1.6]
@@ -154,7 +164,7 @@
       . fircls, fircls1 = constrained ls with pass, stop ripple
     - The Remez or Parks-McClellan algorithm (scipy.signal.remez(), MATLAB
       firpm)
-  . MATLAB filters specify order n, so yield n + 1 coefs
+  . MATLAB FIR filters specify order n, so yield n + 1 coefs
   . LS and Remez can do bandpass (= flat), differentiator, hilbert
   . LS and Remez for large Ntaps (> 1000) can fail, but can be achieved
     using Fourier interpolation a filter with less coefs.
@@ -217,7 +227,7 @@
     - 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
+    - Choose transition region specification in order of 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 [JOS4 B.18]:
@@ -344,12 +354,26 @@
     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]
-  . HT design using two path, all pass IIR [HARRIS 10.7, STREAM 9]
+  . An Efficient Analytic Signal Generator [STREAMLINING 38]. Use -45 and +45
+    degrees phase rotation of a complex bandpass filter that is 0 for f < 0
+    and has rect window with cosine-taper (~= Tukey window, ~= raised cosine)
+    for f > 0. This makes that the real branch filter and the imag branch
+    filter have time - reversed coefficients, which implies that their
+    magnitude responses are identical and the filters can share coefficients.
+    However the filters are no longer symmetrical, due to the phase rotation,
+    so they do not have linear phase anymore, but within the band they have
+    almost linear phase.
+    The filter has taper a at f_low and f_high, so it can be defined for any
+    bandpass. This is useful to combine the analytic signal creation with a
+    specific BPF.
+  . HT design using two path, all pass IIR [HARRIS 10.7, STREAMLINING 9]
 
 
 5) Infinite Impulse Response (IIR) filters and z-transform
 - IIR filters are called recursive and FIR filters are called nonrecursive,
-  however FIR filters can also be recursive
+  however FIR filters can also be recursive [LYONS 6.6.1]
+- IIR filter order is largest of nof zeros or nof poles, so largest of len(b)
+  - 1 or len(a) - 1.
 - Difference equation:
   . Sign of ak coefficients, [LYONS 6.1 uses +ak], [JOS2 5.1, PROAKIS 7.1 use
     -ak] in y[n]. Using -ak in difference equation is preferred, because then
@@ -406,7 +430,7 @@
   . scipy.signal.freqz plots frequency response along unit circle
     z = exp(j w) defines ak for k > 1 as negative
 
-- biquad (= 2nd-order section = sos):
+- biquad (= second-order section = sos):
 
              Y(z)   b0 + b1 z^-1 + b1 z^-2     (z - z0)(z - z1)
       H(z) = ---- = ---------------------- = K ----------------
@@ -419,7 +443,7 @@
   overflow, rounding [LYONS 6.7].
   . Transposed Direct-Form II is preferred, but with Direct-Form I the z^-1
     delays can also be shared between biquad stages [LYONS 6.8.1].
-  . Use biquads (sos).
+  . Use biquads = second-order sections (sos).
   . Do biquad with poles closest to unit circle first, combined with zeros
     that are closest to those poles. Order biquads in either decreasing or
     increasing pole distance fron the unit circle [LYONS 6.8.1]
@@ -434,7 +458,8 @@
       . scipy.signal.residuez(b, a, fs)
   . Bilinear transform IIR design method, maps s-plane to z-plane, this warps
     the analogue frequencies axis to range 0 - fs to avoid aliasing, but does
-    introduce nonlinear distortion [LYONS 6.11].
+    introduce nonlinear distortion [LYONS 6.11]. This is first order
+    approximation of z = exp(sT) --> s = 1/T ln(z):
 
                  1 - z^-1
       s = 2 / Ts --------
@@ -449,15 +474,21 @@
    . Optimized IIR design methods, iteratively adjust coefficients to achieve
      a specified arbitrary frequency response (f_pass, f_stop, r_pass_dB,
      r_stop_dB bands).
+     . band types: lowpass, highpass, bandpass, bandstop
      . scipy.signal.iirdesign(), can do all of:
-       . scipy.signal.cheby2(), ripple in stop band
-       . scipy.signal.cheby1(), ripple in pass band
        . scipy.signal.butter(), maximally flat in pass band, decay in stop band
+       . scipy.signal.cheby1(), ripple in pass band
+       . scipy.signal.cheby2(), ripple in stop band
        . scipy.signal.ellip(N, rp, rs, Wn): Elliptic (Cauer) digital and analog
          filter design. As rp approaches 0, the elliptical filter becomes a
          Chebyshev type II filter (cheby2). As rs approaches 0, it becomes a
          Chebyshev type I filter (cheby1). As both approach 0, it becomes a
          Butterworth filter (butter).
+       . Cannot do bessel (= maximal flat group delay), because bilinear
+         transform preserves the amplitude response but not the group delay.
+     . scipy.signal.iirfilter(), can design IIR filter for given order.
+     . For single biquad bandpass alternatively use signal.iirpeak()
+       For single biquad bandstop alternatively use signal.iirnotch()
 
 6) Discrete Fourier Transform (DFT)
 - The N roots of unity [JOS1 3.12, 5.1, PROAKIS 5.1.3, LYONS 4.3]. Note JOS