From 4df3ea7a95117c2ed6e3ff2dc3c784f142449380 Mon Sep 17 00:00:00 2001 From: Andrew Belt Date: Tue, 12 Jun 2018 03:27:28 -0400 Subject: [PATCH] Add Impulse responses and Brickwall filter sections to DSP --- DSP.md | 65 +++++++++++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 53 insertions(+), 12 deletions(-) diff --git a/DSP.md b/DSP.md index ad7d8d3..9387d52 100644 --- a/DSP.md +++ b/DSP.md @@ -67,7 +67,7 @@ The signal is integrated for a small fraction of the sample time $1/f_{sr}$ to o Analog-to-digital converters (ADCs) convert a digital value to an amplitude and hold it for a fraction of the sample time. A [reconstruction filter](https://en.wikipedia.org/wiki/Reconstruction_filter) is applied, producing a signal close to the original bandlimited signal. -High-quality ADCs may include digital upsampling before reconstruction . +High-quality ADCs may include digital upsampling before reconstruction. [Dithering](https://en.wikipedia.org/wiki/Dither) may be done but is mostly unnecessary for bit depths higher than 16. Of course, noise may also be introduced in each of these steps. @@ -142,9 +142,22 @@ A linear filter is stable (its [impulse response](https://en.wikipedia.org/wiki/ We should now have all the tools we need to digitally implement any linear analog filter response $H(s)$ and vise-versa. +#### IIR filters + +An infinite impulse response (IIR) filter is a digital filter that implements all possible rational transfer functions. +By multiplying the denominator of the rational $H(z)$ definition above on both sides and applying it to an input $x_k$ and output $y_k$, we obtain +$$\sum_{m=0}^M a_m y_{k-m} = \sum_{n=0}^N b_n x_{k-n}$$ +Usually $a_0$ is normalized to 1, and $y_k$ can be written explicitly. +$$y_k = \sum_{n=0}^N b_n x_{k-n} - \sum_{m=1}^M a_m y_{k-m}$$ +This iterative process can be computed for each sample + +For $N, M = 2$, this is a [biquad filter](https://en.wikipedia.org/wiki/Digital_biquad_filter), which is very fast, numerically stable (assuming the transfer function itself is mathematical stable), and a reasonably good sounding filter. + +$$H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}$$ + #### FIR filters -A finite impulse response (FIR) filter is a digital filter with $M = 0$ (a transfer function denominator of 1). For an input $x_k$ and output $y_k$, +A finite impulse response (FIR) filter is a specific case of an IIR filter with $M = 0$ (a transfer function denominator of 1). For an input and output signal, $$y_k = \sum_{n=0}^N b_n x_{k-n}$$ They are computationally straightforward and always stable since they have no poles. @@ -153,17 +166,45 @@ Note that the above formula is the convolution between vectors $y$ and $b$, and $$y \ast b = \mathcal{F}^{-1} \{ \mathcal{F}\{y\} \cdot \mathcal{F}\{b\} \}$$ where $\cdot$ is element-wise multiplication. -While the naive FIR formula above is $O(N^2)$ when processing blocks of $N$ samples, the FFT FIR method is $O(N \log N)$. -A disadvantage of the FFT FIR method is that the signal must be delayed by $N$ samples to produce any output. -You can combine the naive and FFT methods into a hybrid approach with the [overlap-add](https://en.wikipedia.org/wiki/Overlap%E2%80%93add_method) or [overlap-save](https://en.wikipedia.org/wiki/Overlap%E2%80%93save_method) methods. +While the naive FIR formula above requires $O(N)$ computations per sample, the FFT FIR method requires an average of $O(\log N)$ computations per sample when processing blocks of $N$ samples. +A disadvantage of the FFT FIR method is that the signal must be delayed by $N$ samples before any output can be determined. +You can combine the naive and FFT methods into a hybrid approach with the [overlap-add](https://en.wikipedia.org/wiki/Overlap%E2%80%93add_method) or [overlap-save](https://en.wikipedia.org/wiki/Overlap%E2%80%93save_method) methods, allowing you to process smaller blocks of size $M < N$ at a slightly worse $O((N/M) \log M)$ average computations per sample. -#### IIR filters -An infinite impulse response (IIR) filter is a general rational transfer function. By multiplying the denominator of the rational $H(z)$ definition above on both sides and applying it to an input and output signal, we obtain -$$\sum_{m=0}^M a_m y_{k-m} = \sum_{n=0}^N b_n x_{k-n}$$ -Usually $a_0$ is normalized to 1, and $y_k$ can be written explicitly. -$$y_k = \sum_{n=0}^N b_n x_{k-n} - \sum_{m=1}^M a_m y_{k-m}$$ +### Impulse responses -For $N, M = 2$, this is a [biquad filter](https://en.wikipedia.org/wiki/Digital_biquad_filter), a very fast, numerically stable (assuming the transfer function itself is mathematical stable), and reasonably good sounding filter. +Sometimes we need to simulate non-rational transfer functions. +Consider a general transfer function $H(f)$, written in terms of $f$ rather than $s$ for simplicity. +Our original definition of a transfer function for an input $x(t)$ and output $y(t)$ is +$$ Y(f) = H(f) X(f) $$ +and by applying the convolution theorem again, we obtain +$$ y(t) = (h \ast x)(t) = \int_{-\infty}^\infty h(\tau) x(t - \tau) d\tau $$ +where $h(t)$ is the *impulse response* of our filter. -$$H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}$$ +The signal $h(t)$ is the result of processing a [delta function](https://en.wikipedia.org/wiki/Dirac_delta_function) through our filter, since $\delta(t)$ is the "identity" function of convolution, i.e. $h(t) = (h \ast \delta)(t)$. +Clapping your hands or popping a balloon (both good approximations of $\delta$) in a large cathedral will generate a very sophisticated impulse response, which can be recorded and processed in a FIR filter algorithm to reproduce arbitrary sounds as if they were performed in the cathedral. + +Repeating this process in the digital realm gives us the discrete convolution. +$$ y_k = \sum_{n=-\infty}^\infty h_n x_k $$ +Note that $h_n$ is both non-causal (nonzero for negative $t$ or $n$) and infinitely long, which is addressed later. + + +#### Brickwall filter + +An example of a non-rational transfer function is the ideal lowpass filter that fully attenuates all frequencies higher than $f_c$ and passes all frequencies below $f_c$. +Including [negative frequencies](https://en.wikipedia.org/wiki/Negative_frequency), its transfer function is +$$ +H(f) = \begin{cases} + 1 & \text{if } -f_c \leq f \leq f_c \\ + 0 & \text{otherwise} +\end{cases} +$$ +The inverse Fourier transform of $H(f)$ is +$$ h(t) = 2 f_c \operatorname{sinc}(2 f_c t) $$ +where $\operatorname{sinc}(x) = \sin(\pi x) / (\pi x)$ is the [normalized sinc function](https://en.wikipedia.org/wiki/Sinc_function). +Although the impulse response is infinitely long, restricting it to a finite range $[-T, T]$ and shifting it forward by $T$ produces a finite causal impulse response that can be solved by a fast FIR algorithm to produce a close approximation of an ideal brickwall filter. + + +### Window functions + +*Coming soon* \ No newline at end of file