Understanding convolution, 'same' might give wrong result

44 次查看（过去 30 天）
Silke
Silke 2017-11-7
评论： Silke 2017-11-10
I am having a problem with the convolution function. I am trying to convolute two functions, one is bi-exponential decay and the other one a pulse shape. I did it like that:
t = 0:250e-12:1e-6;
a1 = 0.6;
a2 = 0.4;
D1 = 25e-9;
D2 = 350e-9;
filename1 = 'pulse410.txt';
pathname1 = 'D:\My Documents\PulseShapes\';
pulse410 = importfile1([pathname1, filename1]);
[P,I] = max(pulse410(:,1));
pulse410s=spline(pulse410(:,1), pulse410(:,2), t); % pulse410s contains the same amount of data points as t.
t_in = find(t==P);
for ii = t_in:numel(t)
pulse410s(ii)=0;
end
f1 = conv(a1 * exp(-t/D1) + a2 * exp(-t/D2), pulse410s./max(pulse410s), 'same') ;
In the picture below, data 1 is the bi-exponential decay, data 2 is pulse410s, data3 is the convoluted signal. Well, I would expect the convolution of the two signals to look different. Especially, the drop at 0.5e-6 looks weird to me. Is this weird result due to me using the 'same'-option in the convolution?
How can I obtain a good convolution by keeping the amount of data points the same before and after the convolution? Thanks for your help!

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

回答（2 个）

David Goodmanson
David Goodmanson 2017-11-7
Hi Silke,
Yes, this is a direct result of using the 'same' option. If you take a look at the convolution without that option, all 'same' is doing is taking the middle section, between 1/4 and 3/4 of the x axis and throwing away the rest. But some of the rest has needed information.
When you do a convolution of arrays with n1 and n2 points, the result has n1+n2+1 points. Right now you are zerofilling the pulse to give it the same number of points as the exponentials, 4001 points. There is no need to do that. It appears that the pulse has a nonzero extent of about 65 points after interpolation. If you keep an array of that size and convolve it with the exponentials, you will have an array of about 4067 points or so with the same information.
Starting from somewhere in negative time, the exponentials have a sharp rise at t=0 which is real, and then a sharp drop at the end because your time record runs out. The second drop is an artifact. When you do the convolution to get 4067 or whatever points, the part at the beginning is real and the part at the end is not. So if you want the convolution to have the same length as the signal, the best approximation is to take the first 4001 points in the convolution.
  3 个评论
显示 1更早的评论
Eric
Eric 2017-11-8
I think you expect the answer to be the continuous convolution of the two functions, namely the one involving integrals and the like. What MATALAB actually does is a discrete convolution, which is essentially the same thing, but effectively drops the integral term dt from the definition. If you want to approximate the continuous convolution, you can multiply MATLAB's discrete convolution by your time step, namely
f1 = conv(....) .* dt
where dt=diff(t(1:2)) or something of the like. As for all estimates, the finer your time step, the closer your value will be to the continuous value.
Silke
Silke 2017-11-8
Yes, this is actually what I was looking for. But I am not sure if it is really required for my goal. I am quite puzzled myself at the moment. Well, I have put the convolution problem now in a bigger context (fitting a function). If you have time/interest, it would be great if you could have a look at it. https://nl.mathworks.com/matlabcentral/answers/365865-nlinfit-gives-no-result Thank you!

请先登录，再进行评论。

Eric
Eric 2017-11-7
Your convolution looks correct, at least for what you have done, and yes, it is partially because of the 'same' option. The other part of the reason you are seeing the behavior you are is due to your extending the pulse waveform by appending zeros to it.
Let us assume I have a pulse that is 10ms long and my decay is simply e^-t. Here are several different ways to use MATLAB's convolution function:
t = 0:0.001:1;
f0 = conv(exp(-t),[ones(10,1)]); %The default, which uses 'full'
f1 = conv(exp(-t),[ones(10,1)],'same'); %Just the pulse, 'same'
f2 = conv([ones(10,1)],exp(-t),'same'); %Just the pulse, 'same' but inputs reversed
f3 = conv(exp(-t),[ones(10,1); zeros(length(t)-10,1)],'same'); %Pulse appended with zeros, 'same'
f4 = conv(exp(-t),[ones(10,1); zeros(length(t)-10,1)]); %Pulse appended with zeros, 'full'
f5 = conv(exp(-t),[ones(10,1)],'valid'); %Just the pulse, 'valid'
What you will find when comparing f1 and f2 to f0, or f3 to f4, is that the argument 'same' basically takes the center n points from the full convolution, where n is the length of the first vector. This is expected and how it is defined in the documentation. Now, when comparing the plots of f0 and f4, you can see that the 'full' result of f4 is also extended by a number of zeros relative to f0. As such, when you call for 'same' using the f4 method, your center points will be different than if you had not appended the zeros (f0). MATLAB will do any zero padding for you internally, so there is no need to do it yourself.
P.S. There is also the option of 'valid', (f5) which may interest you, which will return only those results where it did not need to use zero-padding. Basically, It will remove the sharp drops at the edges of the convolution. If you are concerned about it, the output vector's length can be obtained based on the lengths of the inputs (so adding zeros will change this). See conv:shape for more info.
  8 个评论
显示 6更早的评论
David Goodmanson
David Goodmanson 2017-11-9
编辑：David Goodmanson 2017-11-9
Hi Silke,
Clearly you have to multiply by the array spacing dx when you are trying to approximate a continuous integral, convolution being one possible example, Int{-inf,inf} f(x)g(y-x) dx. In lots of cases f and g are algebraic functions so you can make the x array twice as fine to improve the accuracy. Twice as many points in the 'active' region doubles the integral but the new dx is half as large and compensates.
In signal processing there is a sample frequency which is taken as a given, you are not playing around with the width of the intervals and within that mileau there is no disadvantage to setting the width to be 1. One can do all the calculations like the one Eric showed above without any reference to what the sampling frequency actually is.
Silke
Silke 2017-11-10
Hi Eric and David,
thanks for your reply. Well, I am actually interested in getting an answer with the correct units, as I am trying to fit measurement data. So I really believe that I need to multiply with dx. Thanks a lot for your help!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by