UCL Department of Phonetics and Linguistics
Introduction to Computer Programming with MATLAB
Lecture 8: Discrete Linear Systems
Objectives
By the end of the session you should:
q understand how an impulse response can be used as an alternative to the frequency response for the characterisation of a discrete linear system
q know how to implement a digital system from knowledge of its impulse response using convolution
q know how to build and apply simple non-recursive filters
q know how to build and apply simple recursive filters
q know how to plot the frequency response of a filter
Outline
1. Characterisation of Linear Systems
A discrete linear system is a digital implementation of a linear time-invariant system. Its input is a vector representing the sampled input signal (we will use x[]), its output is a vector of the same size as the input representing the sampled output signal (y[]).
The system itself has to obey the usual conventions of linearity, particularly that it obeys the principle of superposition, so that f(a+b) = f(a) + f(b), where f() is the processing performed by the system and a and b are arbitrary signals. Linear time-invariant systems can be characterised in the frequency domain by their frequency response. The frequency response of a system tells us how a sinusoidal input is changed in magnitude and phase as a function of its frequency. Although it is possible to implement systems in the frequency domain, it is more typical to implement them in the time domain: sample by sample. To do this we use the time domain equivalent of the frequency response, known as the impulse response. A digital impulse response is just the sequence of values that emerge from a discrete system when a unit impulse is applied to the input (namely: x=[1 0 0 0 0 0 …]).
If we know the impulse response of a discrete linear system (and it is of finite length) we can readily implement the system in a simple way. To understand this observe that any sampled signal is really just a series of scaled impulses whose amplitudes are the instantaneous amplitudes of the original analogue signal and which occur regularly at the sampling instants. Thus if the input signal is just a series of impulses and we know what the system does to impulses and we know the system obeys the principle of superposition then we can calculate the output of the system to any input signal. We do this by triggering an impulse response at each input sample scaled according to the amplitude of the sample. The sum of the overlapping triggered scaled impulse responses is the output signal. Mathematically this is called convolution, and can be expressed as:
Where the summation is over the length of the impulse response. MATLAB provides a conv() function to perform convolution between two vectors. To apply a linear system of impulse response h[] to an input signal x[], we run:
y = conv(h,x);
When conv() is used in this way the output signal is longer than the input signal by the length of the impulse response (it is as if the input was padded with trailing zeros). Here is a simple example:
h=[10 -9 8 -7 6 -5 4 -3 2 -1];
x=[0 1 0 0 0 0 0 0 0 0 ];
y=conv(h,x);
subplot(2,3,1); plot(h);
subplot(2,3,2); plot(x);
subplot(2,3,3); plot(y);
x=[0 1 0 0 0 0 0 1 0 0 ];
y=conv(h,x);
subplot(2,3,4); plot(h);
subplot(2,3,5); plot(x);
subplot(2,3,6); plot(y);
2. Finite Impulse Response (FIR) Filters
To obtain the impulse response of a system specified in the frequency domain we can use an inverse Fourier transform to get a sampled impulse response then truncate or window the impulse response to a suitable length. The inverse Fourier transform is described in lecture 9.
MATLAB provides a function fir1() for generating a truncated impulse response for the implementation of low-pass, high-pass and band-pass filters. The call:
h = fir1(n,w)
generates a truncated impulse response in h[] of length n samples, corresponding to a low-pass filter with a cut-off at a normalised frequency w. In MATLAB filter design functions, frequencies are normalised to half the sampling rate, i.e. a normalised frequency of 0.5 corresponds to one quarter of the sampling rate (don't ask).
Here we generate 1 second of white noise signal at 10,000 samples/sec and low-pass filter it at 2500Hz with a filter of size 100 samples.
x = rand(1,10000);
soundsc(x,10000);
pause;
h = fir1(100,0.5);
y = conv(h,x);
soundsc(y,10000);
High-pass filters are designed with
h = fir1(n,w,'high');
and band-pass filters are designed with
h = fir1(n,[w1 w2]);
where w1 and w2 are the normalised frequencies of the band edges.
In general you need truncated impulse responses that are long enough to capture the majority of the energy in the full impulse response. Filters with 100 coefficients are not uncommon. Filters with sharp edges need longer impulse responses.
3. Infinite length impulse response filters
Finite impulse response filters have a number of good characteristics: they are simple to understand, have good phase response, and are always stable. However they are inefficient. Each output sample requires a convolution sum that is the size of the impulse response. That is for an impulse response of length N, the system implementation will require N multiply-add operations per output sample.
Increased efficiency can be obtained by allowing the linear system calculation to use not only the current and past samples of the input signal, but also past samples from the output signal. Often the calculation performed by a linear system can be simplified if it can re-use some of the calculation that it performed on previous samples. This is called recursive system design, in contrast to finite impulse response filters which can be called non-recursive. Recursive systems are more complex to build, often have messy phase responses, and can be unstable. However they can be significantly more efficient and have impulses responses of infinite length.
A recursive system is specified by two vectors a[] and b[]: the b[] coefficients are convolved with the current and past input samples, while the a[] coefficients are convolved with the past output samples.
Here is a diagram of the idea:
In this figure the input signal x is processed sample-by-sample into the output signal y. At the stage shown, we are processing sample number n. To calculate output sample y(n), the filter multiplies the current and past input samples x(n), x(n-1), x(n-2), x(n-3), …, x(n-k) by the set of b coefficients: b(0), b(1), b(2), b(3), …, b(k); and sums them. The filter then multiplies the past output samples: y(n-1), y(n-2), y(n-3), …, y(n-k) by the a coefficients: a(1), a(2), a(3), …, a(k) and sums them. It then combines these to form y(n), according to this formula:
In MATLAB, this whole process is performed by the filter() function:
y = filter(b,a,x)
where as before, x is the input signal, y the output signal, and where b and a are the coefficients that define the recursive linear system.
The MATLAB signal processing toolbox contains a number of different functions for deigning recursive low-pass, high-pass and band-pass filters. We shall only refer to one simple design, the Butterworth filter. This design generates filters with maximally flat responses in the pass band, although they may not have as steep a cut-off as other designs.
A butterworth low-pass filter can be constructed with
[b,a] = butter(n,w);
Where n is the number of coefficients required in the filter, and w is the normalised cut-off frequency. Typically n is much smaller for recursive filters than for FIR filters, values of 10 to 20 are common.
A butterworth high-pass filter can be constructed from
[b,a] = butter(n,w,'high');
while a band-pass filter from w1 to w2 can be constructed with
[b,a] = butter(n,[w1 w2]);
In this simple example we filter one second of white noise with a low-pass filter at 2500Hz:
x = rand(1,10000);
soundsc(x,10000);
pause;
[b,a]=butter(10,0.5)
y = filter(b,a,x);
soundsc(y,10000);
4. Frequency Response
To display a graph of the frequency response of a filter, use the function freqz(). For a recursive system with coefficients [b,a], this can be called with
freqz(b,a,n,Fs);
The argument n is the number of points to calculate; 512 is a suitable size. The argument Fs allows you to scale the graph to any given sample rate. For a non-recursive system, where only the impulse response is available, use:
freqz(h,1,512,Fs);
By default, freqz() plots both the magnitude response and the phase response. You can plot just the magnitude response (in decibels) with this code:
[h,f] = freqz(b,a,512,Fs);
plot(f,20*log10(abs(h)));
Here h[] stores the complex response and f[] stores the frequency indices.
Reading
Units 4 and 7 of the "Introduction to Digital Signal Processing" course notes provide a more extended introduction to these topics: http://www.phon.ucl.ac.uk/courses/spsci/dsp/
This book is highly recommended for the student wanting to understand DSP in more detail: P.A. Lynn, W. Fuerst, B. Thomas, "Introductory Digital Signal Processing with Computer Applications", John Wiley, 1997.
Exercises
For these exercises, use the editor window to enter your code, and save your answers to files under your account on the central server. When you save the files, give them the file extension of ".m". Run your programs from the command window.
1. Write a program (ex81.m) that asks the user for a filename (WAV file), then loads the file, low-pass filters it at 1000Hz and replays the original and the filtered version.
2. Write a program (ex82.m) that builds a low-pass filter, a high-pass filter and a band-pass filter and then plots their frequency response curves on the same graph. Use frequency edges of 1000 and 3000Hz and a sampling frequency of 10,000Hz.
3. Write a program (ex83.m) to compare the frequency responses of a 40 coefficient non-recursive implementation and a 20 coefficient recursive implementation of a low-pass filter at a normalised cut-off of one quarter of the sampling rate. Interpret the results.
4. (Homework) Write a simple GUI program for interactive filtering. It should allow the user to choose a filter type and cut-off frequency, then display the filter frequency response. It should also be possible to load a WAV file, then display and replay the unfiltered and filtered signals. This is how it might look:
