{"id":297,"date":"2009-11-03T14:56:58","date_gmt":"2009-11-03T19:56:58","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2009\/11\/03\/the-conv-function-and-implementation-tradeoffs\/"},"modified":"2019-10-29T08:03:07","modified_gmt":"2019-10-29T12:03:07","slug":"the-conv-function-and-implementation-tradeoffs","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2009\/11\/03\/the-conv-function-and-implementation-tradeoffs\/","title":{"rendered":"The conv function and implementation tradeoffs"},"content":{"rendered":"
\r\n

A friend from my grad school days (back in the previous millenium) is an electrical engineering professor. Students in his\r\n class recently asked him why the conv<\/tt> function is not implemented using FFTs.\r\n <\/p>\r\n

I'm not on the team responsible for conv<\/tt>, but I wrote back with my thoughts, and I thought I would share them here as well.\r\n <\/p>\r\n

Let's review the basics. Using the typical convolution formula to compute the one-dimensional convolution of a P-element\r\n sequence A with Q-element sequence B has a computational complexity of . However, the discrete Fourier transform (DFT) can be used to implement convolution as follows:\r\n <\/p>\r\n

1. Compute the L-point DFT of A, where .\r\n <\/p>\r\n

2. Compute the L-point DFT of B.<\/p>\r\n

3. Multiply the two DFTs.<\/p>\r\n

4. Compute the inverse DFT to get the convolution.<\/p>\r\n

Here's a simple MATLAB function for computing convolution using the Fast Fourier Transform (FFT), which is simply a fast algorithm\r\n for computing the DFT.\r\n <\/p>

type conv_fft<\/span><\/pre>
\r\nfunction c = conv_fft(a, b)\r\n\r\nP = numel(a);\r\nQ = numel(b);\r\nL = P + Q - 1;\r\nK = 2^nextpow2(L);\r\n\r\nc = ifft(fft(a, K) .* fft(b, K));\r\nc = c(1:L);\r\n\r\n<\/pre>
Note that the code uses the next power-of-two greater than or equal to L, although this is not strictly necessary. The fft<\/tt> function operates in   time whether or not L is a power of two.\r\n   <\/p>\r\n   
The overall computational complexity of these steps is  . For P and Q sufficiently large, then, using the DFT to implement convolution is a computational win.\r\n   <\/p>\r\n   
So why don't we do it?<\/p>\r\n   
There are several technical factors. Let's look at speed, exact computation, and memory.<\/p>\r\n   
Speed<\/b><\/p>\r\n   
One factor is that DFT-based computation is not always<\/b> faster.  Let's do an experiment where we compute the convolution of a 1000-element sequence with another sequence of varying\r\n      length.\r\n   <\/p>
x = rand(1, 1000);\r\nnn = 25:25:1000;\r\n\r\nt_normal = zeros(size(nn));\r\nt_fft = zeros(size(nn));\r\n\r\nfor<\/span> k = 1:numel(nn)\r\n    n = nn(k);\r\n    y = rand(1, n);\r\n    t_normal(k) = timeit(@() conv(x, y));\r\n    t_fft(k) = timeit(@() conv_fft(x, y));\r\nend<\/span><\/pre>
plot(nn, t_normal, nn, t_fft)\r\nlegend({'Normal computation'<\/span>, 'FFT-based computation'<\/span>})<\/pre> 
For sequences y<\/tt> shorter than a certain length, called the cross-over point<\/i>, it's quicker to use the normal computation.\r\n   <\/p>\r\n   
Exact computation<\/b><\/p>\r\n   
A second consideration is whether the computation is subject to floating-point round-off errors and to what degree.  There\r\n      are applications, for example, where the convolution of integer-valued sequences is computed.  For such applications, a user\r\n      would reasonably expect the output sequence to be integer-valued as well.\r\n   <\/p>\r\n   
To illustrate, here's a simple function that computes n-th order binomial coefficients<\/a> using convolution:\r\n   <\/p>
type binom<\/span><\/pre>
\r\nfunction c = binom(n)\r\n\r\nc = 1;\r\nfor k = 1:n\r\n    c = conv(c, [1 1]);\r\nend\r\n\r\n<\/pre>
binom(7)<\/pre>
\r\nans =\r\n\r\n     1     7    21    35    35    21     7     1\r\n\r\n<\/pre>
I wrote a variation called binom_fft<\/tt> that is the same as binom<\/tt> except that it calls conv_fft<\/tt>.\r\n   <\/p>
format long<\/span>\r\nbinom_fft(7)<\/pre>
\r\nans =\r\n\r\n  Columns 1 through 4\r\n\r\n   1.000000000000000   6.999999999999998  21.000000000000000  35.000000000000000\r\n\r\n  Columns 5 through 8\r\n\r\n  35.000000000000000  21.000000000000000   7.000000000000000   0.999999999999996\r\n\r\n<\/pre>
Whoops! What's going on? The answer is that the FFT-based implementation of convolution is subject to floating-point round-off\r\n      error.\r\n   <\/p>\r\n   
I imagine that most MATLAB users would consider the output of binom_fft<\/tt> to be wrong<\/b>.\r\n   <\/p>\r\n   
Memory<\/b><\/p>\r\n   
The last technical consideration I want to mention is memory.  Because of the padding and complex arithmetic involved in the\r\n      FFT computations, the FFT-based convolution implementation requires a lot more memory than the normal computation. This may\r\n      not often be a problem for one-dimensional computations, but it can be a big deal for multidimensional computations.\r\n   <\/p>\r\n   
Final thoughts<\/b><\/p>\r\n   
The technical considerations listed above can all be solved in principle.  The implementation could switch to using the normal\r\n      method for short or integer-valued sequences, for example. And there are FFT-based techniques such as overlap-and-add<\/a> to reduce the memory load.\r\n   <\/p>\r\n   
But the problem can get quite complicated.  Testing floating-point values to see if they are integers, for example, can be\r\n      slow.  Also, I suspect (but have not checked) that a multithreaded implemention of the normal computation will take advantage\r\n      of multiple cores better than a multithreaded FFT-based method. That would change the cross-over point between the two methods.\r\n      And the exact cross-over point will vary from computer to computer, making it likely that our implementation would be somewhat\r\n      slower for some people for some problems.\r\n   <\/p>\r\n   
Overall, my inclination would be to provide FFT-based convolution as a separate function rather than reimplementing conv<\/tt>.  And that's what we've done. See the function fftfilt<\/tt><\/a> in the Signal Processing Toolbox.\r\n   <\/p>\r\n   
Do you disagree with this approach?  Post your comment.<\/p>