diff --git a/applications/lofar2/model/pfb_os/dsp_study_erko.txt b/applications/lofar2/model/pfb_os/dsp_study_erko.txt
index 15d3ef6b3bf0b3d602595bcf695c0bf9b737b1dd..62305fe862a04e1634b6f91943bb8b38774743c1 100644
--- a/applications/lofar2/model/pfb_os/dsp_study_erko.txt
+++ b/applications/lofar2/model/pfb_os/dsp_study_erko.txt
@@ -47,10 +47,11 @@
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)
+ [LYONS 3.14]
+ * z-transform: sum n=0 --> inf, 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) [LYONS 3]
* 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]
@@ -61,14 +62,17 @@
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]
+- FT, DTFT, DFT [LYONS 3.14, 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.
+ - The DFT is equal to the DTFT of a periodic time-domain discrete sequence,
+ so the DFT 'assumes' its input sequence is periodic in time [LYONS 3.14]
+ - 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. For N_dtft > N
+ the DFT approximates the DTFT [LYONS 3.11].
- Transforms are used because the time-domain mathematical models of systems
are generally complex differential equations. Transforming these complex
@@ -179,18 +183,22 @@
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]:
+- Convolution equation [LYONS Eq. 5.6, JOS4 7.2.4]:
N-1
- y(n) = sum h(k)(x(n-k) = h(k) * x(n)
+ 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)
+ x(n) --> x(n-1) --> ... --> x(n-(N-1)) --\
+ | | |
+ h(0) h(1) h(N-1)
+ \----------\-- ... ----------------\--> + --> y(n)
+
+- Difference equation, so b(m) = h(m) for m = 0, 1, ...,N-1:
+
+ y[n] = b[0]x[n] + b[1]x[n-1] + ... + b[N-1] x[n-(N-1)]
- Convolution in time domain is equivalent to multiplication in frequency
domain:
@@ -212,7 +220,7 @@
- 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:
+ the discontinuity in time domain [JOS4 B.18]:
. 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.
@@ -226,6 +234,9 @@
. Design using windowed sinc, because windows are symmetrical [PROAKIS 8.2.2,
LYONS 5.3.2, DSPGUIDE 16]
+- Zero phase signal is a linear phase signal with phase slope 0 and phase 0
+ or pi. It is a zero-centered signal or even signal [JOS1 7.4.4].
+
- Half band FIR filter [LYONS 5.7, HARRIS 8]
. When Ntaps is odd, then half of the coefs are 0.
- Ntaps + 1 is multiple of 4
@@ -238,11 +249,53 @@
. 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_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)
+- Filtering using FFT [LYONS 13.10, PROAKIS 5.3.2, JOS4 8.2.5]
+ . Convolution in time domain <==> multiplication in frequency domain (and
+ vice versa) [LYONS 5.9.2].
+ . rectangular window has sinc transfer function in frequency domain, so
+ applying (multiplying) window in time domain convolves in frequency
+ domain, causing the main pass band lobe to become wider (and the stop
+ band side lobes to become lower = purpose of window).
+ . sampling in time domain cause spectrum to repeat in frequency domain,
+ because sampling is a train of delta pulses with time period Ts and
+ their FT is a train of delta pulses with frequency interval fs
+ [JOS4 B.16].
+ . For L_h >~ 64 filtering in frequency domain is more efficient than
+ convolution [JOS4 8.1.4]
+ . Use overlap save (OLS) method (prepend block with zeros) or overlap add
+ (OLA) method (postpend block with zeros) of input x to fit FFT size. The
+ FFT size is Lx + L_h - 1 to avoid aliasing [JOS4 8.1.2, 8.2.6], where L_x
+ is length of block from x. Use IFFT on X_block(f) H(f). For OLS discard
+ first samples and for OLA add last sample to first samples of next.
+ . Without output overlap (so FFT size N = Lx and n - k index is modulo N)
+ the frequency domain convolution is cyclic. By adding input zeros, the FFT
+ convolution becomes effectively acyclic [JOS4 8.1]. The added zeros cause
+ that the output blocks overlap, to keep the output and input sample rate
+ the same.
+
+- Correlation equation [JOS1 7.2.5] (for comparison with convolution):
+
+ N-1
+ xy(n) = sum conj(x(k)) y(n + k), time shift n is correlation lag
+ k=0
+
+ The cross power spectrum is [JOS1 8.4]:
+
+ HXY(w_k) = DFT_k(xy(n)) = 1 / N conj(X(w_k)) Y(w_k)
+
+- FIR system identification from input-output measurements [JOS1 8.4.5,
+ PROAKIS 12]
+
+ y = h * x <==> H Y
+
+ xy = x cross y <==> conj(X) Y = conj(X) H Y = H |X|^2, so H = Rxy / Rxx
+
4) Hilbert transform (HT) and analytic signal [LYONS 9]
- The Hilbert transform (HT) creates a 90 degrees phase-shifted version of a
@@ -269,7 +322,7 @@
- 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
+ . HT window design: ht[n] * w[n], best use Ntaps is odd for unit group delay.
. 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
@@ -285,15 +338,127 @@
. 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).
+ . Using the DFT [LYONS 9.4.2], like with scipy.signal.hilbert. The
+ scipy.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]
+ . HT design using two path, all pass IIR [HARRIS 10.7, STREAM 9]
-- IIR Hilbert transform [HARRIS 10.7]
-5) Discrete Fourier Transform (DFT)
+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
+- 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
+ they become +ak in H(z), so with same sign as a0. The scipy.signal.freqz
+ also uses -ak for y[n] to have +ak for H(z).
+
+ a0 y[n] = b0 x[n] + b1 x[n-1] + ... + bM x[n-M]
+ - a1 y[n-1] - ... - aN y[n-M]
+
+ M N
+ = sum bk x[n - k] - sum ak y[n - k]
+ k=0 k=1
+
+- Implmentation structure [LYONS 6.6, JOS2 9.1]:
+ . Direct-Form I: B(z) * 1 / A(z)
+
+ x[n] ---v-----> b0 --> + --> + -------------v----> y[n]
+ z^-1 --> b1 --> + + <-- -a1 <-- z^-1
+ z^-1 --> b2 --> + + <-- -a2 <-- z^-1
+
+ . Direct-Form II: 1 / A(z) * B(z), to share z^-1 delay memories, but needs
+ more bits to avoid internal overflow.
+
+ x[n] --> + -------------v-----> b0 --> + --> y[n]
+ + <-- -a1 <-- z^-1 --> b1 --> +
+ + <-- -a2 <-- z^-1 --> b2 --> +
+
+ . Transposed Direct-Form II: reverse all arrows, swap nodes and +, swap
+ input and output. Has advantage that zeros are calculated first, so that
+ problematic large intermediate values from poles are avoided.
+
+ x[n] --v-> b0 -->--+------------v--> y[n]
+ --> b1 --> z^-1 <-- a1 <--
+ --> b2 --> z^-1 <-- a2 <--
+
+- For FIR b = h. For IIR it is not possible to directly derive b, a from h
+ [LYONS 6.1]. Therefor use z-transform [LYONS 6.3]:
+
+ +inf +inf
+ H(z) = sum h[n] z^-n = sum h[n] r^-n exp(-j w n)
+ n=-inf n=-inf
+
+ = DFT(h[n]) for |z| = 1, so on unit circle
+
+ with polar form z = r exp(j w), w = -pi to +pi, w = 0 at z = 1 + 0j.
+
+- Unit delay y[n] = x[n-1] --> Y(z) = z^-1 X(z). Therefore IIR z-transform
+ of difference equation is [LYONS 6.4.1 uses -ak, JOS2 6.5 uses +ak in H(z)]:
+
+ Y(z) M N
+ H(z) = ---- = sum bk z^-k / (a0 + sum ak z^-k)
+ X(z) k=0 k=1
+
+ . 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):
+
+ Y(z) b0 + b1 z^-1 + b1 z^-2 (z - z0)(z - z1)
+ H(z) = ---- = ---------------------- = K ----------------
+ X(z) a0 + a1 z^-1 + a2 z^-2 (p - p0)(p - p1)
+
+. Transfer function to:
+ . scipy.signal.tf2sos(b, a)
+ . scipy.signal.tf2zpk(b, a)
+- Implementation issues for finite word length: coefficient quantization,
+ 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).
+ . 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]
+
+- IIR filter design methods:
+ . Impulse invariance IIR design methods, sampled version of an analogue
+ filter impulse response [LYONS 6.10]. Suffers from aliasing, when fs is
+ not high enough. This typically shows as nonzero frequency response value
+ at fs / 2. Two methods that yield the same:
+ - Using inverse Laplace Transform and the z transform,
+ - Partial fraction expansion into multiple single pole parts.
+ . 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].
+
+ 1 - z^-1
+ s = 2 / Ts --------
+ 1 + z^-1
+
+ w_d = 2 arctan(w_a Ts / 2), with -pi <= w_d <= +pi,
+ f_d = w_d / (2 pi) fs
+
+ . scipy.signal.bilinear(b, a, fs)
+ scipy.signal.bilinear(z, p, k, fs)
+
+ . 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).
+ . 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.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).
+
+6) 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:
@@ -342,8 +507,8 @@
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]:
+- Discrete Fourier Transform of x(n), n, k = 0,1,...,N-1 [LYONS 3, 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
@@ -397,7 +562,123 @@
X(m) = sin(pi * m) / sin(pi * m / K)
~= K * sinc(m) for K = N >~ 10
-6) Multirate processing:
+- Fourier transform theorems [JOS4 B]
+ . Scaling: x(t / a) <==> |a| X(a w)
+ . Shift: x(t - T) <==> exp(-j w T) X(w)
+ . Modulation: x(t) exp(j v t) <==> X(w - v), is dual of shift
+ . Convolution:
+ x * y <==> X Y
+ x y <==> 1 / (2 pi) X * Y
+ . flip(x) <==> flip(X)
+ . d(t) <==> 1, dirac pulse with area 1 at t = 0
+
+ +inf
+ . d_train_P(t) = sum d(t - m P), period P
+ m=-inf
+ <==>
+ +inf
+ d_train_P(f) = 1 / P sum d(f - m / P)
+ m=-inf
+
+ . sampling: x_d(t) = x(t) d_train_Ts(t)
+ <==>
+ X_d(f) = X * d_train_fs(f)
+ +inf
+ = fs sum X(f - k fs)
+ k=-inf
+
+
+7) 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.
+- The m-th windowed data frame of x(n) is [JOS4 8.2.1]:
+
+ x_m(n) = x(n) w(n - m R), n in Z, m is frame index, R is hop size
+ x_m = x SHIFT_mR(w)
+
+ STFT is the DTFT of x_m [JOS4 8], w = w_k [JOS4 9]:
+
+ +inf
+ X_m(w) = sum x[n] w[n - m R] exp(-j w n) = DTFT_w(x SHIFT_mR(w))
+ n=-inf
+
+- STFT spectrogram parameters [JOS4 7.2]:
+ . Window size W determines frequency resolution
+ . FFT size N >= W determines the frequency grid (amount of spectral
+ interpolation beyond the resolution)
+ . Hop size R is number of samples between successive windows, determines
+ resolution in time, R <= W. R = 1 is 'sliding FFT'. R > 1 is downsampling
+ in time. N > W and R < W cause that the output sample rate increases.
+
+- OLA = Overlap Add and uses rectangular window w of length Lx on x and step
+ with Lx through x. The purpose of OLA is to use FFT and implement
+ convolution with impulse response h as multiplication of X(f) H(f) in
+ frequency domain, and then do IFFT [JOS4 8.2].
+- COLA = Constant Overlap Add is like OLA, but applies a window w on x (and
+ no filter h). The COLA constraint for perfect reconstuction of x(n) from
+ x_m(N), and therefore also from IFFT thanks to linearity, is:
+
+ +inf +inf
+ x(n) = sum x_m(n) ==> sum w(n - m R) = 1, for all n in Z
+ m=-inf m=-inf
+
+ . COLA(R) means that the window has COLA for hop size R [JOS4 8.2.1, 8.2.2]:
+ - Any window is COLA(1).
+ - Retangular window is COLA(W), but also COLA(W / k).
+ - Hamming window is COLA(W / 2) and COLA(W / 4)
+ . With COLA the sum X_m(w) = X(w) = DTFT_w(x)
+ . With filtering, then STFT is used like OLA to perform spectral
+ modifications via H(f).
+ . WOLA = Weighted OLA implies using a two windows, one before FFT (like with
+ OLA) and one after IFFT. Typically both windows are equal and the sqrt(w),
+ so that together they are COLA [JOS4 8.6].
+ . COLA contraint in frequency domain [JOS4 8.3.2]. Window transform is zero
+ for all harmonics of the frame rate 2 pi / R:
+
+ W(w_k) = 0, with w_k = 2 pi k / R and for |k| = 1,2,...,R-1
+
+ This is the weak COLA contraint. The strong COLA constraint, that is
+ useful if the signal is modified in the frequency domain, requires that
+ the window is bandlimitted consistent with down sampling by R:
+
+ W(w) = 0, for |w| >= pi / R
+
+ . Convolution of x with d() leaves the signal x unchanged. Convolution of
+ x with a periodic train of d() pulses creates aliases of x:
+
+ w_train_R(n) = w(n) * d_train_R(n) --> aliases of w
+
+ +inf
+ = sum w(n) * d(n - m R), period R
+ m=-inf
+
+ +inf
+ = sum w(n - m R), period R
+ m=-inf
+ <==> ???
+ w_train_R(f) = W(f) d_train_R(f) --> sampling of W
+
+ +inf
+ = W(f) 1 / R sum d(f - m / R), for f = m / R
+ m=-inf
+
+ +inf
+ = 1 / R sum W(m / R), for f = m / R
+ m=-inf
+
+ . Filterbank interpretation of STFT is using flip(w) as impulse response
+ and filter down converted x by -w_k (frequency shift by modulation)
+ [JOS4 9.1.2, HARRIS 6]:
+
+ +inf
+ X_m(w_k) = sum x[n] exp(-j w_k n) w[n - m R]
+ n=-inf
+
+ = [x_k * flip(w)](m), with x_k(n) = x(n) exp(-j w_k n)
+
+
+8) 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,
@@ -420,7 +701,7 @@ Upsampling + LPF = interpolation:
pass band of sin(x)/x [LYONS 10.5.1]
-Appendix 7) Signal operators [JOS1 7.2]
+Appendix 8) 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
@@ -453,7 +734,7 @@ Appendix 7) Signal operators [JOS1 7.2]
. 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
+- ALIAS_L,m(x) = sum x(m + l M), 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]
@@ -474,3 +755,12 @@ https://learning.anaconda.cloud/
https://realpython.com/python-scipy-fft/
https://mne.tools/0.24/auto_tutorials/preprocessing/25_background_filtering.html
+
+Fixed-point versus floating point in antenna data
+
+I do not think having an exponent helps here, because the information is in the mantissa. Suppose you have system noise (sky noise + receiver noise) and strong RFI. In quantisation you want to keep the quantisation noise small with respect to the system noise. This then means that the mantissa needs to fit the RFI and preserve the level system noise. If you would reduce the mantissa width and increase the exponent, then the quantisation noise becomes dominant with respect to the system noise. The point with RFI is that you hope that if you separate the subbands into finer channel, that then some fine channels are RFI free. However, if you have introduced more quantisation noise due to using the exponent instead of a wider mantissa, then you introduced more white quantization noise in the whole subband, so then separating into finer channel does not yield clean fine channels anymore.
+
+For subbands without RFI we need e.g. 4 bit complex subband samples, provided that they are equalized to their nominal level (i.e. compensated for frequency dependent gain of the recever chain).
+For subbands with RFI we need enough extra bits to fit the RFI level, so that we can still hope to find finer channels without RFI within that subband.
+This implies that we would need extra meta information that tells number of bits per subband. Alternative is to use one fixed width that kind of fits all subbands, e.g. 8 bit complex subband values.
+Anyway, I think in each case we would transport fixed point (= integer) values, because transporting exponent bits does not add information.