How to plot the convolution integral of two functions

Hi, I'm wondering how I can plot out the convolution of two integreals in matlab to see if my hand calculated answer is correct. for this example, we convolve e^(ax) and e^(bx), where a and b are random integers (I used 1 and 2 to keep the math simple), but I am unable to plot the two functions or their convolutions. My code is as follows

x = 0:0.5:10;
y = exp(x);
h = exp(2*x);
q = conv(y, h)
plot(x,q)

Any help is appreciated

Answer 1

[Paul]

Hi Anthony,

When using conv to compute a convolution sum to approximate a convolution integral of two signals (not two integrals), keep in mind that:

a) conv assumes that both signals are equal to zero for values of x greater (less) than x(end) (x(1))

b) The output of conv as used below will satisfy numel(q) = numel(y) + numel(h) - 1

c) The convolution sum has to be muliplied by dx when compared to the convolution integral.

Correcting for (b) and (c) yields (but don't forget (a))

x = 0:0.5:10;  dx = 0.5;  % 0.5 might be too large for this problem?
y = exp(x);
h = exp(2*x);
q = conv(y, h);
plot((0:(numel(q)-1))*dx,q*dx)
xlim([0 10])

Also, you can use the Symbolic Math toolbox to compute a symbolic, closed form expression for the convolution integral

syms

symfun

int

Answer 2

[Chunru]

Use symbolic math:

syms x y h tau
%x = 0:0.5:10;
y = exp(-x);        % use -1 and -2 to ensure convolution exist 
h = exp(-2*x);
%q = conv(y, h)
ytau = subs(y, x, x - tau);     
htau = subs(h, x, tau); 
q = int(htau * ytau, tau, 0, inf)   % this assumes x \in [0 \inf]

q = 

fplot(q, [0 10])

Answer 3

[William Rose]

@Anthony Koning,

Let's use e(-x) and e(-2x), since e(+20) is so big.

N=21; deltax=.5; 
x=(0:N-1)*deltax;
xc=(0:2*N-2)*deltax;  %x values for convolution
y = exp(-x);
h = exp(-2*x);
q = conv(y, h);
%top plot
subplot(311), plot(x,y,'-r.',x,h,'-g.'); 
grid on; legend('y','h'); xlabel('x'); title('y(x),h(x)')
%middle plot
subplot(312); plot(xc,q,'-b.'); 
xlabel('xc'), title('conv(y,h)'); grid on

Try it.

There is no graphical way that I know of to illustrate the convolution function. You can understand the value of the convolution integral at one particular overlap value, i.e. you can attempt to illustrate q(xc) at a specific value of xc, by plotting x and a backwards, translated-by-xc version of h(x). This is not enough, however. The value of q(xc) is the sum of the area under the curve y(x)*h(xc-x). Therefore I have also plotted that curve, in black. In the code and plot below, "hft" stands for "h, flipped and translated".

For xc=2:

hflip=flip(h); 
hft=[hflip(end-4:end),zeros(1,N-4-1)];      %h, flipped and translated
%bottom plot
subplot(313), plot(x,y,'-r.',x,hft,'-g.',x,y.*hft,'-k.')
grid on; xlabel('x'), title('y(x), h(2-x), y(x).*h(2-x)'); legend('y','hft','y.*hft');

The sum of all the values on the black curve in the bottom plot should equal the value of q at xc=2 in the middle plot.

fprintf('q(xc=2)=%.2f, sum(y.*hft)=%.2f.\n',q(5),sum(y.*hft));

q(xc=2)=0.32, sum(y.*hft)=0.32.

Does it? Yes. Multiply both values by deltax, if you want to get area under the curve, instead of simple sum.

To get the next point, i.e. q(xc=2.5) on the blue curve in the middle plot, shift the green curve on the bottom plot one point to the right, i.e. put the green peak at 2.5, instead of 2. Then recalculate the black curve on the bottom plot, and add up the black values on the recalculated bottom plot. The shift again to get the next point on the middle plot. And so on, until you get get all the blue points on the middle plot.