Add allpass IIR note and WOLFSOUND link.

applications/lofar2/model/pfb_os/dsp_study_erko.txt

@@ -39,7 +39,7 @@
 #   . 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)
-#
+# * [WOLFSOUND] https://thewolfsound.com/
 
 1) Linear Time Invariant (LTI) system [LYONS 1.6]
 
@@ -170,7 +170,7 @@
     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]:
+    [LYONS 9.4.1, VAIDYANATHAN 2.4.2]:
       Type  Ntaps  Symmetry  H(0) H(fs/2)
       I     Odd    Even      any  any  --> LPF, HPF
       II    Even   Even      any  0    --> LPF
@@ -214,6 +214,13 @@
   domain:
     y(n) = h(k) * x(n) ==> DFT ==> Y(m) = H(m) X(m)
 
+  For DFT this is circular convolution. With suffcient zero padding N >= len(the
+  circular convolution can calculate the linear convolution:
+
+           N-1
+    y[n] = sum h(k) x((n - k) % N)
+           k=0
+
 - 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
@@ -224,9 +231,9 @@
   . 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 = func(fs, fpass, fstop, pb ripple, sb ripple) [HARRIS 3.3,
-      LYONS 5.10.5]:
-           ~= fs / df * Atten(dB) / 22, df = abs(fstop - fpass)
+    - Nof taps [HARRIS 3.3, LYONS 5.10.5, VAIDYANATHAN 3.2]:
+      . Ntaps = func(fs, fpass, fstop, pb ripple, sb ripple)
+             ~= fs / df * Atten(dB) / 22, df = abs(fstop - fpass)
     - 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
@@ -299,6 +306,22 @@
 
   HXY(w_k) = DFT_k(xy(n)) = 1 / N conj(X(w_k)) Y(w_k)
 
+  Correlation is a measure of similarity between two function x(k) and y(k),
+  for different shifts (lag) in time. Autocorrelation can show the periodicity
+  of a signal, because they it has similarity for some k > 0. To prove that
+  correlation can be expressed as convolution use a helper function
+  [WOLFSOUND]:
+
+  xh[n] = sum_k x[n + k] h[k],  with sum_k for k = -inf to +inf
+        = sum_k x[-(-n - k)] h[k]
+        = sum_k x1[-n - k] h[k],  with x1[p] = x[-p]
+        = (x1[n] * h[n])[-n], because convolution equation is [LYONS Eq. 5.6,
+             JOS4 7.2.4]: y(n) = sum_k x(n - k) h(k) = x(n) * h(k)
+        = x[-n] * h[n])[-n], again with x[-p] = x1[p], so correlation can be
+             calculated by convolving the time flipped x and then time flip
+             the result. Beware correlation is not commutative, so xh[n] !=
+             hx[n].
+
 - FIR system identification from input-output measurements [JOS1 8.4.5,
   PROAKIS 12]
 
@@ -353,7 +376,8 @@
     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]
+  . Half band [LYONS 13.1, 13.37]
+  . Multirate Hilbert Transformer [CROCHIERE 6.4]
   . 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)
@@ -429,6 +453,8 @@
 
   . scipy.signal.freqz plots frequency response along unit circle
     z = exp(j w) defines ak for k > 1 as negative
+  . analytic continuation means if H0(z) = H1(z) for z = exp(jw), then
+    H0(z) = H1(z) for all z [VAIDYANATHAN 2.4.3].
 
 - biquad (= second-order section = sos):
 
@@ -490,7 +516,70 @@
      . For single biquad bandpass alternatively use signal.iirpeak()
        For single biquad bandstop alternatively use signal.iirnotch()
 
-6) Discrete Fourier Transform (DFT)
+
+6) Allpass filters [VAIDYANATHAN 3.4, STREAMLINING 9, HARRIS 10, JOS3 2.8,
+   WOLFSOUND]
+- H(z) is allpass if |H(exp(jw))| = c for c > 0
+
+  --> H(exp(jw)) = c exp(j phi(w))
+
+- Trivial examples: H(z) = 1, H(z) = +-z^(-K) (FIR delay line)
+- First-order allpass filter (real coefficient a) [WOLFSOUND, STREAMLINING
+  9.1]. Note using a = 0 yields delay z^(-1), hence this can be regarded and
+  used as a generalized delay element. Via a the -pi / 2 phase shift can be
+  obtained at any normalized frequency.
+
+            a + z^(-1)     1 + a z
+    H(z) = ------------- = --------
+            1 + a z^(-1)    z + a
+
+    y[n] = a x[n] + x[n - 1] - a y[n - 1]
+
+  The phase shift is 0 at f = 0, -pi at f = fNyquist = fs / 2 and -pi / 2 at
+  fb the break or center frequency [WOLFSOUND]. The fb follows from the bilear
+  transform of an analogue all pass filter and is defined via the coeficient a:
+
+    a = (tan(pi fb / fs) - 1) / (tan(pi fb / fs) + 1)
+
+- Second-order allpass filter [WOLFSOUND, STREAMLINING 9.1, HARRIS 10.5.1]:
+
+           -c + d (1 - c) z^(-1) + z^(-2)     1 + d (1 - c) z - c z^2
+    H(z) = -------------------------------- = -----------------------
+            1 + d (1 - c) z^(-1) - c z^(-2)   z^2 + d (1 - c) z - c
+
+  The phase shift is now -pi at the fb break or center frequency and controlled
+  by parameter d:
+
+    d = -cos(2 pi fb / fs)
+
+  The parameter c determines the bandwidth (BW) or steepness of the phase slope
+  through fb:
+
+    c = (tan(pi BW / fs) - 1) / (tan(pi BW / fs) + 1)
+
+  Hence fb and BW can be controlled independently via parameters d an c.
+
+- Consider Type I allpass H(z) and Type II allpass H(z^2), so M = 2, and 1st
+  and 2nd order as building blocks for M-path = two-path recursive allpass
+  filters [STREAMLINING 9.1].
+
+
+- Higher order IIR allpass filter, must have transfer function of the form:
+
+    H(z) = +-z^(-K) A~(z) / A(z)
+
+  where A~ is obtained by reversing the polynomial coefficients of A. This
+  implies that the poles and zeros occur in (conjugate) reciprocal pairs.
+
+
+- Allpass decomposition theorem [VAIDYANATHAN 3.6]:
+    H0(z) = (A0(z) + A1(z)) / 2
+  If H0(z) and H1(z) are power complementary, so |H0|^2 + |H1|^2 = c^2 for all
+  w, then:
+    H1(z) = (A0(z) - A1(z)) / 2
+
+
+7) 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:
@@ -620,7 +709,7 @@
                          k=-inf
 
 
-7) Short Term Fourier Transform (STFT) [JOS4 7, 8]
+8) Short Term Fourier Transform (STFT) [JOS4 7, 8]
 - Purpose of the STFT is to display frequency spectra in time (spectrogram). A
   window is used on input x, so that non-linear and/or time varying spectral
   modifications can be performed.
@@ -710,7 +799,7 @@
              = [x_k * flip(w)](m), with x_k(n) = x(n) exp(-j w_k n)
 
 
-8) Multirate processing:
+9) 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,
@@ -732,8 +821,12 @@ Upsampling + LPF = interpolation:
   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
+
 
-Appendix 8) Signal operators [JOS1 7.2]
+Appendix A) 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