How to Compute Convolution Using FFT in Matlab?

[bjnartowt]

Homework Statement

You have two functions in Matlab (represented as column vectors). Compute their convolution using the fast Fourier transform.

Homework Equations

The Attempt at a Solution

I am having trouble finding a book with this topic. I would like to know where the Matlab-code for such an algorithm would be discussed. It does not seem to be treated in my Giordano/Nakanishi computational text, nor is it in a text on pure Matlab programming. I just need someone to tell me a book.

 

[jbunniii]

Ignoring Matlab for a moment, do you know how to compute a convolution using an FFT, theoretically?

 

[bjnartowt]

Yeah--I know how to theoretically a convolution.See attached.

I also know that "non-theoretically", the domain over which you convolve is not infinite, so you need to zero-pad. I also know you have to "shift" the convolution.

I have written some code already, and little gremlins keep making it not work. Just kidding...but that's what it seems like. I was looking for a text to help me figure out how to do it correctly.

 

[bjnartowt]

I know how to do it, and I have code that, in theory, should work properly. The code is behaviing in such a bizarre manner that I am almost certain I need to consult books about the "theory" of convolving-by-fft in order to grasp what the bug is.

 

[MisterX]

Would you be able to explain what you think is bizarre? The result of ifft(fft(f).*fft(g)) should be the same as cconv(f, g, length(f)) (sizes of f and g must be the same). That is, it will be the same as circular convolution, modulo the length of the vectors.

It will not be the same as conv(f, g), since that would be a vector of length length(f)+length(g)-1.

 

[bjnartowt]

I took two boxcar functions, f=bxc(x,-1,1); g=bxc(x,-2,2); (f has a value of 1 from -1 to 1, and g has value of 1 from -2 to 2; they are 0 everywhere else). The analytic form convolution can be computed, integrate(f(x')*g(x-x')), from x integrated over all domain, yields,

analytic=-(3+xp).*heaviside(-3-xp)+(1+xp).*heaviside(-1-xp)-(-2+xp).*heaviside(2-xp)+xp.*heaviside(-xp);

You have the vectors x = -1:0.01:1; xp = -2-0.005:0.01:2+0.005; where the vector xp is double the length of vector x, because the padding doubles the length of the convolved function.

The program produced the following plot [see attached .png], where "result" is circularly shifted,
plot(xp,circshift(result',200)',xp,circshift(analytic',0)');

-------------------------------------------------------------
clear; clc;
bxc = @ (x,xmin,xmax) heaviside(xmax-x)-heaviside(xmin-x);

x = -1:0.01:1; xp = -2-0.005:0.01:2+0.005;
L = length(x) -1;

% f = sin(pi*x); g = [ zeros(1,L/2) 1 zeros(1,L/2) ];
f=bxc(x,-1,1)*1.2; g=bxc(x,-2,2)*1.2;
analytic=-(3+xp).*heaviside(-3-xp)+(1+xp).*heaviside(-1-xp)-(-2+xp).*heaviside(2-xp)+xp.*heaviside(-xp);

fpad = [ zeros(1,L/2) f zeros(1,1 + L/2) ];
gpad = [ zeros(1,L/2) g zeros(1,1 + L/2) ];

fpadfft = fft(fpad);
gpadfft = fft(gpad);
prodfft = fpadfft.*gpadfft;
result = ifft(prodfft)*0.01;
% result2 = ifft(gpadfft.*gpadfft);

indexi = 1:(2*L + 2);
shiftzero = 1 + mod(indexi -1 +L,(2*L + 2));

plot(xp,circshift(result',200)',xp,circshift(analytic',0)');

 

[bjnartowt]

I am sorry--I see at least one mistake. My analytic convolution was taken for boxcar functions that run to intervals as wide as -2 to 2, but the convolution I effect numerically is over the interval -1 to 1. On remedying this problem, one convolution-function still appears stretched out with respect to the other.

 

[bjnartowt]

Anyway--I don't expect someone to troubleshoot my code. I just need the name of a book in which this is discussed and treated pedagogically.