6 ms·

Understanding the FFT Algorithm

114 points by linux_devil 13y ago

af3 13y ago
Python is too slow for FFT. There is no reason to understand how FFT could be implemented in python, if at the end of the day, C or Fortran will be used [1]. [1] http://www.freeware.vasp.de/FFT/fft-pentium.f http://www.freeware.vasp.de/FFT/fft-pentium.f
- sampo 13y ago
  The title of the post if just "Understanding the FFT Algorithm", and his goal is to explain the algorithm. He just uses Python instead of pseudocode, to explain things. The HN title has misled you, in this respect.
- Wilya 13y ago
  Read the article. Quote for you: > "Though the pure-Python functions are probably not useful in practice, I hope they've provided a bit of an intuition into what's going on in the background of FFT-based data analysis." There's very little to be gained by reading pure Fortran code, especially if you don't know the mathematical tricks/optimisations used beforehand. The Python code, on the other hand, stays (relatively) readable. It's not extremely fast (though it stays decent, thanks to Numpy), but it's not the point.
  - af3 13y ago
    > There's very little to be gained by reading pure Fortran code, especially if you don't know the mathematical tricks/optimisations used beforehand. The Python code, on the other hand, stays (relatively) readable. Don't agree. Why trying to understand implementation in the language that will never be used in the production? And, there is a lot to gain by reading Fortran code.
    - Wilya 13y ago
      Because a good part of the optimizations are language-independent. And you need to understand these optimizations before going further. The computations which are redundant and can be done once and stored afterwards don't change depending on whether you're coding in Python of Fortran. The symmetries that you can exploit to compute less stuff don't either. The math behind is always the same. Once you have done that (and gained two or three order of magnitude in computation time), and you're still not fast enough, then okay, it's time to drop to C/Fortran and start playing bookkeeper with your memory allocations. But not before.
      - af3 13y ago
        > Once you have done that (and gained two or three order of magnitude in computation time), and you're still not fast enough, then okay, You are talking about the language (or better implementation) that is compiled to bytecode. The speed increase that you gain by optimization is of the order of magnitude smaller than by just using language that compiles to machine code. Once you understood how it's done in python, then you need to port it to fortran/c, what is the point? Fortran is not much different from python in syntax.
        darkarmani 13y ago
        > Once you understood how it's done in python, then you need to port it to fortran/c, what is the point? At that point, why would you port it to fortran/c instead of using FFTW? The whole point of the exercise was to learn about and optimize the algorithm. If you can learn the easiest using python, why wouldn't you do your exploration in python, since going from python to your language of choice should be the easiest part of the exercise if you choose to do that.
    - dbaupp 13y ago
      > Why trying to understand implementation in the language that will never be used in the production? It's the underlying mathematics that's most interesting: that's how one goes from O(n^2) to O(n log n). And the mathematics is language independent, so it might as well be demonstrated in a reasonably accessible way.
  - sampo 13y ago
    "There's very little to be gained by reading pure Fortran code, especially if you don't know the mathematical tricks/optimisations used beforehand. The Python code, on the other hand, stays (relatively) readable." Depends. His first Python code was: x = np.asarray(x, dtype=float) N = x.shape[0] n = np.arange(N) k = n.reshape((N, 1)) M = np.exp(-2j * np.pi * k * n / N) return np.dot(M, x) I can write more or less the same in Fortran: n = reshape([(i, i=0,len-1)], shape=[1,len]) k = transpose(n) M = exp(-2 * j * pi * matmul(k,n) / len) res = matmul(M, x) But of course Fortran is a typed language, so I need to explicitly write the types for everything, so my whole function gets longer pure function dft_slow(x, len) result(res) complex, intent(in) :: x(len) integer, intent(in) :: len complex :: n(1,len), k(len,1), M(len, len), res(len) integer :: i n = reshape([(i, i=0,len-1)], shape=[1,len]) k = transpose(n) M = exp(-2 * j * pi * matmul(k,n) / len) res = matmul(M, x) end function dft_slow and so yes, maybe less readable. P.S. One must also define somewhere global constants complex, parameter :: j = (0,1) ! imaginary unit real, parameter :: pi = acos(-1.) for the above function to work. Also everything is by default in single precision here.
    - vonmoltke 13y ago
      I have a sudden urge to find an excuse to write Fortran again.
- weland 13y ago
  The usual approach is to use the foreign function interface and work with FFTW or some other high-performance library.
  - af3 13y ago
    that's correct. AFAIK the point of the OP was to understand how that "high-performance library" works.
deleted 13y ago
[deleted]
James_Duval 13y ago
I don't know Python, so I may be missing something, but: I don't understand how the naive approach can be O(n^2) if it involves matrix multiplication. Surely at best it would be somewhere between O(n^2) and O(n^3), and most likely O(n^3)?
- alexchamberlain 13y ago
  There's no matrix multiplication, just vector dot product. (Yes, you can express it as a matrix multiplication, but nurrrr)
  - weland 13y ago
    I don't think I've seen the dot product there. Don't just throw terms around, it's very confusing for beginners :-). While I may have missed your idea about the dot product, I'm quite sure you can't express a dot product as a matrix multiplication; perhaps you were thinking you can express matrix multiplication in terms of dot products, but I don't see how that answer's OP's question.
    - dbaupp 13y ago
      A conventional (i.e. with vectors over ℝ) dot product is just a matrix multiplication between a 1×n and an n×1 matrix.
      - weland 13y ago
        That hardly qualifies as "you can express dot product as matrix multiplication" in the context of the discussion about complexity. I also think it's backwards from a theoretical standpoint. Dot product isn't a cornercase of matrix multiplication, it has a meaning of its own that can be viewed without even introducing matrices. You can view matrix multiplication as computationally equivalent with a sequence of dot products, but that where stuff ends; otherwise, for instance, matrix multiplication also doesn't have all the properties that the dot product has (e.g. it's not commutative). I know a lot of books take this shortcut, but the two concepts are different enough that, IMO, they do not warrant alexchamberlain's remark.
existencebox 13y ago
Coincidentally, I was just reading through the wiki pages on Fourier transforms the other evening, and having a fantastically hard time wrapping my head around some parts of it. Maybe the wrong place to ask, but if anyone had on hand some good intuitive (potentially directed more towards the programmer than the data scientist) description of the Fourier transform?
- wetmore 13y ago
  http://www.altdevblogaday.com/2011/05/17/understanding-the-fourier-transform/ http://www.altdevblogaday.com/2011/05/17/understanding-the-f...
- JonnieCache 13y ago
  Another approach to understanding fourier transforms at the esteemed BetterExplained: http://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/ http://betterexplained.com/articles/an-interactive-guide-to-... This is more about the actual Fourier Transform operation itself, rather than the FFT algorithm. OTOH, it has whizzy interactive animations. The thing that made it simpler to me was the fact that an FT takes us from an oscilloscope to a spectrum analyser and back again.
- sampo 13y ago
  You have a time series of N points. You could use normal curve fitting to fit a constant function + (N-1) sine and cosine waves, and it turns out there is always a perfect fit, whatever your N points are. Discrete fourier transform is "just" a faster way to compute the same result.
- deleted 13y ago
  [deleted]
- KaeseEs 13y ago
  I don't know any data scientists, or even what one is really, but here's the EE explanation: the Fourier transform maps a signal in the time domain to the frequency domain, and vice versa. So if you have a sine wave in the time domain [say something at 3khz: y = sin(3000x)], the Fourier transform would turn that into a pair of Dirac deltas (spikes that are infinitely thin but have an area of one) at omega = +-3000. Or, if you had a pulse, the Fourier transform would map that to a sinc [sinc(x) = sin(x)/x]. The Fourier transform is related to the Laplace transform, which maps signals in the time domain to those in the "S domain", a domain where locations are normally described with complex coordinates which contains information both on frequency of the periodic components of the signal as well as any transients (starting conditions, in layman's terms).
deleted 13y ago
[deleted]
j2kun 13y ago
If anyone is interested in more of the mathematics, I have a few posts on Fourier Analysis (starting with four lengthy posts here: http://jeremykun.com/primers/ http://jeremykun.com/primers/ ) and a follow-up on actually deriving and computing the FFT here: http://jeremykun.com/2012/07/18/the-fast-fourier-transform/ http://jeremykun.com/2012/07/18/the-fast-fourier-transform/ with a (relatively misguided) application of denoising a sound bite. My implementation is in almost pure python (except for complex number arithmetic) and it shows.
- achompas 13y ago
  Your primers look great! Do you have any sources for people interested in learning more about the various topics?
  - j2kun 13y ago
    I think you mean about Fourier analysis, yes? Brad Osgood's Stanford lectures (available on YouTube) give a very good (read: steady, unassuming) approach and the lecture notes associated with that course are extremely detailed and informative if you want a deeper dive.
    - achompas 13y ago
      Yep, definitely Fourier analysis, I haven't studied it formally at all. Thanks for the pointer to Osgood's lectures. A co-worker recommended Stein's "Fourier Analysis" after a bit of work on my part. Do you have any thoughts on that text?
scast 13y ago
The book by Cormen et al. has a fantastic explanation of the FFT too.
- mtdewcmu 13y ago
  This is true. I looked up CLR's treatment recently, and it was surprisingly readable (for CLR). It describes the FT as transforming a polynomial between its coefficient representation and point-value representation (n point values are sufficient to fully represent a polynomial of order-bound n). To yield the point values, the polynomial is evaluated at the n complex roots of unity, because these numbers have useful properties. My background in math is thin when it comes to things like complex numbers, but the explanation more or less made sense to me. CLR was focused on how the FFT was useful for multiplication, so it didn't really go into all the things I would have liked to know.
coherentpony 13y ago
Oh, he didn't do the bit reversal trick! That's the really cool part :)

114 points by linux_devil 13y ago

af3 13y ago

Python is too slow for FFT. There is no reason to understand how FFT could be implemented in python, if at the end of the day, C or Fortran will be used [1]. [1] http://www.freeware.vasp.de/FFT/fft-pentium.f http://www.freeware.vasp.de/FFT/fft-pentium.f

sampo 13y ago
The title of the post if just "Understanding the FFT Algorithm", and his goal is to explain the algorithm. He just uses Python instead of pseudocode, to explain things. The HN title has misled you, in this respect.
Wilya 13y ago
Read the article. Quote for you: > "Though the pure-Python functions are probably not useful in practice, I hope they've provided a bit of an intuition into what's going on in the background of FFT-based data analysis." There's very little to be gained by reading pure Fortran code, especially if you don't know the mathematical tricks/optimisations used beforehand. The Python code, on the other hand, stays (relatively) readable. It's not extremely fast (though it stays decent, thanks to Numpy), but it's not the point.
- af3 13y ago
  > There's very little to be gained by reading pure Fortran code, especially if you don't know the mathematical tricks/optimisations used beforehand. The Python code, on the other hand, stays (relatively) readable. Don't agree. Why trying to understand implementation in the language that will never be used in the production? And, there is a lot to gain by reading Fortran code.
  - Wilya 13y ago
    Because a good part of the optimizations are language-independent. And you need to understand these optimizations before going further. The computations which are redundant and can be done once and stored afterwards don't change depending on whether you're coding in Python of Fortran. The symmetries that you can exploit to compute less stuff don't either. The math behind is always the same. Once you have done that (and gained two or three order of magnitude in computation time), and you're still not fast enough, then okay, it's time to drop to C/Fortran and start playing bookkeeper with your memory allocations. But not before.
    - af3 13y ago
      > Once you have done that (and gained two or three order of magnitude in computation time), and you're still not fast enough, then okay, You are talking about the language (or better implementation) that is compiled to bytecode. The speed increase that you gain by optimization is of the order of magnitude smaller than by just using language that compiles to machine code. Once you understood how it's done in python, then you need to port it to fortran/c, what is the point? Fortran is not much different from python in syntax.
      - darkarmani 13y ago
        > Once you understood how it's done in python, then you need to port it to fortran/c, what is the point? At that point, why would you port it to fortran/c instead of using FFTW? The whole point of the exercise was to learn about and optimize the algorithm. If you can learn the easiest using python, why wouldn't you do your exploration in python, since going from python to your language of choice should be the easiest part of the exercise if you choose to do that.
  - dbaupp 13y ago
    > Why trying to understand implementation in the language that will never be used in the production? It's the underlying mathematics that's most interesting: that's how one goes from O(n^2) to O(n log n). And the mathematics is language independent, so it might as well be demonstrated in a reasonably accessible way.
- sampo 13y ago
  "There's very little to be gained by reading pure Fortran code, especially if you don't know the mathematical tricks/optimisations used beforehand. The Python code, on the other hand, stays (relatively) readable." Depends. His first Python code was: x = np.asarray(x, dtype=float) N = x.shape[0] n = np.arange(N) k = n.reshape((N, 1)) M = np.exp(-2j * np.pi * k * n / N) return np.dot(M, x) I can write more or less the same in Fortran: n = reshape([(i, i=0,len-1)], shape=[1,len]) k = transpose(n) M = exp(-2 * j * pi * matmul(k,n) / len) res = matmul(M, x) But of course Fortran is a typed language, so I need to explicitly write the types for everything, so my whole function gets longer pure function dft_slow(x, len) result(res) complex, intent(in) :: x(len) integer, intent(in) :: len complex :: n(1,len), k(len,1), M(len, len), res(len) integer :: i n = reshape([(i, i=0,len-1)], shape=[1,len]) k = transpose(n) M = exp(-2 * j * pi * matmul(k,n) / len) res = matmul(M, x) end function dft_slow and so yes, maybe less readable. P.S. One must also define somewhere global constants complex, parameter :: j = (0,1) ! imaginary unit real, parameter :: pi = acos(-1.) for the above function to work. Also everything is by default in single precision here.
  - vonmoltke 13y ago
    I have a sudden urge to find an excuse to write Fortran again.
weland 13y ago
The usual approach is to use the foreign function interface and work with FFTW or some other high-performance library.
- af3 13y ago
  that's correct. AFAIK the point of the OP was to understand how that "high-performance library" works.

deleted 13y ago

[deleted]

James_Duval 13y ago

I don't know Python, so I may be missing something, but: I don't understand how the naive approach can be O(n^2) if it involves matrix multiplication. Surely at best it would be somewhere between O(n^2) and O(n^3), and most likely O(n^3)?

alexchamberlain 13y ago
There's no matrix multiplication, just vector dot product. (Yes, you can express it as a matrix multiplication, but nurrrr)
- weland 13y ago
  I don't think I've seen the dot product there. Don't just throw terms around, it's very confusing for beginners :-). While I may have missed your idea about the dot product, I'm quite sure you can't express a dot product as a matrix multiplication; perhaps you were thinking you can express matrix multiplication in terms of dot products, but I don't see how that answer's OP's question.
  - dbaupp 13y ago
    A conventional (i.e. with vectors over ℝ) dot product is just a matrix multiplication between a 1×n and an n×1 matrix.
    - weland 13y ago
      That hardly qualifies as "you can express dot product as matrix multiplication" in the context of the discussion about complexity. I also think it's backwards from a theoretical standpoint. Dot product isn't a cornercase of matrix multiplication, it has a meaning of its own that can be viewed without even introducing matrices. You can view matrix multiplication as computationally equivalent with a sequence of dot products, but that where stuff ends; otherwise, for instance, matrix multiplication also doesn't have all the properties that the dot product has (e.g. it's not commutative). I know a lot of books take this shortcut, but the two concepts are different enough that, IMO, they do not warrant alexchamberlain's remark.

existencebox 13y ago

Coincidentally, I was just reading through the wiki pages on Fourier transforms the other evening, and having a fantastically hard time wrapping my head around some parts of it. Maybe the wrong place to ask, but if anyone had on hand some good intuitive (potentially directed more towards the programmer than the data scientist) description of the Fourier transform?

wetmore 13y ago
http://www.altdevblogaday.com/2011/05/17/understanding-the-fourier-transform/ http://www.altdevblogaday.com/2011/05/17/understanding-the-f...
JonnieCache 13y ago
Another approach to understanding fourier transforms at the esteemed BetterExplained: http://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/ http://betterexplained.com/articles/an-interactive-guide-to-... This is more about the actual Fourier Transform operation itself, rather than the FFT algorithm. OTOH, it has whizzy interactive animations. The thing that made it simpler to me was the fact that an FT takes us from an oscilloscope to a spectrum analyser and back again.
sampo 13y ago
You have a time series of N points. You could use normal curve fitting to fit a constant function + (N-1) sine and cosine waves, and it turns out there is always a perfect fit, whatever your N points are. Discrete fourier transform is "just" a faster way to compute the same result.
deleted 13y ago
[deleted]
KaeseEs 13y ago
I don't know any data scientists, or even what one is really, but here's the EE explanation: the Fourier transform maps a signal in the time domain to the frequency domain, and vice versa. So if you have a sine wave in the time domain [say something at 3khz: y = sin(3000x)], the Fourier transform would turn that into a pair of Dirac deltas (spikes that are infinitely thin but have an area of one) at omega = +-3000. Or, if you had a pulse, the Fourier transform would map that to a sinc [sinc(x) = sin(x)/x]. The Fourier transform is related to the Laplace transform, which maps signals in the time domain to those in the "S domain", a domain where locations are normally described with complex coordinates which contains information both on frequency of the periodic components of the signal as well as any transients (starting conditions, in layman's terms).

deleted 13y ago

[deleted]

j2kun 13y ago

If anyone is interested in more of the mathematics, I have a few posts on Fourier Analysis (starting with four lengthy posts here: http://jeremykun.com/primers/ http://jeremykun.com/primers/ ) and a follow-up on actually deriving and computing the FFT here: http://jeremykun.com/2012/07/18/the-fast-fourier-transform/ http://jeremykun.com/2012/07/18/the-fast-fourier-transform/ with a (relatively misguided) application of denoising a sound bite. My implementation is in almost pure python (except for complex number arithmetic) and it shows.

achompas 13y ago
Your primers look great! Do you have any sources for people interested in learning more about the various topics?
- j2kun 13y ago
  I think you mean about Fourier analysis, yes? Brad Osgood's Stanford lectures (available on YouTube) give a very good (read: steady, unassuming) approach and the lecture notes associated with that course are extremely detailed and informative if you want a deeper dive.
  - achompas 13y ago
    Yep, definitely Fourier analysis, I haven't studied it formally at all. Thanks for the pointer to Osgood's lectures. A co-worker recommended Stein's "Fourier Analysis" after a bit of work on my part. Do you have any thoughts on that text?

scast 13y ago

The book by Cormen et al. has a fantastic explanation of the FFT too.

mtdewcmu 13y ago
This is true. I looked up CLR's treatment recently, and it was surprisingly readable (for CLR). It describes the FT as transforming a polynomial between its coefficient representation and point-value representation (n point values are sufficient to fully represent a polynomial of order-bound n). To yield the point values, the polynomial is evaluated at the n complex roots of unity, because these numbers have useful properties. My background in math is thin when it comes to things like complex numbers, but the explanation more or less made sense to me. CLR was focused on how the FFT was useful for multiplication, so it didn't really go into all the things I would have liked to know.

coherentpony 13y ago

Oh, he didn't do the bit reversal trick! That's the really cool part :)