# Computing a convolution using FFT

I have two sequences of the same length, $(x_i), i=1, 2, \ldots, N$ and $(y_i), i=1, 2, \ldots, N$ and a function $K(t) = -t \times \exp(-t^2 / 2)/ \sqrt{2 \pi}$.

I need to compute the following quantity for each $m=1, 2, \ldots, N$:

$\sum_{j=1}^N K(x_m - x_j) \times y_j$

which is a tad slow when done directly (I need it when $N = O(10^4)$ ). I know this can be much improved on with the use of FFT, however it is something I never really worked with. Could anyone suggest any links or a way to rewrite it in FFT form?

I'd be very grateful for any help :)

-
Apply the FFT two both sequences, multiply the results pointwise, then transform back. –  dls Mar 16 '12 at 16:11
Well that's the general idea. Except that to do it for each $m$ is no time gain. I was guessing the clever thing to do is to approximate the function $\sum_{j=1}^N K(s - x_j) y_j$ using FFT, and this bit I'm not sure how to do. Besides, the FFT theorem convolution applies directly to convolutions like $\sum a_{i-j} b_j$, so to do the straightforward thing I would need to recompute the sequence $K(x_m-k_j)$ for each $m$, which would eat away any time gain..? –  Maciej Mar 16 '12 at 16:20
I've not heard of approximating a function by its Fourier transform. You have $$(K \ast y) (x_m) = \sum_{j = 1}^N K(x_m - x_j) y_j$$ and $$\widehat{K \ast y} (x_m) = \hat{K}(x_m) \hat{y_m}$$ You can use FFT to compute these: do the multiplication and then do the inverse transform for each $m$. I don't understand how to use the Fourier transform to approximate a function. The other thing you could think about is to numerically approximate it: compute only a few of the function values and then use interpolation. –  Matt N. Mar 17 '12 at 8:32
I'm not sure if this is precisely what you are talking about but naive convolution is $O(n^2)$, while the FFT method is only $O(n \log n)$ (it is 2 FFTs, which are $O(n \log n)$, and one pointwise multiply, which is $O(n)$, so the process is $O(n \log n)$ overall.) –  dbaupp May 23 '12 at 8:34

Here you have a good link.

Basically you need to compute the FFT of each signal individually, multiply the spectrums and the do the inverse FFT of the resulting sectrum. That will be the result of the convolution.

Good luck!

-

Need to account for wrap around issues

1) Shift the K sequence such that it is centered in the middle of the array call K1 of lenght n1 2) Store the K1 array in a new array K2 with double the size of K1 n2 = 2Xn1 such that K2(1 to n1 ) = complex(K1,0) and K2( n1+1 to n2 ) = 0 where K2 is complex 3) Next shift the K2 array over to the left by n1/2-1 this will center the maximum of the K function in the first bin of the K2 array 4) Store the y data in an array y2 of of lenght n2 such that y2(1 to n1 ) = complex(y,0) and y2(n1+1 to n2) = 0 where y2 is complex

5) calculate the FFT of K2, i.e., FFT_K2 = FFT( K2, 1 ) 1 => exp(iwt) 6) calculate the FFT of y2, i.e., FFT_y2 = FFT( y2, 1 ) 1 => exp(iwt) 7) Let FFT_d2 = CONJUGATE(FFT_K2) * FFT_y2 8) Inverse FFT, i.e., d2 = FFT( FFT_d2, -1 ) -1 => exp(-iwt) 9) convolution = REAL( ds(1 to n1) ) 10) this assumes the exp(-iwt) has the 1/n factor

-