Simple digital filter design

06 Sep 2019 - tsp

What are digital filters?

Digital filters employ the idea of simulating traditional analog filter responses with digital calculations. The basic idea is to sample a signal with an analog digital converter (ADC) and perform the calculations numerically. This allows a more flexible setup (for example on can change the filter response at runtime) as well as a more cheap setup since less analog components are used. These filters are normally limited only by Nygists sampling theorem that says that one cannot process signals with a frequency above half the sampling frequency after discretization as well as available processing power. Note that you should normally employ at least an analogous lowpass filter in front of your ADCs to remove sampling artifacts from high frequency components.

Digital filters are normally specified as a filter function that can be applied to the input signal. This specification is done in frequency space as a function $H(\omega)$ which is simply multiplied to the input signal $x(\omega)$. This leads to an output signal

$y(\omega) = x(\omega) H(\omega)$

Since calculation in frequency space would require one to transform the input signal into frequency space, multiply it and then perform a backtransformation which would not allow a continuous evaluation of the response one can use the convolution (since it’s known that the timedomain convolution of two signals equals the frequency domain multiplication):

$y(t) = x(t) * H(t) = \int_{-\infty}^{\infty} x(\tau) H(t-\tau) d\tau$

Later on we will do the convolution in discretized space - this can be done by moving from the integral (an infinite sum) to a discretized sum:

$y(n) = x(n) * H(n) = \sum_{\tau} x(\tau) H(n - \tau)$

Recalling the (discrete) fourier transform

First we have to recall how fourier transformation worked. First we define our quantities:

A function in time domain $f(t)$ where $t \in \mathfrak{R}$ is the continuous time and $f(t)$ can be any arbitrary amplitude.
A function in frequency domain $f(\omega)$ with $\omega \in \mathfrak{R}$ being the continous normalized frequency.

Note that we have defined normalized frequency since the results can be used that way in a general way. For generalized normalization we have to know the sampling frequency $f_s$ at which our digital system works. Due to the Nyquist sampling theorem we know the highest frequency we can work with is $\frac{f_s}{2}$. Since the normal definition of the angular frequency would be $\omega = 2 \pi f$ we normalize using $\omega = 2 \pi \frac{f}{f_s}$. This yields a useable normalized frequency range from $-\pi$ to $\pi$.

The continous transformation between frequency and time domain is done via the fourier transformation. This is based on the assumption that each signal can be represented by an infinite series of sine and cosine functions:

$f(\omega) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} f(t) e^{-i 2 \pi \omega t} dt$

The factor $\frac{1}{\sqrt{2\pi}}$ compensated for our usage of angular frequency instead of frequency itself. To transform a signal back from frequency domain to time domain we use the inverse transformation:

$f(t) = \int_{-\infty}^{\infty} f(\omega) e^{-i 2 \pi \omega t} d\omega$

Discrete transformations

Since we use a digital system and don’t have a continuous signal available we could use discrete transformations. This will be especially interesting if we don’t transform basic signals and want to calculate these numerically or if we want to do a frequency anlysis or digital synthesis of a signal. To move from continuous to discrete systems we simply introduce the sampling number $n$ (which is - using the sample count per unit time $N$ - is related to continuous time via $t = \frac{1}{N} n = \delta t n$). The same way we introduce the wave vector $k$ as counterpart to the continuous frequency $\omega$ - with definition $\omega_0 = \frac{2\pi}{N}$ and $\omega = \omega_0 k$.

Now we can simply substitute the integrals (which are just infinite sums anyways) by discrete sums. This will of course lead to known discretization error:

$f(k) = \sum_{n=0}^{N-1} f(n) e^{-i k \omega_0 n}$

$f(n) = \frac{1}{N} \sum_{k=0}^{N-1} f(k) e^{i k \omega_0 n}$

Designing filters by ideal response functions

Low pass filter

The ideal response for a lowpass filter is to filter all frequencies whose absolute value lies above a given threashold - the cutoff frequency $\omega_c$. We can use the impulse response of our filter:

[ H(\omega) = \left\{ \begin{align} 1, \mid\omega\mid \leq \omega_c \\ 0, else \end{align} \right. ]

This is the response of the filter that we expect when putting in an impulse $\delta(t)$ into out filter. Since we know that the fourier transform of $\delta(t)$ just leads to a constant:

$\int_{-\infty}^{\infty} \delta(t) e^{-i 2\pi \omega t} dt = e^{0} = 1$

Since we want to perform the convolution in time domain we can simply perform the backtransformation:

$H(t) = \int_{-\infty}^{\infty} H(\omega) e^{i2\pi \omega t} d\omega = \int_{-\omega_c}^{\omega_c} e^{i 2\pi t \omega} d\omega$ $ = \frac{e^{i 2 \pi t \omega}}{2i \pi t} \mid_{-\omega_c}^{\omega_c} = \frac{1}{\pi t} \frac{e^{i 2\pi t \omega_c} - e^{-i 2\pi t \omega_c}}{2i}$

Using Eulers formula $e^{i\phi} = \cos(\phi) + i \sin(\phi)$ and it’s presentation $sin(\phi) = \frac{e^{i\phi} - e^{-i\phi}}{2i}$ we can conclude that

$H(t) = \frac{1}{\pi t} \sin(2 \pi t \omega_c)$

As we can see the function extends into both the positive and negative time domain. This means this filter is not causal since its output depends on future values. Such filters can be calculated exactly on previously calculated signals but not on samples filtered in real time.

Using $f_c = \frac{1}{2} f_s$ i.e. $\omega_c = 2 \pi \frac{f}{f_s} = \pi$ we could calculate the following response function:

Time domain representation of an ideal lowpass filter

Since - in applications that are not working with previously recorded data - we are using only old known values we have to make the filter causal. To do this we simply drop all values that depend on positive times. We also have to limit ourself to a fixed number of simples to we cut the coefficients after a given number of samples. The filter coefficients will simply be extracted from our transformed function at the given number of timesteps. This leads to a finite impulse response (FIR) filter. Note that such a construction will lead to errors in the filter function.

Bandpass filter

Since we will construct our high pass filter (due to some mathematical constraints) as a bandpass filter we will first look at the filter coefficients for a bandpass filter. Again we start with our filter specification in frequency space:

[ H(\omega) = \left\{ \begin{align} 0, \mid\omega\mid < \omega_c \\ 1, \omega_c \leq \mid\omega\mid \leq \omega_N \\ 0, \mid\omega\mid > \omega_N \end{align} \right. ]

Again we simply perform the inverse fourier transformation:

$H(t) = \int_{-\infty}^{\infty} H(\omega) e^{i 2\pi \omega t} d\omega = \int_{-\omega_N}^{-\omega_C} e^{i2\pi\omega t} d\omega + \int_{\omega_C}^{\omega_N} e^{i2\pi\omega t} d\omega$

Again this leads to an expression that we can transform into a sum of sines:

$H(t) = \frac{1}{\pi t} \left[ -\frac{e^{i2\pi \omega_c t} - e^{-i2\pi \omega_c t}}{2i} + \frac{e^{i2\pi \omega_N t} - e^{-i2\pi \omega_N t}}{2i} \right]$

$H(t) = \frac{1}{\pi t} \left[ -\sin(2\pi \omega_c t) + \sin(2 \pi \omega_N t) \right]$

Time domain representation of an ideal bandpass filter

As before we see that the function is not causal (i.e. it contains terms depending on future events). This can be solves as above introducing the same artifacts as above. A simple FIR filter can be constructed by dropping all elements depending on future events and sampling the function on a finite interval.

High pass filter

If we try to construct the highpass filter as above we get into a mathematical problem:

[ H(\omega) = \left\{ \begin{align} 1, \mid\omega\mid \geq \omega_c \\ 0, \omega < \omega_c \end{align} \right. ]

$H(t) = \frac{e^{-i2\pi t \omega_c t} - \overbrace{e^{-i 2 \pi t \infty}}^{\to 0}}{i 2 \pi t} + \frac{\overbrace{e^{i 2 \pi t \infty}}^{\to \infty} - e^{i 2 \pi \omega_c t}}{i 2 \pi t}$

As we see the second term diverges. Even when using the same trick as before we get into the problem of having to evaluate our sine at infinity. Since the limit of the sine doesn’t exist (it oscilates between $-1$ and $+1$) we cannot use this term to evaluate our filter expression. Of course we could argue that this term is not relevant because it’s representing the dependence of our filter on future events so the simple windowing method of simply cutting events from the future might resolve that problem. But there is another way that we can use to construct our filter in case we want to apply it to real-time sampled data: Make it a bandpass filter. Since we know that we are designing a time discrete digital filter we also know that there is an upper sampling frequency that we can depend on. That way we can set the upper frequency to the Nyquist frequency (i.e. half the sampling frequency) and use the above equation for bandpass filters to construct our digital highpass.