# convolution of signals using matlab

134 views (last 30 days)

Show older comments

The file ‘ekg.txt’ is a human electrocardiogram (EKG) recording which is corrupted with a lot of noise. It was recorded over 30 seconds with sampling frequency = 200 Hz.

a) Load this file into Matlab or python and make a nice plot of it as a function of time on the horizontal axis.

numberOfSeconds = 30

samplingRate = 200 % Samples per second.

numberOfSamples = samplingRate * numberOfSeconds

t = linspace(0, numberOfSeconds, numberOfSamples);

v = load('ekg.txt','-ascii');

% Plot:

plot(t, v);

b) In linear systems analysis, we often use convolution to apply a “filter” (another name for a system’s impulse response) to an input signal. Create your own 1D causal moving average filter of length 25, and convolve it with the noisy EKG signal. Plot the original and resulting signal on the same plot in a 2x1 grid using the subplot command. Compare the two signals. What did the filter do?

windowWidth = 25;% Moving average filter length

smoothy = movmean(v,windowWidth);

plot(t,smoothy);

w = conv(smoothy,v);

plot(w)

windowWidth = 25; % Moving average filter length

kernel = ones(windowWidth,1) / windowWidth;

out = filter(kernel, 1, v);

plot(t,out);

k = conv(out,v,'same');

plot(k);

subplot(3,1,1);

plot(t, v);

subplot(3,1,2);

plot(t,smoothy);

subplot(3,1,3);

plot(w);

subplot(3,1,1);

plot(t, v);

subplot(3,1,2);

plot(t,smoothy);

subplot(3,1,3);

plot(k);

Can i get help on the following questions?

c) Did the filter in b) delay the original signal? By how much? Can you modify your moving average filter so that it still has length 25, but doesn’t delay the signal? Is it still causal?

d) Reverse the order of the convolution you did in c) (i.e., convolve the input with the filter instead of the filter with the input). What happens? What convolution property does this demonstrate?

e) Create a filter to delay the signal by 5 seconds. Convolve with the signal and plot the resulting signal. Write a simple convolution expression for this delayed signal y(t), using x(t) as the input sound file and the mathematical expression for your kernel.

f) Create a time-reversed version of the original signal. Plot and describe the result. Can you invent a filter to create the time-reversed signal by convolution? Why or why not?

##### 0 Comments

### Accepted Answer

Image Analyst
on 28 Sep 2022

@Jiby your code is not correct. To create a causal filter, one which depends on only the current and prior values, you'll need an assymmetric kernel, not just a uniform box filter of all 1's. Remember the convolution flips the kernel so to get a causal filter you'll need to have a kernel that is zeros in the first half and 1's in the second half. Like 12 zeros, then 13 ones.

Next, if you're using conv to do a full convolution (i.e. not using the 'same' option) you get an array that is as long as the sum of the signal length plus kernel length. I don't believe this "delays" the signal -- it merely has a different element that is the time=0 point. Instead of being the first element, it's half a kernel window width into the output vector. So for a 25 element long kernel t=0 happens at index 13 and index 1 is a negative time. I wouldn't say that the signal is delayed if you use a causal filter, it's just that the signal for the t=0 time is a little bit into the array and it has different/filtered values. The signal might look like it was delayed depending on what your kernel is. Let's take the simple case of a delta function, so that convolution gives you back the same array. If the kernel was 25 long, 12 0's then a 1 then 12 more 0s, then the 1 is in the middle of the kernel and there will be no delay whatsoever. However remember that the t=0 time will now be at the 13th element of the output instead of at the first element of the output. Now if you put the 1 not in the middle but to the left or right of the middle it will look like the signal has been delayed or advanced since you'll just get the same signal as the original but starting either before or after element 13 (before or after the t=0 point).

I know it's confusing and I'm sure you'll have to reread this several times, but I hope it helps.

##### 5 Comments

Image Analyst
on 29 Sep 2022

### More Answers (1)

Paul
on 28 Sep 2022

Edited: Paul
on 29 Sep 2022

Hi Jiby,

With respect to part b), you may want to consider:

b) In linear systems analysis, we often use convolution to apply a “filter” (another name for a system’s impulse response) to an input signal. Create your own 1D causal moving average filter of length 25, and convolve it with the noisy EKG signal. Plot the original and resulting signal on the same plot in a 2x1 grid using the subplot command. Compare the two signals. What did the filter do?

The question asks to define a 1D causal filter, apply it to the data, and analyze the results

Here, we're applying movmean to the data. But if you look at the doc page, and/or try some simple examples, you'll see that movmean, as used in the Question, is a non-causal operation. Also, one has to be careful about how it treats samples that are off the edges of the input data. Again, check the doc page and try some simple examples.

windowWidth = 25;% Moving average filter length

smoothy = movmean(v,windowWidth);

plot(t,smoothy);

Even if smoothy is the desired result, it's not clear why that result would again be convolved with the input.

w = conv(smoothy,v);

plot(w)

This part looks very reasonable. Did it come from filter? I think it works as is, but in Matab it's more common to represent a polymomial as a row vector, not a column vector.

windowWidth = 25; % Moving average filter length

kernel = ones(windowWidth,1) / windowWidth;

out = filter(kernel, 1, v);

plot(t,out);

Of course, the problem asks to use convolution, which leads us here. Unfortunately, the 'same' option isn't what you want here. Instead, check out the 'full' option, keeping in mind the result will be longer than the input data (conv 'full' assumes zero-padding as needed), so you'll have to trim the result.

k = conv(out,v,'same');

plot(k);

I suggest you try to work this problem with a much smaller windowWidth, say 3, and a simpler input vector, say v = 1:10, so you can easily compare the results from these different code snippets with a straightforward, "by-hand" calculation.

##### 3 Comments

Paul
on 2 Oct 2022

Edited: Paul
on 3 Oct 2022

If it's still worthwhile to continue the discussion with a simple example ....

Consider a simple filter with output y[n] being the average of the current and previous four inputs:

y[n] = 0.2 * (u[n] + u[n-1] + u[n-2] + u[n-3] + u[n-4])

This filter is causal because the output does not depend on future inputs. Its impulse, or unit pulse, response is found by inspection:

h1[n] = 0.2, 0 <= n < =4

h1[n] = 0, otherwise

Define this filter in a struct:

h1.vals = 0.2*ones(1,5);

h1.n = 0:4;

where it's understood that h1.vals = 0 for discrete-time points n that are not stored in h1.n.

figure

subplot(511)

stem(h1.n,h1.vals);

ylabel('h1[n]')

Let's define a simple input sequence to work with

u[n] = 1+n, 0 <= n <= 9

u[n] = 0, otherwise

u.n = 0:9;

u.vals = 1 + u.n;

subplot(512)

stem(u.n,u.vals);

ylabel('u[n]')

Because h1[n] and u[n] are both finite duration signals, we can use conv to compute their convolution sum (with default option 'full')

y1[n] = h1[n] * u[n]

The values of y1[n] are simply

y1.vals = conv(h1.vals,u.vals);

The corresponding discrete-time points are

y1.n = h1.n(1) + u.n(1) + (0:(numel(y1.vals)-1));

Verifty for one time that h1 does, in fact, perform the desired averaging, and then plot

[y1.vals(y1.n==6) mean(u.vals(find(u.n==6))+(-4:0))]

subplot(513)

stem(y1.n,y1.vals)

ylabel('y1[n]')

Now, let's define a different, 5-point average filter, where the output, y[n], is the average of u[n] and the two preceding and following input samples

y[n] = 0.2*(u[n+2] + u[n+1] + u[n] + u[n-1] + u[n-2])

This filter is clearly non-causal because the output does depend on future inputs. Again, by inspection, its unit pulse response is

h2[n] = 0.2, -2 <= n <= 2

h2[n] = 0, otherwise

As must be the case for a non-causal filter, the unit pulse response is non-zero for some values of n < 0

h2.vals = 0.2*ones(1,5);

h2.n = -2:2;

subplot(514)

stem(h2.n,h2.vals)

ylabel('h2[n]')

Note that h2.vals == h1.vals, but the discrete-time points that correspond to those values are not the same. Again, use conv to compute the values of the output and set the corresponding discrete-time points

y2.vals = conv(h2.vals,u.vals); % obviously, these are the same values as y1.vals

y2.n = h2.n(1) + u.n(1) + (0:(numel(y2.vals)-1));

Verify for a single point and plot

[y2.vals(y2.n==6) mean(u.vals(find(u.n==6))+(-2:2))]

subplot(515)

stem(y2.n,y2.vals)

ylabel('y2[n]')

xlabel('n')

set(get(gcf,'children'),'XLim',[-3 15])

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!