Continuous Convolution Using 3 Different Methods in MATLAB Gives Different Results

Question

I'm trying to understand why the results for the convolution of two functions in MATLAB is different when I'm trying different methods. As an example suppose that my functions are sin(x) and cos(x). The first approach is using the conv() command in MATLAB. The second approach is to calculate it directly using the definition of convolution. The third method is using the convolution theorem and calculating it using fft() and ifft(). The code is attached here. I have used sin(x) and exp(-x) just as an example. I want to do the same for any two arbitrary functions.

x = linspace(0,10,100);
dx = x(2)-x(1);
t = -100:0.01:100;

y = conv(sin(x),exp(-x),'same')*dx;
figure(1)
plot(x,y);

for i=1:length(x)
    zz = sin(t).*exp(-x(i)-t);
    zz = trapz(t,zz);
    z(i) = zz;
end
figure(2)
plot(x,z)

ww = fft(sin(x)).*fft(exp(-x));
w = ifft(ww).*dx;
figure(3)
plot(x,w)

As you can see the sizes of the two inputs are equal. The results from these three approaches are very different. What should I do in order for the 3 approaches to give the same answer? I would appreciate if anyone explained the issues here.

Answer 1

Let’s look at the three methods implemented in the OP:

x = linspace(0,10,100);
dx = x(2)-x(1);
t = -100:0.01:100;

y = conv(sin(x),exp(-x),'same')*dx;

Here the signals you are convolving are defined over the range [0,10], and assumed to be zero outside this range.

for i=1:length(x)
    zz = sin(t).*exp(-x(i)-t);
    zz = trapz(t,zz);
    z(i) = zz;
end

Here the signals are defined over the range [-100,100] for the sine wave, and a larger range for the exponent (the actual range shifts, such that no assumptions are ever being made about what happens outside the range): [-110,100].

Note that the exponential function grows exponentially... at the lower end of the range used its value is 2.6881e+43, a value that dwarfs everything else and dictates the results.

There's also an error: exp(-x(i)-t) should read exp(-(x(i)-t)).

ww = fft(sin(x)).*fft(exp(-x));
w = ifft(ww).*dx;

Here the signals are defined over the range [0,10] as in the first case, but are assumed to repeat outside this range. This is equivalent to convolving [sin(x),sin(x),sin(x)] with [exp(-x),exp(-x),exp(-x)], whereas the first case is convolving [zeros(size(x)),sin(x),zeros(size(x))] with [zeros(size(x)),exp(-x),zeros(size(x))].

You can make the first and last cases equal like this:

a = [zeros(size(x)),sin(x),zeros(size(x))];
b = [zeros(size(x)),exp(-x),zeros(size(x))];
y = conv(a,b,'same');
w = ifft(fft(a).*fft(ifftshift(b)));
plot(y)
hold on
plot(w)

The ifftshift is used here to move the origin of the signal b to the first bin, where the FFT expects it. Without this, the result will be shifted by half the signal length. If you shift a in the same way, w must be shifted in the opposite way, so it's simpler to just not shift it.

These are both discrete-time convolutions. Sampling theory says that, for two band-limited signals, convolving then sampling is the same as first sampling and then convolving, and interpolation of the sampled signal can return us the continuous one. But this is true only if we could sample the functions until infinity, which we can't. Note also that the exponential function is not band-limited. So there are two problems here when relating the discrete and continuous convolutions: the sampling domain and the choice of function.

So let's instead take two signals where we can make these comparisons meaningful. For the first one, a, we can again take a sine wave. For the second one, b, we'll take the derivative of the Gaussian, which is nearly band-limited and is limited in time as well.

x = linspace(-10, 10, 101); % odd number of samples so we sample the origin
dx = x(2) - x(1);

a = @(x) sin(x);
s = 3*dx;
b = @(x) -exp(-x.^2/(2*s^2)) / (sqrt(2*pi)*s^3) .* x;

y = conv(a(x), b(x), 'same') * dx;

t = linspace(-15, 15, 400); % some larger domain, and higher sampling rate, though this is not necessary!
z = zeros(size(x));
for i = 1:length(t)
    z(i) = trapz(t, a(t) .* b(t(i)-t));
end

w = ifft( fft(a(x)) .* fft(ifftshift(b(x))) ) * dx;

plot(x, y, 'DisplayName','conv()')
hold on
plot(t, z, 'DisplayName','"continuous"')
plot(x, w, 'DisplayName','FFT')
legend('-DynamicLegend');

Note how the output of the three forms of convolution are nearly identical in a region in the middle of the plot. Near the edges of the sampling domain they start to differ (the size of this region is given by the size of b, which is why we chose a compact function for it).

[Note also that z is really just a discrete convolution with a higher sampling rate, really, considering that trapz is the same as sum except the two samples at the ends are halved.]