Writing a Discrete Fourier Transform program

I would like to write a DFT program using FFT. This is actually used for very large matrix-vector multiplication (10^8 * 10^8), which is simplified to a vector-to-vector convolution, and further reduced to a Discrete Fourier Transform.

May I ask whether DFT is accurate? Because the matrix has all discrete binary elements, and the multiplication process would not tolerate any non-zero error probability in result. However from the things I currently learnt about DFT it seems to be an approximation algorithm?

Also, may I ask roughly how long would the code be? i.e. would this be something I could start from scratch and compose in C++ in perhaps one or two hundred lines? Cause actually this is for a paper...and all I need is that the complexity analyis is $\mathcal{O}(n \log n)$ and the coefficient in front of it doesn't really matter :) So the simplest implementation would be best. (Although I did see some packages like kissfft and FFTW, but they are very lengthy and probably an overkill for my purpose...)

• The example on FFTW's tutorial page is 6 lines long. – Scott Jul 31 '14 at 18:28
• Is the matrix pure ones and zeros? In that case you should probably not go for the DFT approach. If you have fixed-point values you will most likely not get them back exactly through the DFT approach. If you want that exactly there are some alternative related transforms which can keep the integer property. – Oscar Aug 1 '14 at 12:19

When applied correctly, the DFT yields exactly the same result as convolution, it is no approximation (neglecting rounding errors due to limited precision of floating point calculations, of course).

To calculate the circular convolution of vectors $x$ and $y$ (both of length $N$): $$x \circledast y = \operatorname{IDFT}_N\left[\operatorname{DFT}_N[x]\cdot\operatorname{DFT}_N[y]\right]$$ To calculate the linear convolution first append $N-1$ zeros to both vectors yielding the zero-padded vectors $\tilde x$ and $\tilde y$. Then: $$x * y = \operatorname{IDFT}_{2N-1}\left[\operatorname{DFT}_{2N-1}[\tilde x]\cdot\operatorname{DFT}_{2N-1}[\tilde y]\right]$$

To implement the [I]DFT I highly recommend the FFTW library. I promise that it will take you less time to learn its usage than to write your own implementation. For the special case of $x$ and $y$ containing only $1$ and $0$ their might be a more efficient algorithm, though. You would have to look into the Cooley-Tukey algorithm that is used for the FFT and see if you can simplify it for your needs. Whether it gets simpler will also depend on how you define computational complexity - if you define it in terms of additions and multiplications a custom FFT implementation probably won't buy you much. But if you go down to the hardware implementation level you can save some multipliers in the first butterfly stage by replacing them with multiplexers. Anyway, it is well known that the FFT algorithm complexity is of order $N\log N$ as you say. I would like to think that you don't have to show this again.

If you'd like to play around with IFFT and FFT I recommend Octave (an open source Matlab clone).

• why was this answer down voted? @Deve is correct. – JayInNyc Aug 2 '14 at 15:28
• @JayInNyc I'd be interested in that as well. If I made a mistake I'd be glad it was pointed out to me so that I can improve my answer. – Deve Aug 4 '14 at 12:48
• I particularly liked that you distinguished b/w open and circular convolution. Most people don't get that, as important as it is. – JayInNyc Aug 7 '14 at 0:11
• Thank you very much! I have managed to learn FFTW and integrated it into my code. However, there is still one the RAM-efficiency problem...although it works accurately and fast, even the least RAM-intensive float-accuracy would take oevr 20GB of RAM for my dataset of 10^9... – Mike Wong Aug 13 '14 at 20:17
• I was thinking would there be some way to reduce the data node size? I mean currently for input it uses float* (with r2c scheme) and for output it's fftwf_complex which is float*, so this is quite inefficient for only 0 and 1 input...although r2c uses all-real property to halve the usage, still it doesn't make use of binary property...would there be some way to input bool*?? (Or do I have to try and modify the code...since I'm on Windows I'm using the dll libraries, and to modify I'd have to get the gcc code and recompile it all over again...which would be unimaginable right now for me...) – Mike Wong Aug 13 '14 at 20:21