Time-domain filtering in Matlab

 

Time domain filtering in Matlab using “convolution”.

 

This tutorial is about the Matlab “conv” command.

To find out what “conv” does, you could try reading the Matlab online help.

 

Try typing

help conv

and check out the explanation.

(You may find it easiest to copy the command lines in italics to the Matlab comand window using copy and paste.)

However, unless you know a lot about engineering maths, you will probably find the Matlab “help” much more “convoluted” than “helpful”. In the unlikely event that you don’t have that problem, carry on with the rest of these exercises (otherwise, either see whether one of the people you fancy in this class looks lost and could benefit from your insights and your charm, or just go to the beach for a little while).

Now, “convolving” is something that we do to “signals”, and in Matlab signals are vectors. It’s often a good idea to play with very little artificial signal vectors first to understand what the signal processing commands do, before launching them onto big “real” signals.

 

In Matlab, let’s make a toy “data” vector d with the following values:

d=[2 1 2 1];

then try the following

conv(d,1)

What do you get?

conv(d,0.5)

What do you get?

conv(d,-1)


Convolving a vector with a “scalar” (i.e. a simple number or vector with just one element), “scales” the vector. (Now you know why the call them scalars).

 

Now for something a little more mysterious. Try each of these in turn:

conv(d,[1 0 0 0 0 0])

conv(d,[0 1 0 0 0 0])

conv(d,[0 0 1 0 0 0])

conv(d,[0 0 0 1 0 0])

conv(d,[0 0 0 0 1 0])

conv(d,[0 0 0 0 0 1])

 

What’s going on here? Think of a vector that is all zeros except for a single one as a “unit impulse”. Now, if we “convolve” a signal with a unit impulse, then we don’t change the signal (except, perhaps, by “zero padding” the end of it), but we may delay it (effectively by sticking zeros in front of it).

 

So can you guess what would happen if you try this?

conv(d,[0 0 2 0 0 0])


Think about it first, then try it out to see if you were correct.

What you should have found is that you can combine the scaling and the delaying properties of the conv command, you should have signal twice as big as d, delayed by to samples (or “taps”).

 

Before we can understand what all of this has to do with filtering we need one more thing.

 

Try this:

a=[1 0 0 0 0 0]

b=[0 0 0 0 0 1]

conv(d,a)

conv(d,b)


Nothing new here, we just assigned our “impulses” of varying delay to variables

But now, if I do

a+b

what do I get?

 

And what, therefore, of

conv(d,a+b)


What you should see is that a+b has two impulses, and that the convolution therefore contains two copies of the signal d, one with zero delay, and one delayed by six taps.

 

If we make the delay smaller and smaller, we can make the tow copies of the signal “collide”, or superimpose.

 

Try each of these in turn:

conv(d,[1 0 0 0 0 1])

conv(d,[1 0 0 0 1 0])

conv(d,[1 0 0 1 0 0])

conv(d,[1 0 1 0 0 0])

conv(d,[1 1 0 0 0 0])


You should see that, as the delayed copies of the signal “collide”, they add.

In summary, convolution combines scaling, delaying and adding in one handy package.

 

Now what has all of this to do with filtering?

Let’s consider a relatively simple example: the moving average filter. Imagine I want to obtained a filtered signal f by two-point moving average filtering on vector
d= [d(1) d(2) d(3) ...]. What this means I want the first point of f to be the average of the first two points of d, i.e.  f(1)=d(1)/2 + d(2)/2. Similarly I want f(2)=d(2)/2+d(3)/2 and so on until f(n)=d(n)/2+d(n+1)/2. If you think about it, you may notice that f is simply the sum of two copies of d, with the second copy delayed by one tap, and both copies scaled by one half, and rather than doing this one element at a time, we let conv do it for us

 

f=conv(d,[0.5 0.5])

 

will give you the two point moving average filtered (smoothed) version of d.

 

Now, supposedly, this smoothing can help reduce “noise”.

 

Let’s try this. First, let’s make a nice sinusoidal signal:

s=sin(pi/10:pi/10:10*pi)*3;

plot(s)

legend([{'clean'}]);


Now we contaminate the signal with some ugly gaussian noise

n=randn(size(s));

sn=s+n;

hold on;

plot(sn,'r')

legend([{'clean'},{'noisy'}]);


Will smoothing with conv recover the clean signal? Let’s try it:

plot(conv(sn,[0.5 0.5]), 'g')

legend([{'clean'},{'noisy'},{'filtered'}]);


Well, looking at the plot it is clear that the two-point moving average filter improved matters only a tiny bit. How about a “more aggressive” filter?

Try five-point moving average

plot(conv(sn,[0.2 0.2 0.2 0.2 0.2]), 'm')

legend([{'clean'},{'noisy'},{'filtered'},{‘heavily filtered’}]);


(By the way, for n-point moving average filters longer than n=5, writing them out explicitly becomes tedious, and you will find it easier to write ones(1,n)/n rather than the explicit [1/n 1/n 1/n ..... 1/n])

 

That has removed a fair bit of the noise (by no means all of it, but it’s not straight-forward to do much better, which is why it’s good to try to have as little noise as possible to begin with!) but if you look carefully you will see that the filtered signal is actually time shifted relative to the original signal (and it is a little longer too!). Electrical engineers would say that our five tap moving average filter had a “group delay” of two taps, and given the earlier exercises with delayed impulses you may have a bit of an intuition as to where these delays come from. These delays can be a nuisance. Imagine you wanted to measure the times when the peaks occur. Filtering would make the peaks cleaner, but would delay them. You could get rid of the delays by throwing away the first n/2 points of the n-point moving average filtered signal.

 

Matlab has actually a clever variant on conv called filtfilt which avoids the delays by running the filtering “twice”, once “forward” then “backward”. The filtering is then twice as “effective” too though. Try this:

 

clf;

hold on

plot(s)

plot(filtfilt([0.2 0.2 0.2 0.2 0.2],1,sn), 'm')


There is no delay and the filtered signal is quite smooth, but it’s actually beginning to be reduced in amplitude. The filter is starting to eat into the signal as well as the noise. Perhaps five-point moving average is too much if applied twice.

 

How much filtering is “too much”?

There is no single correct answer to this, and to get a proper appreciation of how much a particular filter will affect which frequencies requires an understanding of signal processing that is beyond these simple introductory exercises.

 

You can nevertheless get a bit of an intuitive understanding if you consider the following.

Earlier you saw that the convolution f=conv(s,a), where a is n taps long, means that f(i)=s(i)*a(1)+s(i+1)*a(2)+s(i+2)*a(3)+...+s(i+n-1)*a(n).  So f(i) is the dot product or the inner product of the vectors a and s(i:i+n-1). If you have studied linear algebra you may remember that the dot product measures the “projection” of s(i:i+n-1) onto a. Even if you haven’t studied linear algebra, if you think about it you will probably see that f(i) is effectively a measure of how well the signal segment  s(i:i+n-1) is correlated with the filter a. So the i-th time point of the filtered signal f(i) tells us how well the original signal segment s(i:i+n-1) matches the waveform of the filter a. In the examples we have seen so far, our moving average filters were all straight line segments. Clearly, fast fluctuations in a signal match a straight line more poorly than slow fluctuations, so the straight line segments we used in running average filtering are examples of low pass filters. They are a bit like flat “irons”, used to iron out small wrinkles in the signal. And if you think about them that way, it should become intuitively clear that:

 

an n-point moving average filter will strongly suppress fluctuations whose frequencies have a wavelength less than n, but will have very little effect on frequencies with wavelengths much larger than n.

 

Convolution can be used not just for low pass filtering. If we make a filter that contains lots of rapid fluctuations (high frequencies) then it will amplify the rapid fluctuations and suppress the low ones.

 

Let’s try that. Make an n-point “high pass filter”:

 

a=[1 -1 0 1]/3;


The filter a now oscillates as fast as possible at the current sampling rate.

 

What will it do:

clf;

hold on;

plot(sn,'r')

plot(filtfilt(a,1,sn), 'g')

legend([{'noisy signal'},{'high passed'}]);

hold off


What you should see is that the amplitude of the sinusoid is much reduced, but the rapid fluctuations “sort of” remain (although they are not unaffected!)

 

How did I arrive at a=[1 -1 0 1]/3? A blind guess, which is not something you’d ever do if you “seriously” filtered signals for scientific research (this was merely a “proof of principle”). To make “real” filters that you can use to high-pass, band-pass or low-pass filter “real” signals, Matlab offers you a handy command called “fir2”.

 

Let’s illustrate its use with some real auditory cortex data. (You may need to download the auditory.mat file below and copy it to your Matlab working directory)

load auditory 

Fs=12000; % sample rate

period=1/Fs;

v=WBNa_0_1current(1,:);

time=(1:numel(v))*period;

plot(time,v)


Just eyeballing the signal, we see large fluctuations on a scale of about 1/20th of a second, probably due to synaptic inputs, as well as a small amount of high frequency noise.

 

Say we want to low-pass filter this signal so that all frequencies below 1 kHz are unaffected, but higher frequencies are suppressed. The fir2 command will give us such a filter, but it needs two vectors, A and F, which tell it what Amplitudes we want at which Frequencies. Unfortunately we can’t just specify F in kHz, but have to give it as a fraction of the “Nyquist” frequency, which is equal to half our sample rate. In this example, the sampling rate is 12000, so if we want to keep our frequencies up to 2kHz flat, then we want A=1 for all 0<F<2000/6000 and A=0 for 2000/6000<F<1. For theoretical reasons, it is not a good idea to try to make the filter too “steep”, so rather than saying F=[0 2000/6000 2001/6000 1]; A=[1 1 0 0]; we would be better advised to use

 

F=[0 1950/6000 2100/6000 1]; A=[1 1 0 0];


Now let’s make a 20 point “finite impulse response” filter for these values and use it to filter the data

 

flt=fir2(20,F,A);

v_filtered=filtfilt(flt,1,v);

plot(time,v); hold on; plot(time, v_filtered, 'r'); hold off;


Use the zoom command in the figure window to plot in.

 

Experiment with different cut-off frequencies.

How might you use the fir2 command to make a high-pass filter?

How might you make a band-pass filter?

 

Try it!

AttachmentSize
auditory.mat1.04 MB