Normalised frequency

Note that the sampling frequency

is lost and needs to be stored additionally so that the analogue signal can be reconstructed:

What is the frequency range of our digital signal? What is the maximum frequency in digital terms? Fig 1 shows that the max frequency is \(0.5\) because we need two samples to represent a “wave”. The frequency range \(0\ldots 0.5\) is called the normalised frequency range. What is the relation between normalised frequency and sampling rate?

Let’s have a look at the analog signal \(X_{a}(t)\) with the frequency \(F\):

and we sample it only at:

Then the digital signal is:

Now, we can define the normalised frequency as:

which has its max value at \(0.5\) which represents one period of a sine wave within two samples ¹

Nyquist frequency

Recall that the max value for the normalised frequency \(f\) is \(0.5\) and that:

with \(n\) as an integer because we are sampling.

What happens above \(f{\gt}0.5\)? Imagine \(f = 1\)

which gives us DC or zero frequency. The same is true for \(f = 2, 3, 4, \ldots \). We see that above \(f=0.5\) we never get higher frequencies. Instead they will always stay between \(0\ldots 0.5\) for the simple reason that it is not possible to represent higher frequencies. This will be discussed later in greater detail.

The ratio \(F/F_s\) must be lower than \(0.5\) to avoid ambiguity or in other words the maximum frequency in a signal must be lower than \(\frac{1}{2} F_s\). This is the Nyquist frequency.

If there are higher frequencies in the signal then these frequencies are “folded down” into the frequency range of \(0\ldots \frac{1}{2} F_s\) and creating an alias of its original frequency in the so called “baseband” (\(f=0\ldots 0.5\)). As long as the alias is not overlapping with other signal components in the baseband this can be used to downmix a signal. This leads to the general definition of the sampling theorem which states that the bandwidth \(B\) of the input signal must be half of the sampling rate \(F_s\):

The frequency \(\frac{1}{2}F_s\) is called the Nyquist frequency.

What do we do if the signal contains frequencies above \(\frac{1}{2} F_s\)? There are two ways to tackle this problem: The classical way is to use a lowpass filter (see Fig. 2A) which filters out all frequencies above the Nyquist frequency. However this might be difficult in applications with high resolution A/D converters. Alternatively one can use a much higher sampling rate to avoid aliasing. This is the idea of the sigma delta converter which operates at sampling rates hundred times higher than the Nyquist frequency.

Sampling theorem

Is it possible to reconstruct an analogue signal from a digital signal which contains only frequencies below the Nyquist frequency?

where \(F_{\mbox{max}}\) is max frequency in the signal which we represent by sine waves:

The analogue signal \(x_{a}(t)\) can be completely reconstructed if:

with

The problem is that \(g(t)\) runs from negative time to positive time and as we see later is a-causal so that this cannot be implemented for real but approximations of \(g(t)\) are possible and are analogue lowpass filters which smooth out the step like outout of an digital to analogue converter.

Periodic signals

Periodic signals can be composed of sine waves :

where \(F_1\) is the principle or fundamental frequency and \(k \neq 1\) are the harmonics with strength \(c_k\). Usually \(c_1\) is set to \(1\). This is a Fourier series with \(c_{k}\) as coefficients. For example an ECG has a fundamental frequency of about 1Hz (60 beats per minute). However, the harmonics give the ECG its characteristic peak like shape.

How do we find the coefficients \(c_k\)?

For simple cases there are analytical solutions for \(c_k\), for example for square waves, triangle wave, etc.

What are the properties of \(c_{k}\)?

or

Proof: with the handy equation…

we get

How are the frequencies distributed? Let’s have a look at the frequency spectrum of a periodic signal: \(P_{k} = \mid c_{k} \mid ^{2}\)

• There are discrete frequency peaks

• Spacing of the peaks is \(\frac{1}{T_1} = F_1\)

• Only the positive frequencies are needed: \(c_{-k} = c_{k}^{*}\)

A-periodic signals

In case nothing is known about \(X(t)\) we need to integrate over all frequencies instead of just the discrete frequencies.

Consequently, the frequency spectrum \(X(F)\) is continuous.

Discrete time Fourier Transform (DTFT)

The signal \(x(n)\) is discrete whereas the resulting frequency spectrum is considered as continuous (arbitrary signals).

• Analysis or direct transform:

  $$\begin{equation}X(\omega )=\sum _{n=-\infty }^{\infty } x(n) e^{-j\omega n}\tag{37}\end{equation}$$

• Synthesis or inverse transform:

  $$\begin{align} \tag{38} x(n) & = & \frac{1}{2\pi } \int _{-\pi }^{\pi } X(\omega ) e^{j\omega n} d\omega \\ \tag{39}& = & \frac{1}{2\pi } \int _{-0.5}^{0.5} X(f) e^{j 2\pi f n} df \end{align}$$

  note the range here. \(f\) is the normalised frequency.

The effect of time domain sampling on the spectrum (Sampling theorem)

What effect has time domain sampling on the frequency spectrum ² ? Imagine we have an analog spectrum \(X_a(F)\) from a continuous signal \(x(t)\). We want to know how the spectrum \(X(F)\) of the discrete signal \(x(n)\) looks like.

The signal is represented by \(x(n)\) so that we can equate the Fourier transforms of the sampled spectrum \(X(F)\) and of the analogue spectrum \(X_a(F)\).

Obviously \(X(f)\) must be different to accommodate the different integration ranges. The trick is now to divide the integral on the right hand side of Eq. 41 into chunks to make it compatible to the range on the left hand hand side.

Remember that the normalised frequency is \(f=F/F_s\) which allows us to change the integration to analogue frequency on both sides:

and now we divide the right hand side into chunks of \(F_s\) which corresponds to the integration range on the left hand side.

This gives us now an equation for the sampled spectrum:

This equation can now be interpreted. In order to get the sampled spectrum \(X(F)\) we need to make copies of the analog spectrum \(X_a(F)\) and place these copies at multiples of the sampling rate \(F_s\) (see Fig. 5). This illustrates also the sampling theorem: if the bandwidth of the spectrum is wider than \(F_s/2\) then the copies of the analogue spectrum will overlap and reconstruction would be impossible. This is called aliasing. Note that it is not necessary bad that the spectrum of the analogue signal lies within the range of the so called “base band” \(-F/2 \ldots F/2\). It can also lie in another frequency range further up, for example \(-F/2 + 34 \ldots F/2 +34\) as long as the bandwidth does not exceed \(F_s/2\). If it is placed further up it will automatically show up in the baseband \(-F/2 \ldots F/2\) which is called “fold down”. This can be used for our purposes if we want to down mix a signal.

With the insight from these equations we can create a plot of how analogue frequencies map onto sampled frequencies. Fig 6 shows how the analogue frequencies \(F_\textrm {analogue}\) map on the normalised frequencies \(f_\textrm {sampled}\). As long as the analogue frequencies are below \(F_s/2\) the mapping is as usual as shown in Fig 6A. Between \(F_s/2\) and \(F_s\) we have an inverse mapping: an increase in analogue frequency causes a decrease in frequencies. Then, from \(F_s\) we have again an increase in frequencies starting from DC. So, in general if we keep a bandlimited signal within one of these slopes (for example from \(F_s \ldots F_s + 1/2 F_s\) as shown in Fig 6B) then we can reproduce the signal.

This leads us to the generalised Nyquist theorem: if a bandpass filtered signal has a bandwidth of \(B\) then the minimum sampling frequency is \(F_s = 2B\).

Discrete Fourier Transform (DFT)

So far we have only sampled in the time domain. However, on a digital computer the Fourier spectrum will always be a discrete spectrum.

The discrete Fourier Transform (DFT) is defined as:

where \(N\) is the number of samples in both the time and frequency domain.

The inverse discrete Fourier Transform (IDFT) is defined as:

What is the effect in the timedomain of this discretisation ³ ? We start with the continuous Fourier transform and discretise it into N samples in the frequency domain:

Let’s subdivide the sum into chunks of length \(N\):

We note the following:

• Ambiguity in the time domain

• The signal is repeated every N samples

Practically this repetition won’t show up because the number of samples is limited to \(N\) in the inverse transform. However, for operations which shift signals in the frequency domain it is important to remember that we shift a periodic time series. If we shift it out at the end of the array we will get it back at the start of the array.

Fig. 8 shows an example of a DFT. It is important to know that the spectrum is mirrored around \(N/2\). DC is represented by \(X(0)\) and the Nyquist frequency \(F_s/2\) is represented by \(X(N/2)\). The mirroring occurs because the input signal \(x(n)\) is real. This is important if one wants to modify the spectrum \(X(F)\) by hand, for example to eliminate 50Hz noise. One needs to zero two elements of \(X(F)\) to zero. This is illustrated in this python code:

import scipy as sp
yf=sp.fft(y)
# the sampling rate is 1kHz. We've got 2000 samples.
# midpoint at the ifft is 1000 which corresponds to 500Hz
# So, 100 corresponds to 50Hz
yf[99:101+1]=0;
# and the mirror
yf[1899:1901+1]=0;
#
yi=sp.ifft(yf);

This filters out the 50Hz hum from the signal \(y\) with sampling rate 1000Hz. The signal yi should be real valued again or contain only very small complex numbers due to numerical errors.

Properties of the DFT

• Periodicity:

  $$\begin{align} \tag{54} x(n + N) & = & x(n) \\ \tag{55}X(k + N) & = & X(k) \end{align}$$

• Symmetry: if \(x(n)\) real:

  $$\begin{equation}x(n)\mbox{ is real} \Leftrightarrow X^*(k) = X(-k) = X(N-k)\tag{56}\end{equation}$$

  This is important when manipulating \(X(k)\) by hand.

• Time Reversal:

  $$\begin{align} \tag{57} x(n) & \leftrightarrow & X(k) \\ \tag{58}x(-n)& \leftrightarrow & X(N - k) \end{align}$$

• Circular convolution:

  $$\begin{equation}X_{1}(k) X_{2}(k) \leftrightarrow x_1(n) * x_2(n)\tag{59}\end{equation}$$

  with

  $$\begin{equation}x_{3}(m) = \sum _{n = 0}^{N - 1} x_{1}(n) x_{2}(m - n)\tag{60}\end{equation}$$

More useful equations are at Proakis and Manolakis ( 1996 , pp.415 ) .

Problems with finite length DFTs

where \(x\) is the input signal and \(w(n)\) represents a function which is \(1\) from \(n=0\) to \(n=L-1\) so that we have \(L\) samples from the signal \(x(n)\).

To illustrate the effect of this finite sequence we introduce a sine wave \(x(n) = \cos \omega _0 n\) which has just two peaks in a proper Fourier spectrum at \(-\omega \) and \(+\omega \).

The spectrum of the rectangular window with the width \(L\) is:

The resulting spectrum is then:

Because the spectrum of \(X(\omega )\) consists of just two delta functions the spectrum \(X_{3}(\omega )\) contains the window spectrum \(W(\omega )\) twice at \(-\omega \) and \(+\omega \) (see Fig. 9). This is called leakage. Solutions to solve the leakage problem? Use Windows with a narrower spectrum and with less ripples (see FIR filters).

Fast Fourier Transform

We can rewrite the DFT (Eq. 48) in a slightly more compact form:

with the constant:

The problem with the DFT is that it needs \(N^2\) multiplications. How can we reduce the number of multiplications? Idea: Let’s divide the DFT in an odd and an even sequence:

which gives us with the trick \(W_N^{2mk} = W_{N/2}^{mk}\) because of the definition Eq. 65.

\(F_e\) and \(F_o\) have both half the length of the original sequence and need only \((N/2)^2\) multiplication, so in total \(2 \cdot (N/2)^2 = \frac{N^2}{2}\). Basically by dividing the sequence in even and odd parts we can reduce the number of multiplications by 2. Obviously, the next step is then to subdivide the two sequences \(F_{e}(k)\) and \(F_{o}(k)\) even further into something like \(F_{ee}(k), F_{eo}(k), F_{oe}(k)\) and \(F_{oo}(k)\).

In general the recipe for the calculation of the FFT is:

\(W_{L}^{k}\) is the phase factor in front of the odd sequence. This is continued until we have only two point (\(N=2\)) DFTs (see Eq. 64):

A two point DFT operates only on two samples which can represent only two frequencies: DC and the Nyquist frequency which makes the calculation trivial. Eq. 72 is an averager and Eq. 73 is basically a differentiator which gives the max output for the sequence \(1,-1,1,-1,\ldots \). Fig. 10 illustrates how to divide the initial sequence to arrive at 2 point DFTs. In other words the calculation of the full DFT is done by first calculating \(N/2\) 2 point DFTs and recombining the results with the help of Eq. 71. This is sometimes called the “Butterfly” algorithm because the data flow diagram can be drawn as a butterfly. The number of complex multiplications reduces in this approach to \(N \log _2 N\) which is actually the worst case scenario because many \(W_{L}^{k}\) usually turn out to be \(1,-1,j,-j\) which are just sign inversions or swaps of real and imaginary parts. A clever implementation of this algorithm will be even faster.

In summary the idea behind the FFT algorithms is to divide the sequence into subsequences. Here we have presented the most popular radix 2 algorithm. The radix 4 is even more efficient and there are also algorithms for divisions into prime numbers and other rather exotic divisions. However, the main idea is always the same: subsample the data in a clever way so that the final DFT becomes trivial.

How to characterise filters?

1. Impulse Response

   $$\begin{align} \tag{88} x(t) & = & \delta (t) \qquad \leftarrow \mbox{delta pulse} \\ \tag{89}h(t) & = & y(t) \qquad \leftarrow \mbox{impulse response} \end{align}$$

   The filter is fully characterised by its impulse response \(h(t)\)

2. Transfer function The Laplace transform of the impulse response is called Transfer Function. With the argument \(j\omega \) we get the frequency response of the filter. What does the frequency response tell us about the filter? The absolute value of the

   $$\begin{equation}|H(i \omega )|\tag{90}\end{equation}$$

   gives us the amplitude or magnitude for every frequency (compare the Fourier transform). The angle of the term \(H(j\omega )\) gives us the phase shift:

   $$\begin{equation}\phi = \arg \left(H(i \omega ) \right)\tag{91}\end{equation}$$

   of the filter. In this context the group delay can be defined as:

   $$\begin{equation}\tau _{\omega } = - \frac{d \phi (\omega )}{d\omega }\tag{92}\end{equation}$$

   which is delay for a certain frequency \(\omega \). In many applications this should be kept constant for all frequencies.

Generally FIR filters do not perform a convolution operation

Remember that we have derived the FIR filters from the analogue domain by filtering a signal (Eq.105). The timedomain representation of this is a discrete and causal convolution operation (using the sampled Laplace transform). However, the problem is that this runs to infinite time (see Eq 110) and would require an infinite number of delay steps and thus infinite time to complete. So the FIR filter at best performs an approximation of the convolution operation but there are serious problems with the finite number of taps (\(N\)) which requires a technique called windowing which moves it even further away from a convolution operation. On the other hand if the impulse response is strictly limited in time one can use this property for example for matched filters. However, the bottomline is that sweeping statements such as that “FIR filters perform a convolution” are generally wrong if nothing is known of the impulse response.

FIR filter implementations

No matter the design process the implementation is always the same: we need a delay line for the incoming signal and then weight the delayed outputs by the different coefficients and for that reason are called (water-) “taps” (see Eqs. 107 & 108). We are now presenting different implementations.

• C++: This is a simple example of a filter which stores the values in a simple linear buffer bufferFIR which stores the delayed values. The coefficients are stored in coeffFIR.

    float filter(float value) {
      // shift
      for (int i=taps-1;i>0;i--)  {
        bufferFIR[i]=bufferFIR[i-l];
      }
      //store new value
      bufferFIR[0]=value;
      //calculate result
      for (int i=0;i<taps;i++)  {
        output +=bufferFIR[i]*coeffFIR[i];
      }
    return output;
    }
  

• Python: Here, the FIR filter is implemented as a class which receives the FIR filter coefficients in the constructor and then filters a signal sample by sample in the function filter:

  class FIR_filter:
      def __init__(self,_coefficients):
          self.ntaps = len(_coefficients)
          self.coefficients = _coefficients
          self.buffer = np.zeros(self.ntaps)
  
      def filter(self,v):
          self.buffer = np.roll(self.buffer,1)
          self.buffer[0] = v
          return np.inner(self.buffer,self.coefficients)
  

  which again processes the signal sample by sample. It uses the numpy “roll” command to shift the samples and then the inner product to calculate the weighted sum between the buffer and the coefficients.

  If one wants to filter a whole array one can use Python’s lfilter command:

  import scipy.signal as signal
  y = signal.lfilter(h,1,x)
  

  This filters the signal x with the impulse response h. Note that this operation is on an array and thus a-causal.

More sophisticated code can be found in Teukolsky et al. ( 2007 ) . This book is strongly recommended for any C programmer who needs efficient solutions.

Fixed point FIR filters

These filters receive integer numbers as input, perform integer multiplications/additions and their outputs are integer as well. Thus, these filters do not require a floating point unit on a processor.

Fig. 18 shows a fixed point FIR filter. The input \(x(n)\) is an integer variable with I bit integers, the accumulator is an integer variable with A bits and the output as well (usually the same as the input in terms of bit width).

In contrast to a floating point FIR filter we need to scale up the coefficients so that they use full the integer range to avoid quantisation errors. For example if the coefficients of \(h(n)\) range from \(-0.75\) and \(+0.75\) and we have signed 16 bit integers then the scaling factor is \(2^W, W=15\).

However, the accumulator A which collects the data needs to have more bits because it receives scaled input values at I bits precision and these multiplied by factor \(2^W\). If we have \(M\) taps then the additional bits we need is \(\log _2 M\). The total number of bits we need in the accumulator in the worst case are:

However, this is the worst case scenario because if the gain of the FIR filter is below one then the summations by the \(M\) taps will only create temporary overflows because integer numbers are cyclic in their representation. In case the gain of the FIR filter is below one this can be relaxed:

The actual bitwidth of the accumulator is usually the next integer size available and also makes sure that in case the gain goes slightly over one in an unexpected case that the filter still works. For example if I has \(16\) bits the accumulator has probably \(32\) bits.

Constant group delay or linear phase filter

So far the FIR filter has no constant group delay which is defined by:

This means that different frequencies arrive at the output of the filter earlier or later. This is not desirable. The group delay \(\tau _\omega \) should be constant for all frequencies \(\omega \) so that all frequencies arrive at the same time at the output \(y\) of the filter.

A constant group delay can be achieved by restricting ourselves to the transfer function:

where \(B(\omega )\) is real and the phase is only defined by the exponential. The phase of this transfer function is then trivially \(\omega \tau +\phi \). The group delay is the derivative of this term which is constant.

Eq. 114 now imposes restrictions on the impulse response \(h(t)\) of the filter. To show this we use the inverse Fourier transform of Eq. 114 to get the impulse response:

After some transformations we get:

where \(b(n)\) represents the Fourier coefficients of \(B(\omega )\). Since \(B\) is real the coefficients \(b(n)\) have the property \(b(n)=b^*(-n)\). Thus we get a second impulse response:

Now we can eliminate \(b\) by equating Eq. 116 and Eq. 117 which yields:

With the shift \(\tau \) we have the chance to make the filter “more” causal. We can shift the impulse response in positive time to get \(h(n)\) zero for \(n{\lt}0\). In a practical application we shift the impulse response by half the number of delays. If we have a filter with \(M\) taps we have to delay by \(\tau =M/2\).

The factor \(e^{2i\phi }\) restricts the values of \(\phi \) because the impulse response must be real. This gives us the final FIR design rule:

where \(k=0\) or \(1\). This means that the filter is either symmetric or antisymmetric and the impulse response has to be delayed by \(\tau =M/2\).

Window functions

So far we still have a infinite number of coefficients for for the FIR filter because there’s no guarantee that the impulse response becomes zero after \(M/2\) delays. Thus, we have to find a way to truncate the response without distorting the filter response.

The standard technique is to multiply the coefficients with a window function which becomes and stays zero at a certain coefficient \(n{\gt}N\) so that Eq. 106 need not to run to infinity:

1. Rectangular window: truncating the impulse response. Problem: we get ripples in the frequency-response. The stop-band damping is poor

2. Triangular or Bartlett window: greatly improved stop-band attenuation.

3. Hanning and Hamming window: standard windows in many applications.

   $$\begin{equation}w(n) = \alpha - (1-\alpha ) \cos \left(\frac{2\pi n}{M}\right)\tag{121}\end{equation}$$

   • Hamming: \(\alpha = 0.54\)

   • Hanning: \(\alpha = 0.5\)

4. Blackman window:

   $$\begin{equation}w(n) = 0.42 + 0.5 \cos \left(\frac{2 \pi n}{M}\right) + 0.08 \cos \left( \frac{4 \pi n}{M} \right)\tag{122}\end{equation}$$

5. Kaiser window: control over stop- and passband. No closed form equation available.

To illustrate how window functions influence the frequency response we have taken an impulse response of a lowpass filter (\(f_c=0.1\)) and applied different window functions to it (Fig. 19).

Note that the higher the damping the wider the transition from pass- to stopband. This can be seen when comparing the Blackman window with the Hamming window (Fig. 19). For the lowpass filter this seems to be quite similar. However, for a bandstop filter the wider transition width might lead actually to very poor stopband damping. In such a case a Hamming window might be a better choice.

Python code: impulse response from the inverse DFT - The frequency sampling method

Imagine that we want to remove \(50Hz\) from a signal with sampling rate of \(1kHz\). We define an FIR filter with 100 taps. The midpoint \(N/2=50\) corresponds to \(500Hz\). The \(50Hz\) correspond to index \(5\).

f_resp=np.ones(100)
# note we need to add "+1" to the end because the end of
# the range is not included.
f_resp[4:6+1]=0
f_resp[94:96+1]=0
hc=np.fft.ifft(f_resp)
h=np.real(hc)
# this is from index 0 to index 49 on the left
# and on the right hand side from index 50 to index 99
h_shift[0:50]=h[50:100]
h_shift[50:100]=h[0:50]
h_wind=h_shift*hamming(100)

To get a nice symmetric impulse response we need to shift the inverse around \(50\) samples.

FIR filter design from ideal frequency response – The analytical way

For many cases the impulse response can be calculated analytically. The idea is always the same: define a function with the ideal frequency response

and perform an inverse Fourier transform to get the impulse response \(h(n)\).

We demonstrate this with a lowpass filter:

Use the inverse Fourier transform to get the impulse response:

With these handy equations:

we get for the filter:

This response is a-causal! However we know that we can shift the response by any number of samples to make it causal! And we have to window the response to get rid of any remaining negative contribution and to improve the frequency response.

Highpass, bandstop and bandpass can be calculated in exactly the same way and lead to the following ideal filter characteristics:

• Lowpass with cutoff frequency \(\omega _c=2\pi f_c\):

  $$\begin{equation}h(n) = \left\{ \begin{array}{ll} \frac{\omega _c}{\pi } & \mbox{for $n=0$} \\ \frac{1}{\pi n} \sin (\omega _c n) & \mbox{for $n\neq 0$} \end{array} \right. \label{idealLP}\tag{132}\end{equation}$$

• Highpass with the cutoff frequency \(\omega _c=2\pi f_c\):

  $$\begin{equation}h(n)= \left\{ \begin{array}{ll} 1-\frac{\omega _c}{\pi } & \mbox{for $n=0$} \\ -\frac{1}{\pi n} \sin (\omega _c n) & \mbox{for $n\neq 0$} \end{array} \right.\tag{133}\end{equation}$$

• Bandpass with the passband frequencies \(\omega _{1,2}=2\pi f_{1,2}\):

  $$\begin{equation}h(n)= \left\{ \begin{array}{ll} \frac{\omega _2-\omega _1}{\pi } & \mbox{for $n=0$} \\ \frac{1}{\pi n}(\sin (\omega _2 n)-\sin (\omega _1 n)) & \mbox{for $n\neq 0$} \end{array} \right.\tag{134}\end{equation}$$

• Bandstop with the notch frequencies \(\omega _{1,2}=2\pi f_{1,2}\):

  $$\begin{equation}h(n)= \left\{ \begin{array}{ll} 1-\frac{\omega _2-\omega _1}{\pi } & \mbox{for $n=0$} \\ \frac{1}{\pi n}(\sin (\omega _1 n)-\sin (\omega _2 n)) & \mbox{for $n\neq 0$} \end{array} \right.\tag{135}\end{equation}$$

See Diniz ( 2002 , p.195 ) for more impulse responses.

Here is an example code for a bandstop filter which fills the array h with the analytically calculated impulse response:

f1 = 45.0/fs
f2 = 55.0/fs
n = np.arange(-200,200)
h = (1/(n*np.pi))*(np.sin(f1*2*np.pi*n)-np.sin(f2*2*np.pi*n))
h[200] = 1-(f2*2*np.pi-f1*2*np.pi)/np.pi;
h = h * np.hamming(400)

After the function has been calculated it is windowed so that the cut off is smooth. It’s not very elegant as it’s causes a division by zero first and then the coeffecient at 200 is fixed. Do it better!

Design steps for FIR filters

Here are the design steps for an \(M\) tap FIR filter.

1. Get yourself an impulse response for your filter:

   1. Create a frequency response “by hand” just by filling an array with the desired frequency response. Then, perform an inverse Fourier transform (see section None).

   2. Define a frequency response analytically. Do an inverse Fourier transform (see section None).

   3. Use an analogue circuit, get its impulse response and use these as filter coefficients.

   4. Dream up directly an impulse response (for example, averagers, differentiators, etc)

2. Mirror the impulse response (if not already symmetrical)

3. Window the impulse response from an infinite number of samples to \(M\) samples.

4. Move the impulse response to positive time so that it becomes causal (move it \(M/2\) steps to the right).

FIR filter design with Python’s high level functions

The “firwin” command generates the impulse response of a filter with m taps and also applies a default window function. For example:

from scipy.signal import firwin
h = firwin(m,2*f)

generates a lowpass FIR filter with the normalised frequency f. Note the factor two because in scipy the normalisation is not the sampling rate but the nyquist frequency. To be compatible with the math here the easiest approach is to multiply the frequency by 2. With this command one can also design high pass, bandstop and bandpass filters. Type help(firwin) which provides examples for all filter types and which parameters need to be set.

Signal Detection: Matched filter

How can I detect a certain event in a signal? With a correlator. Definition of a correlator?

How to build a correlator with a filter? Definition of filtering?

NB. h - integration runs backwards. However we used an int forward! \(h(t) : = r(T - \tau )\), only valid for \(0\ldots T\).

for \(t: = T\) we get:

In order to design a detector we just create an impulse response \(h\) by reversing the template \(r\) in time and constructing an FIR filter with it.

How to improve the matching process? Square the output of the filter!

Adaptive FIR Least Mean Squares (LMS) filters

So far the coefficients of the FIR filter have been constant but now we are going to allow them to change while the FIR filter is operating so that it can learn its own coefficients. The idea is to use an error signal \(e(n)\) to tune the coefficients \(h(n)\) of the FIR filter (see Fig. 20A) by comparing its output \(y(n)\) with a desired output \(d(n)\):

and tuning the FIR filter in a negative feedback loop till the average of \(e(n)\) is zero.

We need to derive how the error \(e(n)\) can change the FIR filter coefficients so that the error is actually minimised. This can be expressed as a so called “gradient descent” which minimises the squared error \(\frac{1}{2} e(n)^2\) because both a positive and a negative error are equally bad:

where \(\Delta h_m(n)\) is the change of the coefficient \(h_m(n)\) at every time step: \(h_m(n+1) = h_m(n) + \Delta h_m(n)\). \(\mu \ll 1\) defines how quickly the coefficients \(h_m\) change at every time step and is called the “learning rate”. Note here the change in notation of the FIR filter coefficients \(h_m\) which are changing much slower than the sampled signals \(e(n), d(n)\) and \(y(n)\) and can still be seen as constant for the realtime filtering operation. Thus, we have \(n\) for the sample by sample processing as before and the index \(m\) of \(h_m\) for the very slowly changing FIR filter coefficients. To gain some intuition why Eq. 143 mimimises the error \(e(n)\) we look at what happens if we increase the FIR coefficient \(h_m\) a little bit, for example, by the small amount of \(|\epsilon | \ll 1\): \(h_m := h_m + \epsilon \). Then we can observe two cases:

1. The squared error \(e(n)^2\) increases so we need to decrease \(h_m\) as it makes it worse.

2. The squared error \(e(n)^2\) decreases so we need to increase \(h_m\) as it makes it better.

This also works in the same way if we decrease \(h_m\) by a small amount. Thus from both directions this always minimises the error. This is basically a carrot & stick approach where the coefficients \(h_m\) are being “rewarded” if they minimise the error and “punished” if they make the error larger. This approach is called “Least Mean Squares” (LMS) as it minimises the squared error on average. It’s also known from neural networks where the \(h_m\) are called the “weights” of a neuron (see delta rule).

Eq 143 now needs to be turned into a form which can directly run in software by solving its partial derivative by inserting Eqs. 142 and 110.

Note that \(x(n-m)\) emerges because of the chain rule when partially differentiating the output \(y(n)\) which depends on the sum of the delay lines of the FIR filter (see Eq 110). Eq. 147 is now our “learning” rule which can simply be applied to the FIR filter as showm in Fig 20B.

Now we have the recipe of an adaptive filter where the FIR filter minimises its own error \(e(n)\) in a negative feedback loop and generating an output \(y(n)\) which is as closely as possible to the desired signal \(d(n)\). The signal \(y(n)\) is often called “the remover” as it cancels out the signal \(d(n)\). However, not everything in \(d(n)\) is cancelled out but only what is correlated with the input \(x(n)\). This means that the error signal \(e(n)\) will not be zero but will contain any signal components which are not correlated with \(x(n)\). This means that the error signal actually is also the cleaned up output signal of the adaptive filter.

Imagine you want to remove \(50\) Hz mains from a signal

then one would provide \(50\) Hz mains via

so that then the filter learns to remove the powerline interference. One might think that the filter will remove everything from \(d(n)\) because it will try to minimise the error \(e(n)\) but only the 50 Hz in the contaminated signal will be removed because only signal components which are correlated between \(x(n)\) and \(d(n)\) will be removed. The cleaned up signal will be identical to the error signal:

Introduction

As an instructional example we take a very simple analogue filter and sample it. This can then be generalised to more complex situations.

We define a first order filter which can be implemented, for example, as a simple RC network:

where its Laplace transform is a simple fraction:

Now we sample Eq. 151:

and perform a Laplace transform on it:

which turns into a z-transform:

Consequently the analogue transfer function \(H(s)\) transfers into \(H(z)\) with the following recipe:

Thus, if you have the poles of an analogue system and you want to have the poles in the z-plane you can transfer them with:

This also gives us a stability criterion. In the Laplace domain the poles have to be in the left half plane (real value negative). This means that in the sampled domain the poles have to lie within the unit circle.

The same rule can be applied to zeros

together with Eq. 159 this is called “The matched z-transform method”. For example \(H(s)=s\) turns into \(H(z)=1-z^{-1}e^{0T}=1-z^{-1}\) which is basically a DC filter.

Determining the data-flow diagram of an IIR filter

To get a practical implementation of Eq. 157 we have to see it with its input- and output-signals:

We have to recall that \(z^{-1}\) is the delay by \(T\). With that information we directly have a difference equation in the temporal domain:

This means that the output signal \(y(nT)\) is calculated by adding the weighted and delayed output signal \(y([n-1]T)\) to the input signal \(x(nT)\). How this actually works is shown in Fig. 21. The attention is drawn to the sign inversion of the weighting factor \(e^{-bT}\) in contrast to the transfer function Eq. 158 where it is \(-e^{-bT}\). In general the recursive coefficients change sign when they are taken from the transfer function.

General form of IIR filters

A filter with forward and feedback components can be represented as:

where \(B_k\) are the FIR coefficients and \(A_l\) the recursive coefficients. Note the signs of the recursive coefficients are inverted in the actual implementation of the filter. This can be seen when the function \(H(z)\) is actually multiplied with an input signal to obtain the output signal (see Eq. 164 and Fig. 21). The “1” in the denominator represents actually the output of the filter. If this factor is not one then the output will be scaled by that factor. However, usually this is kept one.

In Python filtering is performed with the command:

import scipy.signal as signal
Y = signal.lfilter(B,A,X)

where B are the FIR coefficients, A the IIR coefficients and X is the input. For a pure FIR filter we just have:

Y = signal.lfilter(B,1,X)

The “1” represents the output.

IIR filter topologies

In real applications one would create a chain of 2nd order IIR filters. This has numerous advantages:

1. The optimisation can be done within these simple structures, for example omitting array operations completely and coding the 2nd order filters in assembly while keeping the higher level operations in C or C++.

2. Control of stability issues: high order recursive systems might become unstable and it is very hard to predict when this might happen. On the other hand a chain of 2nd order filters can be made stable with ease. One can focus on stability in 2nd order systems which is well understood and manageable.

3. Complex conjugate pole pairs are generated naturally in analogue filter design and directly translate to 2nd order IIR structures. Analogue design is usually described as 2nd order structures (=LCR) so that any transform from analogue to digital just needs to be of 2nd order!

Fig. 22 shows the most popular filter topologies: Direct Form I and II. Because of the linear operation of the filter one is allowed to do the FIR and IIR operations in different orders. In Direct From I we have one accumulator and two delay lines whereas in the Direct Form II we have two accumulators and one delay line. Only the Direct Form I is suitable for integer operations.

A Python class of a direct form II filter can be implemented with a few lines:

class IIR_filter:
    def __init__(self,_num,_den):
        self.numerator = _num
        self.denominator = _den
        self.buffer1 = 0
        self.buffer2 = 0

    def filter(self,v):
        input=0.0
        output=0.0
        input=v
        output=(self.numerator[1]*self.buffer1)
        input=input-(self.denominator[1]*self.buffer1)
        output=output+(self.numerator[2]*self.buffer2)
        input=input-(self.denominator[2]*self.buffer2)
        output=output+input*self.numerator[0]
        self.buffer2=self.buffer1
        self.buffer1=input
        return output

Here, the two delay steps are represented by two variables buffer1 and buffer2.

In order to achive higher order filters one can then just chain these 2nd order filters. In Python this can be achieved by storing these in an array of instances of this class.

Fixed point IIR filters

Fig. 23 shows the implementation of a fixed point IIR filter. It’s a Direct Form I filter with one accumulator so that temporary overflows can be compensated. The original floating point IIR coefficients are scaled up by factor \(2^w\) so that they max out the range of the integer coefficients. After the addition operation they are shifted back by \(w\) bits to their original values.

For example if the largest IIR coefficient is 1.9 and we use 16 bit signed numbers (max number is 32767) then one could multiply all coefficients with \(2^{14}\).

Then one needs to assess the maximum value in the accumulator which is much more difficult than for FIR filters. For example, resonators can generate high output values with even small input values so that it’s advisable to have a large overhead in the accumulator. For example if the input signal is \(16\) bit and the scaling factor is \(14\) bits then the signal will certainly occupy \(16+14=30\) bits. With an accumulator of 32 bits that gives only a headroom of 2 bits so the output can only be 4 times larger than the input. A 64 bit accumulator is a safe bet.

Filter design based on analogue filters

Most IIR filters are designed from analogue filters by transforming from the continuous domain (\(h(t),H(s)\)) to the sampled domain (\(h(n),H(z)\)).

A list of popular analogue transfer functions being used for IIR filter design:

• Butterworth: All poles lie on the left half plane equally distributed on a half circle with radius \(\Omega = 2\pi F\) which is shown in Fig. 24. The Butterworth filter is by far the most popular filter. Here are its properties:

  • monotonic frequency response

  • only poles, no zeros

  • the Poles have analytical solution and very easy to calculate

  • no constant group delay but usually acceptable deviation from a strict constant group delay for many applications.

• Chebyshev Filters:

  $$\begin{equation}|H(\Omega )|^{2} = \frac{1}{1 - \varepsilon ^{2} T_{N} (\Omega /\Omega _{p})}\tag{166}\end{equation}$$

  where T = Chebyslev polynomials. These filters have either ripples in the stop or passband depending on the choice of polynomials. As with Butterworth the polynomials have analytical solutions for the poles and zeros of the filter so that their design again is straightforward.

• Bessel Filter:

  • Constant Group Delay

  • Shallow transition from stop- to passband

  • No analytical solution for the poles/zeros.

How to transform the analogue lowpass filters into highpass, bandstop or bandpass filters? All the analogue transfer functions above are lowpass filters. However this is not a limitation because there are handy transforms which take an analogue lowpass transfer function or poles/zeros and then transforms them into highpass, bandpass and stopband filters. These rules can be found in any good analogue filter design book and industry handouts, for example from Analog Devices.

Bilinear Transform: transforming the analoge transfer function into a digital one: As a next step these analogue transfer functions \(H(s)\) need to be transformed to digital transfer functions \(H(z)\). This could be done by the matched z-transform. However, the problem with these methods is that they map frequencies 1:1 between the digital and analogue domain. Remember: in the sampled domain there is no infinite frequency but \(F_s/2\) which is the Nyquist frequency. This means that we never get more damping than at Nyquist \(F_s/2\). This is especially a problem for lowpass filters where damping increases the higher the frequency. The solution is to map all analogue frequencies from \(0 \ldots \infty \) to the sampled frequencies \(0\ldots 0.5\) in a non-linear way:

IIR filter design steps:

1. Choose the cut-off frequency of your digital filter \(\omega _{c}\).

2. Calculate the analogue cutoff frequency \(\omega _{c} \rightarrow \Omega _{c}\) with Eq. 171

3. Choose your favourite analogue lowpass filter \(H(s)\), for example Butterworth.

4. Replace all \(s\) in the analogue transfer function \(H(s)\) by \(s = \frac{2}{T} \frac{z - 1}{z + 1}\) to obtain the digital filter \(H(z)\)

5. Change the transfer function \(H(z)\) so that it only contains negative powers of \(z\) (\(z^{-1}, z^{-2}, \ldots \)) which can be interpreted as delay lines.

6. Build your IIR filter!

Time or frequency domain? The transforms from analogue to digital alter both the temporal response and the frequency response. However the different transforms have different impact. The bilinear transform guarantees that the digital filter uses the whole frequency range of the analogue filter from zero to infinity which is the best solution for frequency domain problems. However, if one needs to reproduce the temporal response an analogue filter as faithfully as possible then the matched z transform is best because it's based (as the name suggests) on the impulse invariance method which aims to preserve the temporal behaviour of the filter (exact for pole-only transfer functions).

Adaptive IIR filter: The Kalman filter

Often signals are contaminated by high frequency noise where the spectrum of the noise is changing or not known. The cutoff of a fixed low pass filter is not known or might be changing so that we need to adapt the cut-off continuously. One example is a Kalman filter which maximises the predictability of the filtered signal. It’s an adaptive lowpass filter which increases the predictability of a signal because it smoothes it.

The parameter \(a\) determines the cutoff frequency. Frequency Response :

Let’s re-interpret our low-pass filter in the time domain:

We would like to have the best estimate for \(y(n)\)

We need to minimise \(p\) which gives us equations for a and b which implements a Kalman filter.

The role of poles and zeros

Transfer functions contain poles and zeros. To gain a deeper understanding of the transfer functions we need to understand how poles and zeros shape the frequency response of \(H(z)\). The position of the poles also determines the stability of the filter which is important for real world applications.

We are going to explore the roles of poles and zeros first with an instructional example which leads to a 2nd order bandstop filter.

Zeros

Poles

Stability A transfer function \(H(z)\) is only stable if the poles lie inside the unit circle. This is equivalent to the analog case where the poles of \(H(s)\) have to lie on the left hand side of the complex plane (see Eq. 159). Looking at Eq. 183 it becomes clear that \(r\) determines the radius of the two complex conjugate poles. If \(r{\gt}1\) then the filter becomes unstable. The same applies to poles on the real axis. Their the real values have to stay within the range \(-1 \ldots +1\). In summary: poles generate resonances and amplify frequencies. The amplification is strongest the closer the poles move towards the unit circle. The poles need to stay within the unit circle to guarantee stability. Note that in real implementations the coefficients of the filters are limited in precision. This means that a filter might work perfectly in python but will fail on a DSP with its limited precision. You can simlate this by forcing a certain datatype in python or you write it properly in C.

Poles and zeroes combined: design of a pure digital IIR notch filter

Identifying filters from their poles and zeroes Proakis and Manolakis ( 1996 , pp.333 ) has an excellent section about this topic and we refer the reader to have a look. As a rule of thumb a digital lowpass filter has poles where their real parts are positive and a highpass filter has poles with negative real part of the complex poles. In both cases they reside within the unit circle to guarantee stability.