convolution of signals using matlab

67 次查看（过去 30 天）
Jiby
Jiby 2022-9-27
编辑： Paul 2022-10-3
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?

请先登录，再回答此问题。

采纳的回答

Image Analyst
Image Analyst 2022-9-28
@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 个评论
显示 3更早的评论
Jiby
Jiby 2022-9-28
Ran in:
How can i create assymetric filter?
p=1:12; % define space
q=13:25;
delta2 = [ones(size(p)) zeros(size(q))];
delta2
delta2 = 1×25
1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
stem(delta2);
Image Analyst
Image Analyst 2022-9-29
Yes, that would be an asymmetric filter. If you use conv, the signal will be a blurred version of the future values (opposite of causal). Remember conv() flips your kernel.

请先登录，再进行评论。

更多回答（1 个）

Paul
Paul 2022-9-28
编辑：Paul 2022-9-29
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 个评论
显示 1更早的评论
Paul
Paul 2022-9-29
You're welcome.
I don't understand the discussion by @Image Analyst regarding the need for an assymetric kernel. I believe the task in part (b) can be accomplished by a very simple filter (as you originally tried), as illustrated in the doc page for filter in the very first example. However, the output from filter will only be as long as the input sequence, which would be an incomplete result (because there are additional, non-zero values in the output of the filter). Besides you're asked to use convolution. So you need to find the impulse response of the filter, which should be very easy for this problem (and you've already done), and then use conv with 'full' option to get the complete response you seek. But again, I don't see how any zero-padding, or assymetical impulse response, comes into play for this problem.
Paul
Paul 2022-10-2
编辑：Paul 2022-10-3
Ran in:
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))]
ans = 1×2
5 5
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))]
ans = 1×2
7.0000 7.0000
subplot(515)
stem(y2.n,y2.vals)
ylabel('y2[n]')
xlabel('n')
set(get(gcf,'children'),'XLim',[-3 15])

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by