diff --git a/applications/lofar2/model/pfb_os/dsp_study_erko.txt b/applications/lofar2/model/pfb_os/dsp_study_erko.txt
index 8764090c5f10f2efcb576ae4118041dce674f418..b414cff3692cc3b38e52b85f096ead5463dcb0b2 100644
--- a/applications/lofar2/model/pfb_os/dsp_study_erko.txt
+++ b/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