Understanding convolution, 'same' might give wrong result

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!

Answer 1

[David Goodmanson]

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.

[Silke]

Hi David,

Thanks for your help. I am still struggling with the convolution. I have done now this:

[P,I] = max(pulse410(:,1));

   t_in = find(t==P);
   t_x = t(1:t_in);
   pulse410s=spline(pulse410(:,1), pulse410(:,2), t_x);

f1x = conv(a1 * exp(-t/D1) + a2 * exp(-t/D2), pulse410s./max(pulse410s)) ;

But still, I do not get the result.

In the figure, data1 is the bi-exponential decay, data2 is the pulse, and data4 the convoluted signal.

I would not expect the convolution of two signals starting at 1 ending up at ~12. Is my expectation wrong, or is there still something wrong with the code.

[Eric]

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]

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!

Answer 2

[Eric]

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.

[Silke]

Dear Eric, thank you for your detailed reply. I also tried the valid function on my data, but I noticed that I then lose the raise of the function after 0, so the signal peaks at 0.

I still do not have it clear how Matlab does the convolution, as I would expect different results (see my comment above).

[David Goodmanson]

Hi Silke,

I don't know if you are using the visual picture of convolution (ignore this if so) but I find it's quite helpful.

[1] take the functions and flip one of them left-to-right. I suggest flipping the pulse function since it's narrow anyway, but it doesn't matter. The two functions definitely do not have to have the same number of elements.

[2] put the two functions side-by-side so that they do not overlap, with the flipped one on the left. (If the functions contain zeros, those can not overlap either).

[3] start sliding the functions across each other, and take the sum of the products of the overlap elements. The first element of the convolution is when they overlap by one element. The last element of the convolution is when they overlap by one element when you are almost done sliding them by each other.

[4] As Eric has alluded to, adding a bunch of zeros to either end of either function does not change any of the nonzero convolution values, because multiplying by a bunch of extra zeros and adding that in doesn't do anything. But it does change the length of the convolution array, because convolution starts when they overlap, even if all that is overlapping is zeros for either or both of them.

[Silke]

Hi David,

thanks for your reply. Yes, how you describe the convolution is how I learned it roughly 15 years ago (and not having used it since). I agree that the number of elements does not have the same for both functions, but the time-step needs for sure be equal. I think I just got overexcited while programming, so instead of just making the time-step equal, I made the number of elements the same. Also, I did not really get when it is necessary to multiply the convolution by the timestep, and when this is not necessary. Or better to say, why is it not always necessary to multiply by the timestep?

[Eric]

When a function is sampled, you move into the sampling domain where the distance between samples is 1. In other words, the classic discrete convolution is actually

<<latex-codecogs-com-gif-latex-w-n--space---space---sum_-m----infty-----infty--space-u-m--space-v-n-m--space.dm>>

but since dm is always 1, it disappears. The limits of inifinity disappear when you work with pulses, or finite vectors, so you are left with MATLAB's convolution explanation,

<<latex-codecogs-com-gif-latex-w-k----sum_-j-u-j-v-k-j-1.>>

which is basically just taking the convolution where u and v are valid, and assuming zeros/ignoring summation terms where u or v is needed but does not exist. (Note: MATLAB is not actually zero-padding, but by limiting what it adds, it has the same effect as if MATLAB had padded with zeros when/where needed.)

[Silke]

Yes, but do I have to multiply with the timestep then to get the "real" amplitude of the convolution, or not?

[Eric]

It depends on your application and what meaning you want your convolution to have, since you are basically asking about whether to multiply the entire answer by a scalar. The timestep you are referring to is actually a conversion between time domain and sample domain (e.g. 1ms/sample)

Add the time-sample conversion if:

• You want your convolution to have units involving seconds rather than samples (because of the math or just because)
• You are ultimately trying to estimate the continuous convolution's value through MATLAB

You do not need to convert if:

• Your math is okay with or actually expects you to be using sample space
• You only care about the form/shape of the answer (i.e. the graph) and not the exact values
• You have no preference

Based on the fact that you were not expecting the values you got, it seems you might have actually wanted the continuous answer, but I don't know what you need it for. Look at your application and decide. If you are okay with sample space, just keep the conversion in mind for the future if you ever want to compare the two answers again.

[David Goodmanson]

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]

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!