-
Eric Kooistra authoredEric Kooistra authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
dsp_study_erko.txt 20.98 KiB
###############################################################################
# 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]
- Window types table [HARRIS 3.2, PROAKIS 8.2.2]
- Kaiser beta and side lobe levels figure [HARRIS 3.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(m) for Ndtft = M >~ 10
- Properties of rectangular window with M points from [JOS4 3.1.2]:
. M determines transition band width, window shape determines stop band
side lobe level [HARRIS 8.2]
. Zero crossings at integer multiples of 2 pi / M = Ndtft / M [LYONS Eq.
3.45]
. Main lobe width is 4 pi / M between zero crossings.
. Main lobe cutoff frequency is about fpass ~= 0.6 / M, where fpass is
positive frequencies -6 dB power BW [EK]
. 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, MATLAB]
. h_ideal(n) = 2 f_cutoff sinc(2 f_cutoff n), n in Z
- f_cutoff is normalized cutoff frequency with fs = 1
- f_cutoff is positive frequencies -6 dB power BW
. cutoff frequency remarks:
- [Scipy firwin] defines f_c relative to fNyquist = fs / 2 instead of fs,
so need to specify f_c = f_cutoff / 2.
- [Scipy firls, remez] define fpass relative to fs, so specify fc = f_cutoff.
- 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 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.
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]:
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
. Note coefs = flip(h) is important for anti-symmetric
- Band Pass Filter (BPF):
. h_bp[k] = h_lp[k] * s_shift[k]
= h_lp[k] * cos(2 pi fc k)
= h_lp[k] * cos(k pi / 2), for half band fc = fs / 4,
series [1, 0, -1, 0]
- High Pass Filter (HPF):
. h_hp[k] = h_lp[k] * s_shift[k]
= h_lp[k] * cos(2 pi fc k)
= h_lp[k] * cos(k pi), for fc = fs / 2,
series [1, -1]
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
. Ntaps must be odd if a passband includes fNyquist = fs / 2 [scipy firwin].
. Kaiser window based design [HARRIS 3.2, Fig. 3.8 for beta]:
- 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 = 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
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:
. 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.
. -1 order no decay implies impulses at end-points of impulse response.
- 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, HARRIS 8]
. When Ntaps is odd, then half of the coefs are 0.
- Ntaps + 1 is multiple of 4
- Gain at fs / 4 is 0.5
. Symmetrical frequency response about fc = fs / 4, so fpass + fstop = fs / 2
. Design:
- window: h_lp[n] * w[n] [HARRIS 8.3]
- remez:
. same weigths for stop and pass band [HARRIS 8.4]
. 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_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)
4) Hilbert transform (HT) and analytic signal [LYONS 9]
- The Hilbert transform (HT) creates a 90 degrees phase-shifted version of a
real signal. The HT can be defined as a filter with impulse response ht,
. x_ht[t] = xi[t] = ht[t] * xr[t].
- Analytical signal: xa = xr + j xi
- E.g.: xr = cos(wt) --> xi = sin(wt) --> xa = exp(jwt)
. HT(w) = j, 0, -j for f < 0, 0, > 0
- Analytic signal is zero for negative frequencies, and has only the
positive frequencies of xr, but with 2x magnitude. Therefore Xa(w) is
also called a one sided spectrum (of Xr(w)).
. ht[t] = IFT(HT(w)), yields:
ht[n] = 0, when n is even
= 2 / (pi n Ts), when n is odd
- Purpose analytic signal:
. measuring instantaneous characteristics (magnitude sqrt(xr^2 + xi^2),
phase atan(xi / xr), frequency d/dt(phase)) of time signals, e.g. for AGC,
envelop, demodulation
. frequency translation.
- Anti-symmetric response:
. Linear phase, with discontinuity of pi at f = 0, as is the purpose of HT
. Type IV, Ntaps even --> HT(0) = 0, HT(fs/2) = any
. Type III, Ntaps odd --> HT(0) = 0, HT(fs/2) = 0 --> so BPF, advantages:
- 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
. 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
- h_lp[k] is Ntaps LPF with two sided BW equal to BPF BW
- h_bp[k] = h_lp[k] * exp(2 pi fcenter / fs (k - D)), mixer
= h_cos[k] + j h_sin[k]
. k = 0 ... Ntaps - 1
. D = (Ntaps - 1) / 2, but same for both I and Q, so even Ntaps is as feasible
. h_cos is symmetrical, h_sin is anti-symmetrical
. if BPF fcenter = fs / 4 then mixer = exp(pi / 2 (k - D)), so then half
of h_cos and h_sin coefs are 0, alternately
. 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).
. Half band [13.1, 13.37]
- IIR Hilbert transform [HARRIS 10.7]
5) 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(m) for K = N >~ 10
6) 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]
Appendix 7) 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