Interesting debate. I don't disagree with you but I have to make a few comments:)
It could be argued that numbers, 1, 2, 3 etc don't exist. They are just a mathematical short hand, a mental short cut, a figment of our imaginations. "3" is just short hand for "the successor of 2", "2" is just the short hand for "the successor of "1", "1" is just to represent "a thing" (whatever it may be) as opposed to "0" which is "not to have the thing".
We learn about numbers at an early age and live with them for such a long time that they take on a reality like the nose in front of your face. As my young son said when I was counting ducks in a picture book with him "No dad, that's not three ducks, that's just another duck". I could not fault his logic.
Clearly there are no numbers inside a computer. Just ethereal patterns of transistors switched on or off that we interpret as numbers. Then there are the networks of other transistors that shuffle those patterns around in ways that we like to imagine are operations like addition, subtraction etc.
Still, this view of things rarely gets us to any practical result.
If you have a number of numbers, representing a signal, you can see, there is already something confusing: number of numbers
I hope that "number of numbers" is not too confusing or it's going to make handling arrays a bit tricky. Now, a "complex number of numbers" gets interesting:)
Fourier works on a computer, but mostly you do NOT get what you want. You get what you want, if there is a certain problem from which you know, it can be solved using FT, but this is not usual. Only if the bunch of problems you want to solve consists of problems, where FT is a way to solve, you will be successful.
This true BUT..I could say "Addition works on a computer, but mostly you do NOT get what you want". Point being that addition is defined mathematically to work in ways that a computer mostly cannot handle. Even just thinking about integers, a mathematical "+" is supposed to work for all integers, positive and negative, form zero to plus/minus infinity. A computer being a finite machine can only handle a vanishingly small part of the mathematical concept of addition. That's before we get onto addition of real numbers, vectors, matrices etc.
Of course you have to be aware of your machines limits, it's number range, it's memory capacity etc then you can solve real problems in some limited domain.
But, what you don't say and what most people are not aware of: There is one additional condition: If you have n data points, there is a data point n+1 equal to data point 0, n+2 = 2 and so on, over all the times. Only then a FT makes sense.
I disagree, if I have n data points I have n data points, that it. For sure Bean's Fourier transform algorithm knows nothing more than those n points and it seems to do quite well with out any deep thoughts about infinitely repeating sequences.
If this would not be the case, you were never allowed to talk about "frequencies".
This is clearly not true. If I have n samples and they go like -1, +1, -1, +1, -1.... I can easily see n/2 alternations or cycles of the pattern -1, +1 in my sample set. Given that my sample set represents a length of time I can say I have a frequency.
And in this way it is possible to make sense of the Fourier transform without bringing in the concept of infinity.
The point is: there are no frequencies in between those harmonics.
Yep, think I said that earlier when discussing the 1.0Hz and 1.3Hz issue.
It is just this: lets say, the function is a sinusoid of one period in the interval. Than the fourier transform gives a spectrum 1, 0, 0.... 0, 0 , 1.
Actually this not quite right. Let's say you have a single cycle of a sine wave, amplitude 1.0, represented in 16 samples and you perform the discrete Fourier transform as per the usual definition LINK HERE then the result will be 0.0 in the first output bin (The DC level), 0.5 in the second frequency bin (1 cycle/sampling time) and then 0.5 in the last frequency bin.
Why is that? Well during the transform we are performing a multiply of the input sine by a test frequency sine. When the test frequency is the same as our input we have sin * sin or sin squared. The average of which is 0.5.
So half the energy of our input sine wave has shown up at a frequency of 1 and half the energy has turned up at a frequency of -1 !!!
Anyway, I'm just quibbling here:)
We can not decide, if the result comes from a real waveform or from a wrong selection of the interval.
I agree. And I think it is a very important point. If anyone is going to use a Fourier transform as performed by a computer they must be aware of it's limitations. Otherwise it can lead to all sorts strange looking results and confusion.
OK, this really is a interesting discussion. Arguments going forward and backward, but today I am running out of time. I will try to set up another cmap and try do discuss point for point.
Very interesting, because we sometimes are convergent and sometimes divergent and sometimes we don't know, what we are.
2) I have just had an insight into the workings of the FAST Fourier transform.
A few decades back I saw my first listing of a program that performed a Fast Fourier Transform. It was in BASIC for the Atari ST. No explanation was given for how in worked. "No problem" I thought as I had understood something of the Fourier transform back in university. How wrong I was.
On reading the listing I found it started off with a weird thing. One by one it took all the elements of the input array, and moved them into an output array, such that the elements index in the output array was formed by taking its index in the input array and reversing all the bits.
So element at %00000000 was written to %00000000, element at %00000001 was written to %10000000, element at &00000010 was written to %01000000. Etc etc.
The rest of the program then operated on that reordered array of samples.
Needless to say I could not fathom why one would want to do such a crazy nuts thing and the rest of the code did not help.
However today I awoke from my afternoon nap with a nagging thought in my mind, a realization, another little step forward in the journey to understand the FFT. Given that I should have seen this decades ago I am now convinced of my dummy status.
Part 1 - Bit Reversal.
Apologies to all FFT gurus for whom this is all very dull. Here we go:
1) Imagine we want to add up 8 numbers S0, S1, S2, S3 .... S7.
2) Say we had some function to do this. Call it SUM(). Then we could write:
result = SUM(S0, S1, S2, S3, S4, S5, S6, S7)
4) But just to be weird we could use SUM on all the even position numbers separately, then all the odd position numbers, then add the two results. So we write:
This would obviously give the same overall summation.
5) To be more weird we could do that again for each of our two SUM(...)s. Taking first the odd positioned S's and then the even positioned S's. We write:
Well, blow me down, the binary representation of the final order is the mirror image, bit reversal, of the starting order in binary! And so it will be even if we have bigger sets of numbers.
And that is why the famous FFT bit reversal or "butterfly" is so, well, famous. Slaps head very hard.
Part 2 - Why is the important?
If you look at Bean's Fourier Spin code example you will see that it has two nested repeat loops each performing as many repetitions as there are input samples. For each of N output frequencies there is is N multiplications of input samples by sin and cos. That means the execution time goes up by N squared and rapidly becomes very slow for large sample sets.
Let's call the multiplications by sin and cos the "guts" of the algorithm.
For 8 samples there are 8 times 8 "guts" to perform = 64 "guts"
For 1024 samples there are 1024 * 1024 "guts" to perform = 1048576 "guts"
This is not practical. That's a lot of multiplications.
Turns out that for the Fourier transform you can split the job up into halves. Take the transform of each half and then add the results. Using every even sample for the first half and every odd sample for the second half. Much as we did for the SUM process above.
Again, as we did for sum above, you can split the halves into two, take the transform of those and add them.
This can go on recursively splitting the samples in half, transforming them and adding the results.
Why is that important?
Well if you start with 16 samples you have 16 * 16 = 256 "guts" to perform
If you can spilt that into two sets of 8 samples you only have (8 * 8) + (8 * 8) = 128 "guts"
If you can split that into 4 sets of 4 samples it's (4*4) + (4*4) + (4*4) + (4*4) = 64 "guts"
And so on. By splitting the thing up we get rid of the "N squared" execution time problem.
Given that the transform is split by "odds" and "evens" at each level, as in my SUM game above, clearly the same "bit reversal" re-odering of the samples is required.
As you see I'm still miles away from "getting" the FFT or being able to code one. But a little light has come on:)
I'm still not clear on how the job is split into halves or exactly what you have to do to handle the halves correctly.
P.S. Does anyone else have the same problems with this that I do?
@Heater, Thanks that is a great explaination of the "Bit Reversal" process.
What has me stumped is these "Bufferfly" operations. I have no idea what they are. Can you shed some light on them as well ?
"
In the context of fast Fourier transform algorithms, a butterfly is a portion of the computation that combines the results of smaller discrete Fourier transforms (DFTs) into a larger DFT, or vice versa (breaking a larger DFT up into subtransforms). The name "butterfly" comes from the shape of the data-flow diagram in the radix-2 case, as described below.[1] The same structure can also be found in the Viterbi algorithm, used for finding the most likely sequence of hidden states.
Most commonly, the term "butterfly" appears in the context of the Cooley–Tukey FFT algorithm, which recursively breaks down a DFT of composite size n = rm into r smaller transforms of size m where r is the "radix" of the transform. These smaller DFTs are then combined with size-r butterflies, which themselves are DFTs of size r (performed m times on corresponding outputs of the sub-transforms) pre-multiplied by roots of unity (known as twiddle factors). (This is the "decimation in time" case; one can also perform the steps in reverse, known as "decimation in frequency", where the butterflies come first and are post-multiplied by twiddle factors. See also the Cooley–Tukey FFT article.)"
Thanks Heater, this is surely the most simple way to understand why bit reversal works. And this can show, why symmetry is such an important concept. You could say "bit reversal" and think about how to implement this with an algorithm on a computer. But you also could place a mirror right to the index 0, 1, 2, 3, 5 and you will see in the mirror the new index. But only if you represented the original index as a binary number with two symbols. Now think, that you have a next generation optical or photonic propeller (oppro or phopro) and you have leds or quantum lasers, showing the original index in binary, you mirror this light to photo transistors and this give you the index reordered in real-time. This will be a completely different implementation of the same principle. But, isn't that, what a actual optical instrument like a microscope or the eye of an insect does?
Never stop thinking. Even if you are told to.
What has me stumped is these "Butterfly" operations.... Can you shed some light on them as well ?
Earlier I said something about the "blind leading the blind" but I'm starting to get a glimmer of light re: FFT.
Let' play with that SUM function I described earlier. SUM(S0, S1, S2, S3, S3, S5, S6, S7)
Let's imagine that we actually implemented that SUM function in some language and that our implementation was recursive.
In pseudo code it might look like:
REAL SUM (LIST s)
BEGIN
IF LENGTH(S) = 1
RETURN (FIRST(S)) ; Only one element, nothing to sum, so just return it.
ELSE
r1 = SUM (EXTRACT_EVEN_ELEMENTS(s))
r2 = SUM (EXTRACT_ODD_ELEMENTS(s))
RETURN (r1 + r2)
ENDIF
END
Here we assume a language that can handle lists of our S type things and we have functions to extract the odd and even position elements from lists.
We might call this function like:
s = LIST(S0, S1, S2, S3, S4, S5, S6, S7)
result = SUM(s)
Now, if we think about how that call proceeds, we see that a SUM over 8 elements calls SUM over 4 elements which calls SUM over 2 elements, which calls SUM over 1 element which just returns the element.
Graphically we can see the calls proceeding as in the attached picture where I have put the lists involved in each call to SUM into boxes and added lines to track the progress of individual elements.
Now:
1) If you look at that picture you will see why the term "butterfly" is used. The way the paths of elements criss-cross starts to look like butterflies.
2) It seems that doing this recursive process performs the exact same swapping around of elements as the bit reversal algorithm.
Turns out that the fourier transform can be written recursively much like our SUM function. Except that each call of the function actually returns a list of values. Each of those values is in two parts, a sin part and a cosin part. Also it is required to multiply one of the "sub DFT" by a "fiddle factor".
What this means is that if one were to implement a recursive DFT then no bitreversal step is required. The ways the calls progress just does that for you:)
The bit reversal step is required when implementing the DFT non-recursively as is normally done.
You will see it looks much like my SUM function above.
You will also see that the part many people casually refer to as the "fiddle factor" is actually the guts of the thing where the multiplying by sin and cos happens. After all the rest of the function does no calculating at all.
I'm going to have a go at implementing that recursive FT pseudo code as directly as possible in C++ with list objects as it's the closest thing to a Fast Fourier Transform I almost understand at the moment.
Here is an interesting piece of trivia. The Fast Fourier Transform as we know it is attributed to J.W. Cooley and John Tukey as they developed it for computers in the 1960's. Back then Cooley said:
someday radio tuners will operate with digital processing units. I have heard this suggested with tongue in cheek, but one can speculate.
Trying to gain a deeper understanding, I am still far from having it. Now, what happens to the original signal when doing this bit-shift?
Suppose, we have a sinusoid with the fundamental frequency. Taking all the even samples first and the odd later, we compress a full oscillation to the first and to the second half, so we double the frequency of the signal! Doing this again with the two halves, we quadruplicate frequency, so if we do this often enough, he have the highest possible frequency!
What happens, when the origin was the highest frequency, so every period only is samples at + and - amplitude, then sorting odd even will bring 1 to the first half of the interval, -1 to the second. Further steps will change nothing, as all values in the first half are equal.
This remark is just to have a picture in the windings of the brain, what's going on. I see Chip in his video on speech generation, how he demonstrates the phase shift. Maybe, he already knows, how to explain FFT.
Since you asked about the Aha moment...
By combining Beans document which to me explains what we are trying to do combined with heaters description of the butterfly which to me explains how to do the FFT in a practical sense combined with the little I know and can remember, I now have had that Aha insight.
I could never really understand the FFT even though I had all the math 40 years ago and have actually programmed the FFT (well OK, actually copied and used) in several projects over the years. I feel now at least I have a chance at understanding, so well done to all but don't stop because there are some tantalizing hints about more to come.
That approach splits the DFT and calls itself with the two halves recursively. Recursion stops when the is only one element in the set. This means it needs no bit revesal stage that the optimized, "in place" implementations use. As I pointed out with my SUM() game, the way the sample lists are split and passed around does all the reordering anyway.
Suppose, we have a sinusoid with the fundamental frequency..... we have the highest possible frequency!
I can't answer that but already I can see we are looking at things the opposite way around again:)
At the top level a DFT of 8 samples can have a max frequency of 4.
Split to two DFTs of 4 samples each can have a max frequency 2
Spilt to four DFT's of 2 each can have a max frequency of 1
Spilit to 8 DFT's of a single sample each can only have a max frequency of 0 or DC !
Which is exactly what the pdeudo code on wiki does, for a 1 sample DFT it returns a single value, a DC level.
Glad that butterfly thing helped a bit, I was afraid it would just confuse more.
I could never really understand the FFT even though I had all the math 40 years ago and have actually programmed the FFT (well OK, actually copied and used) in several projects over the years.
Exactly, same here. So near but so far, very frustrating.
@heater
Well as you pointed out, Beans sample program has multiple loops inside loops resulting in lots of operations in the inner loop. And then this little butterfly comes along and magically reduces the inner loop count not by a little but by a significant amount. And Viola the FFT on a slow computer suddenly becomes very doable.
I am amazed that these two young graduates had the insight to see this reordering that would make the FFT practical at a time when "computers" really were slow and microprocessors were not even a sparkle in Chips eye. Well, we really owe them a debt because all the noise our games make and audio and signal processing (and much more) is based on this algorithm.
Maybe you should wear some head gear (to protect the wall of course).
Now, what happens to the original signal when doing this bit-shift?
I think I have an answer to your question about what happens when the samples are split into ever smaller halves.
As you say it seems that if you take all the odd samples of a sine wave and put them into a half of the list then you appear to have doubled the frequency.
BUT. I'm just looking at that recursive FFT pseudo code on wiki (reproduced below). Here we can see that at each recursive call level the list is chopped up into odds and evens which end up in the first and second halves of a new list. List X goes in, list Y comes out.
Look what happens to the list length N. At each recursive call level it is divided by two.
That N is used in the "exp(...)" gobbledygook, the so called "twiddle factors"
Well, apologies to the math shy, the twiddle factor is actually our sin and cos work where we multiply samples by our "test" sin/cos waves because:
exp(−2πik/N) = cos(2πk/N) - i*sin(2πk/N)
See that dividing N by 2 actually doubles the frequency of the test sin and cos.
Conclusion: Whilst splitting the samples into ever smaller halves we are also doubling the frequency of the sin and cos waveforms we are using as tests for the presence of frequencies. So it all works out:)
Is this a good bit of hand waving or not? Or even correct logic?
Y0,...,N−1 ← ditfft2(X, N, s): DFT of (X0, Xs, X2s, ..., X(N-1)s):
if N = 1 then
Y0 ← X0 trivial size-1 DFT base case
else
Y0,...,N/2−1 ← ditfft2(X, N/2, 2s) DFT of (X0, X2s, X4s, ...)
YN/2,...,N−1 ← ditfft2(X+s, N/2, 2s) DFT of (Xs, Xs+2s, Xs+4s, ...)
for k = 0 to N/2−1 combine DFTs of two halves into full DFT:
t ← Yk
Yk ← t + exp(−2πi k/N) Yk+N/2
Yk+N/2 ← t − exp(−2πi k/N) Yk+N/2
endfor
endif
Here, ditfft2(X,N,1), computes Y=DFT(X) out-of-place by a radix-2 DIT FFT, where N is an integer power of 2 and s=1 is the stride of the input X array. X+s denotes the array starting with Xs.
Notice how the pseudo code for the recursive FFT looks remarkably similar to the recursive SUM() pseudo code I posted above. It divides into odds and even position elements in the same way, Thus stirring up the sample order in the same way as the bit-reversal stage used in non-recursive (in place) implementations.
Of course the FFT calls return a new list rather than a single sum value. Then each element of one half of the list is multiplied by the "exp gobbledygook" (Our sin and cos test frequencies) prior to adding it to the output.
This is really interesting, I hope I find some time to get involved. I've done a lot with Fourier in optics (because a lens is a Fourier transformer and light can be represented with an angular spectrum etc etc) and I'm pretty happy with it however the F in the FFT always did confuse me. I'm another guy who just ripped it out of numerical recipes in C or worse, just uses the FFT command in Matlab.
One aspect to consider of the final form of this discussion is the Path of Discovery. The target audience is those that "need to use" FFT/DFT but do not necessarily have a deep math background. But, we notice the "A-ha!" moments tend to come after "side trips" to more complicated references.
The ultimate "simple" explanation might benefit from including the links to these more complicated side references, in a similar order as found in the thread.
I'm guessing new explorers will have to read the final explanation several times, going a bit deeper each pass. So far I've read through several times; its still watering dry ground, but getting better.
This is really good work, I appreciate it and I'm sure other do and more will.
The "F" in "FFT" is starting to take on a whole other meaning for me:)
OK. I think we have uncovered a dirty little secret here.
What we have is a world full of software engineers, electronics engineers, etc etc, loaded up with appropriate degrees, and diplomas in maths, physics, electronics, computer science etc. Many of them use the Fast Fourier Transform from time to time in all kind of contexts without ever knowing how it works or how to make one. And probably without wanting to admit it, perhaps thinking that everyone else "gets it" and they don't want to appear dumb. The secret is, everyone else does not get it either:)
This is kind of odd given the ubiquity of the FFT.
Sure we all use micro-processors without knowing how to design one, and we may use other maths packages that we could not create. But I think this is a little different. The maths is all out there in the open, source codes implementing FFT can be found in many places. The original Cooley-Tukey paper is one of the most sited ever. The thing is not actually that big.
Let's do an informal survey: Stand up everyone here who fits my description of the typical engineer above. (Heater stands up). Hmm... that's a lot.
The way I see it is that at one end there is the impenetrable maths with it's integrals, exponentials, imaginary numbers, Euler identities, complex roots of unity etc ect. It is possible to understand the FT from that and get through the typical questions they might ask of a physics/electronics undergraduate. But at the same time not have a real feel for what is actually going on and for sure being a long way from the "F" in "FFT". That's me.
On the other hand there is the source code you can find everywhere. This is totally unintelligible unless you know all the little tricks that went into creating it and so does not help in understanding at all.
Then there are all the explanations of how to put the "F" into "FFT" and get from the maths to the code. Most of those seem to be written by people who either don't really want to tell you, haven't really understood the thing just regurgitated what they saw elsewhere, or do really know what they are doing but are so familiar with it all that they can't imagine the state of mind of their dumbfounded audience.
Is understanding the FFT like some kind of mental illness? Those who don't have it can't possibly imagine what it is like to be a sufferer. Those who do have it can't possibly explain what it's like:)
Well, yes, mathematics has often been described as "self analysis", you can't just learn it like a bunch of historical facts. No one can explain it to you, although they might give hints and clues. You have to see it yourself and convince yourself of the truth of things.
OK, I'm putting by head protection back on....there must be a way to explain the "F" in "FFT" this without Euler.
One aspect to consider of the final form of this discussion is the Path of Discovery.
Yep, essential, that's what I'm interested in. Especially as I haven't trodden much of the path yet.
The target audience is those that "need to use" FFT/DFT but do not necessarily have a deep math background.
Yep, but not just use, it's those that have a nagging itch to know what is going on.
But, we notice the "A-ha!" moments tend to come after "side trips" to more complicated references.
Side trips, yes. I'm hoping those side trips can be simple thought experiments based on common knowledge that lead step by step to an A-ha moment. Like my attempt at showing how bit reversal comes about naturally from the SUM() game. That was part of the magic of Beans CORDIC paper, start with some very familiar ideas, explore them a bit, then put them into the final picture.
So far I've read through several times; its still watering dry ground, but getting better.
If you have any sticking points in Beans document you should speak up. Perhaps Bean could smooth the way in a later version. Often one forgets where one got stuck when one gets to the end of the journey.
OK heater, I'm standing with you and agree with your assessment of the "dirty secret" but the good news is we are getting close.
I think at this point it is time for a first draft, a target if you will, that is going to take a few broadsides maybe even get sunk.
I would suggest this obvious layout:
A statement of what we are doing and why
Then hit em with Beans graphs and description.
Include some pseudo code (and perhaps some spin in an appendix to pay homage to the prop)
Next we state explicitly the problem ie the magnitude of the computation in a practical sense.
Now we hit em with Heaters description of the butterfly
put the two together and expand the guts (pseudo code?)
Now introduce recursion and show the connection to the butterfly
a summary complete with pseudo code and Spin and if someone is crazy enough Assembly
Ok guys let the shooting begin.
Ray
Just read through the thread and I now feel we need to expand the lead in to include some digital sampling and introduce the digital sampling concept along with enough trig to show how to express the sample points on a graph. this can be further expanded when we get to the "guts" part.
I've done it, reached enlightenment, tasted nirvana, come to the end of a long road. My guru status is assured (In my mind anyway). I am not a dummy, I am a free man. The world is mine. Is that Dr Evil you hear laughing? No, it's me:)
Well, OK, what actually happened is that I managed to:
a) Understand enough of the maths of the DFT to convince myself I know how it works. No, feel how it works.
b) Understand enough of the Cooley-Tukey idea to know why splitting the DFT up into halves repeatedly works, and how to do it.
c) Understand how the pseudo code on wikipedia for the recursive splitting works and that it does actually do what is said of it.
d) Now the biggy, actually implement that pseudo code as a real program in C++.
And it works, it does just as Cooley and Tukey and Gauss said it would! The chain is complete.
Calming down a bit, this may be the end of a long road, it's taken me over twenty years to get this far. But it may also be the start of another long road. However I'm savouring the view from here with a large can of beer.
Ahead is:
1) How to rewrite this in C with out using the short hand provided by the "complex" class of C++.
2) This is recursive and requires separate input and output buffers. How to get rid of the recursion and do all the work "in-place" to save memory.
3) Creating a Spin FFT and then perhaps a PASM FFT.
And worst of all, how to create a paper describing how this all works, I still don't see any escape from the heavy maths for the un-mathematical audience. Which is what this is all about. Perhaps "Fast Fourier Transform For Advanced Dummies" where all the hairy maths is illuminated in minute detail but keeping the scope as narrow as possible
Attached is my effort and a plot showing firstly 1024 samples of a sine wave followed by 512 frequency bins containing a peek at frequency = 16.
By the way someone has already posted an FFT in PASM here a while back, forget who it was, just wondering why he has not piped up on this thread.
All right since I shoved my foot in the door ready or not here is a text document which I think we can use as the launch point. There is much still to be done - have fun
Heater,
I just saw your post and while I'm not quite at your place yet I can see your campfire and hear the fireworks.
The text document I just posted does not pretend to get too deep into the math . I think my target audience is an intelligent engineer, like all of us here, who because of need or pride just has to go a few steps deeper without climbing all the stairs in the ivory tower. My intention is kinda like I treat my car. I like knowing how it works but grease on my hands is icky and not for me.
I made a mistake re: the definition of "butterfly". I have to offer my apologies for any confusion and try to make amends.
Bean posted a question about the "butterfly". I posted an explanation, a revelation for me at the time, of how the magical bit-reversal stage works in the common FFT implementations. I used the SUM() game as an example and referred to it as "the butterfly".
This is not common meaning of "butterfly" in FFT circles.
Turns out that whilst the samples are shuffled around by the bit-reversal stage of a common FFT, producing "butterfly" like patterns, as the actual work progresses the data is being constantly re-shuffled such that it is in the correct order when done. It is this re-shuffling, during the work stage, that is commonly referred to as "the butterfly".
What can I say, sorry, I was was misled in that both shufflings are called butterfly but the later use is far more common in everything I have read since. For example see the flow chart in chapter 12 of the DSP Guide: http://www.dspguide.com/ch12/2.htm
Well it looks like I tripped on something on the way to your campfire.
There will be a slight delay while this dummy reattaches the missing parts.
I had not realized that the decomposition had to be undone so this is another deeper aha moment for me. I suppose that if I had gotten some grease on my hands by actually trying to implement the algorithm I would have realized I was missing something. But that is what this process is all about.
so yippie - i think.
Post edit - time passes -
So as a result of my misunderstanding I have realized two important things.
first I realized that when we finish the reordering and have the scrambled array that we have shifted from the time domain to the frequency domain. Recall from Beans diagrams that if we have an Nhz waveform and multiply it by an Nhz waveform the result is N or in other words nothing has to be done to treat it as in the frequency domain.
The second thing I realize with your help of course, is that the reverse reordering also takes place in the frequency domain such that the final array is back in the correct sequence and is the frequency spectrum, and that the butterfly (where the guts live) takes place in the frequency domain.
Your reference to chapter 12 of the Scientists and engineers guide is what did the trick for me.
My head hurts but it feels good. I'll modify the first draft text shortly...
Comments
Interesting debate. I don't disagree with you but I have to make a few comments:)
It could be argued that numbers, 1, 2, 3 etc don't exist. They are just a mathematical short hand, a mental short cut, a figment of our imaginations. "3" is just short hand for "the successor of 2", "2" is just the short hand for "the successor of "1", "1" is just to represent "a thing" (whatever it may be) as opposed to "0" which is "not to have the thing".
We learn about numbers at an early age and live with them for such a long time that they take on a reality like the nose in front of your face. As my young son said when I was counting ducks in a picture book with him "No dad, that's not three ducks, that's just another duck". I could not fault his logic.
Clearly there are no numbers inside a computer. Just ethereal patterns of transistors switched on or off that we interpret as numbers. Then there are the networks of other transistors that shuffle those patterns around in ways that we like to imagine are operations like addition, subtraction etc.
Still, this view of things rarely gets us to any practical result.
I hope that "number of numbers" is not too confusing or it's going to make handling arrays a bit tricky. Now, a "complex number of numbers" gets interesting:)
This true BUT..I could say "Addition works on a computer, but mostly you do NOT get what you want". Point being that addition is defined mathematically to work in ways that a computer mostly cannot handle. Even just thinking about integers, a mathematical "+" is supposed to work for all integers, positive and negative, form zero to plus/minus infinity. A computer being a finite machine can only handle a vanishingly small part of the mathematical concept of addition. That's before we get onto addition of real numbers, vectors, matrices etc.
Of course you have to be aware of your machines limits, it's number range, it's memory capacity etc then you can solve real problems in some limited domain.
I disagree, if I have n data points I have n data points, that it. For sure Bean's Fourier transform algorithm knows nothing more than those n points and it seems to do quite well with out any deep thoughts about infinitely repeating sequences.
This is clearly not true. If I have n samples and they go like -1, +1, -1, +1, -1.... I can easily see n/2 alternations or cycles of the pattern -1, +1 in my sample set. Given that my sample set represents a length of time I can say I have a frequency.
And in this way it is possible to make sense of the Fourier transform without bringing in the concept of infinity.
Yep, think I said that earlier when discussing the 1.0Hz and 1.3Hz issue.
Actually this not quite right. Let's say you have a single cycle of a sine wave, amplitude 1.0, represented in 16 samples and you perform the discrete Fourier transform as per the usual definition LINK HERE then the result will be 0.0 in the first output bin (The DC level), 0.5 in the second frequency bin (1 cycle/sampling time) and then 0.5 in the last frequency bin.
Why is that? Well during the transform we are performing a multiply of the input sine by a test frequency sine. When the test frequency is the same as our input we have sin * sin or sin squared. The average of which is 0.5.
So half the energy of our input sine wave has shown up at a frequency of 1 and half the energy has turned up at a frequency of -1 !!!
Anyway, I'm just quibbling here:)
I agree. And I think it is a very important point. If anyone is going to use a Fourier transform as performed by a computer they must be aware of it's limitations. Otherwise it can lead to all sorts strange looking results and confusion.
Very interesting, because we sometimes are convergent and sometimes divergent and sometimes we don't know, what we are.
1) I have today decided that I am a dummy.
2) I have just had an insight into the workings of the FAST Fourier transform.
A few decades back I saw my first listing of a program that performed a Fast Fourier Transform. It was in BASIC for the Atari ST. No explanation was given for how in worked. "No problem" I thought as I had understood something of the Fourier transform back in university. How wrong I was.
On reading the listing I found it started off with a weird thing. One by one it took all the elements of the input array, and moved them into an output array, such that the elements index in the output array was formed by taking its index in the input array and reversing all the bits.
So element at %00000000 was written to %00000000, element at %00000001 was written to %10000000, element at &00000010 was written to %01000000. Etc etc.
The rest of the program then operated on that reordered array of samples.
Needless to say I could not fathom why one would want to do such a crazy nuts thing and the rest of the code did not help.
However today I awoke from my afternoon nap with a nagging thought in my mind, a realization, another little step forward in the journey to understand the FFT. Given that I should have seen this decades ago I am now convinced of my dummy status.
Part 1 - Bit Reversal.
Apologies to all FFT gurus for whom this is all very dull. Here we go:
1) Imagine we want to add up 8 numbers S0, S1, S2, S3 .... S7.
2) Say we had some function to do this. Call it SUM(). Then we could write:
4) But just to be weird we could use SUM on all the even position numbers separately, then all the odd position numbers, then add the two results. So we write:
This would obviously give the same overall summation.
5) To be more weird we could do that again for each of our two SUM(...)s. Taking first the odd positioned S's and then the even positioned S's. We write:
6) And weirder still we could split it up again:
7) But now our SUM function is not doing any adding at all, it just returns its argument, so we could write:
8) But, coming back to our senses, we realize we have a function to do that, it's called SUM, so we write: Which of course is the same as where we started in 2) but with the parameters in a funny order.
Now, let's look at the order of the arguments we started with and the order we ended up with in both decimal and binary:
Notice anything odd?
Well, blow me down, the binary representation of the final order is the mirror image, bit reversal, of the starting order in binary! And so it will be even if we have bigger sets of numbers.
And that is why the famous FFT bit reversal or "butterfly" is so, well, famous. Slaps head very hard.
Part 2 - Why is the important?
If you look at Bean's Fourier Spin code example you will see that it has two nested repeat loops each performing as many repetitions as there are input samples. For each of N output frequencies there is is N multiplications of input samples by sin and cos. That means the execution time goes up by N squared and rapidly becomes very slow for large sample sets.
Let's call the multiplications by sin and cos the "guts" of the algorithm.
For 8 samples there are 8 times 8 "guts" to perform = 64 "guts"
For 1024 samples there are 1024 * 1024 "guts" to perform = 1048576 "guts"
This is not practical. That's a lot of multiplications.
Turns out that for the Fourier transform you can split the job up into halves. Take the transform of each half and then add the results. Using every even sample for the first half and every odd sample for the second half. Much as we did for the SUM process above.
Again, as we did for sum above, you can split the halves into two, take the transform of those and add them.
This can go on recursively splitting the samples in half, transforming them and adding the results.
Why is that important?
Well if you start with 16 samples you have 16 * 16 = 256 "guts" to perform
If you can spilt that into two sets of 8 samples you only have (8 * 8) + (8 * 8) = 128 "guts"
If you can split that into 4 sets of 4 samples it's (4*4) + (4*4) + (4*4) + (4*4) = 64 "guts"
And so on. By splitting the thing up we get rid of the "N squared" execution time problem.
Given that the transform is split by "odds" and "evens" at each level, as in my SUM game above, clearly the same "bit reversal" re-odering of the samples is required.
As you see I'm still miles away from "getting" the FFT or being able to code one. But a little light has come on:)
I'm still not clear on how the job is split into halves or exactly what you have to do to handle the halves correctly.
P.S. Does anyone else have the same problems with this that I do?
What has me stumped is these "Bufferfly" operations. I have no idea what they are. Can you shed some light on them as well ?
Bean
http://en.wikipedia.org/wiki/Butterfly_diagram
"
In the context of fast Fourier transform algorithms, a butterfly is a portion of the computation that combines the results of smaller discrete Fourier transforms (DFTs) into a larger DFT, or vice versa (breaking a larger DFT up into subtransforms). The name "butterfly" comes from the shape of the data-flow diagram in the radix-2 case, as described below.[1] The same structure can also be found in the Viterbi algorithm, used for finding the most likely sequence of hidden states.
Most commonly, the term "butterfly" appears in the context of the Cooley–Tukey FFT algorithm, which recursively breaks down a DFT of composite size n = rm into r smaller transforms of size m where r is the "radix" of the transform. These smaller DFTs are then combined with size-r butterflies, which themselves are DFTs of size r (performed m times on corresponding outputs of the sub-transforms) pre-multiplied by roots of unity (known as twiddle factors). (This is the "decimation in time" case; one can also perform the steps in reverse, known as "decimation in frequency", where the butterflies come first and are post-multiplied by twiddle factors. See also the Cooley–Tukey FFT article.)"
Never stop thinking. Even if you are told to.
Earlier I said something about the "blind leading the blind" but I'm starting to get a glimmer of light re: FFT.
Let' play with that SUM function I described earlier. SUM(S0, S1, S2, S3, S3, S5, S6, S7)
Let's imagine that we actually implemented that SUM function in some language and that our implementation was recursive.
In pseudo code it might look like:
REAL SUM (LIST s) BEGIN IF LENGTH(S) = 1 RETURN (FIRST(S)) ; Only one element, nothing to sum, so just return it. ELSE r1 = SUM (EXTRACT_EVEN_ELEMENTS(s)) r2 = SUM (EXTRACT_ODD_ELEMENTS(s)) RETURN (r1 + r2) ENDIF ENDHere we assume a language that can handle lists of our S type things and we have functions to extract the odd and even position elements from lists.We might call this function like:
Now, if we think about how that call proceeds, we see that a SUM over 8 elements calls SUM over 4 elements which calls SUM over 2 elements, which calls SUM over 1 element which just returns the element.
Graphically we can see the calls proceeding as in the attached picture where I have put the lists involved in each call to SUM into boxes and added lines to track the progress of individual elements.
Now:
1) If you look at that picture you will see why the term "butterfly" is used. The way the paths of elements criss-cross starts to look like butterflies.
2) It seems that doing this recursive process performs the exact same swapping around of elements as the bit reversal algorithm.
Turns out that the fourier transform can be written recursively much like our SUM function. Except that each call of the function actually returns a list of values. Each of those values is in two parts, a sin part and a cosin part. Also it is required to multiply one of the "sub DFT" by a "fiddle factor".
What this means is that if one were to implement a recursive DFT then no bitreversal step is required. The ways the calls progress just does that for you:)
The bit reversal step is required when implementing the DFT non-recursively as is normally done.
I would refer you to the pseudo code for the DFT on this wiki page:http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm
You will see it looks much like my SUM function above.
You will also see that the part many people casually refer to as the "fiddle factor" is actually the guts of the thing where the multiplying by sin and cos happens. After all the rest of the function does no calculating at all.
I'm going to have a go at implementing that recursive FT pseudo code as directly as possible in C++ with list objects as it's the closest thing to a Fast Fourier Transform I almost understand at the moment.
Back to peering into the darkness....
Edit: Fixed up the diagram a bit.
someday radio tuners will operate with digital processing units. I have heard this suggested with tongue in cheek, but one can speculate.
Of course this has been done even on the Propeller with Phil Pilgrim's radio receiver on the Propeller project:
http://forums.parallax.com/showthread.php?105674-Hook-an-antenna-to-your-Propeller-and-listen-to-the-radio!-%28New-shortwave-prog&highlight=listen+radio
Suppose, we have a sinusoid with the fundamental frequency. Taking all the even samples first and the odd later, we compress a full oscillation to the first and to the second half, so we double the frequency of the signal! Doing this again with the two halves, we quadruplicate frequency, so if we do this often enough, he have the highest possible frequency!
What happens, when the origin was the highest frequency, so every period only is samples at + and - amplitude, then sorting odd even will bring 1 to the first half of the interval, -1 to the second. Further steps will change nothing, as all values in the first half are equal.
This remark is just to have a picture in the windings of the brain, what's going on. I see Chip in his video on speech generation, how he demonstrates the phase shift. Maybe, he already knows, how to explain FFT.
By combining Beans document which to me explains what we are trying to do combined with heaters description of the butterfly which to me explains how to do the FFT in a practical sense combined with the little I know and can remember, I now have had that Aha insight.
I could never really understand the FFT even though I had all the math 40 years ago and have actually programmed the FFT (well OK, actually copied and used) in several projects over the years. I feel now at least I have a chance at understanding, so well done to all but don't stop because there are some tantalizing hints about more to come.
Ha! Me to. I swear this thing is going to drive me insane, if it's not too late already.
"Now, what happens to the original signal when doing this bit-shift?"
I'm going to forget the bit reversal thing for a while, I'm just now trying convince myself I understand the recursive DFT as presented in the derivation and pseudo code on wikipedia here http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm
That approach splits the DFT and calls itself with the two halves recursively. Recursion stops when the is only one element in the set. This means it needs no bit revesal stage that the optimized, "in place" implementations use. As I pointed out with my SUM() game, the way the sample lists are split and passed around does all the reordering anyway.
I can't answer that but already I can see we are looking at things the opposite way around again:)
At the top level a DFT of 8 samples can have a max frequency of 4.
Split to two DFTs of 4 samples each can have a max frequency 2
Spilt to four DFT's of 2 each can have a max frequency of 1
Spilit to 8 DFT's of a single sample each can only have a max frequency of 0 or DC !
Which is exactly what the pdeudo code on wiki does, for a 1 sample DFT it returns a single value, a DC level.
Still banging head against wall....
By the way this is an interesting article: http://www.librow.com/articles/article-10
There he plays a game like my SUM() game but writing out an 8 sample DFT mathematically.
I believe everything we need to understand is in that paper, if you are happy with all that e,i,pi stuff. Just now it's making my brain melt.
Glad that butterfly thing helped a bit, I was afraid it would just confuse more.
Exactly, same here. So near but so far, very frustrating.
Well as you pointed out, Beans sample program has multiple loops inside loops resulting in lots of operations in the inner loop. And then this little butterfly comes along and magically reduces the inner loop count not by a little but by a significant amount. And Viola the FFT on a slow computer suddenly becomes very doable.
I am amazed that these two young graduates had the insight to see this reordering that would make the FFT practical at a time when "computers" really were slow and microprocessors were not even a sparkle in Chips eye. Well, we really owe them a debt because all the noise our games make and audio and signal processing (and much more) is based on this algorithm.
Maybe you should wear some head gear (to protect the wall of course).
There is and interesting historical titbit on the history of FFT in this paper http://www.cs.dartmouth.edu/~rockmore/cse-fft.pdf
I think I have an answer to your question about what happens when the samples are split into ever smaller halves.
As you say it seems that if you take all the odd samples of a sine wave and put them into a half of the list then you appear to have doubled the frequency.
BUT. I'm just looking at that recursive FFT pseudo code on wiki (reproduced below). Here we can see that at each recursive call level the list is chopped up into odds and evens which end up in the first and second halves of a new list. List X goes in, list Y comes out.
Look what happens to the list length N. At each recursive call level it is divided by two.
That N is used in the "exp(...)" gobbledygook, the so called "twiddle factors"
Well, apologies to the math shy, the twiddle factor is actually our sin and cos work where we multiply samples by our "test" sin/cos waves because:
exp(−2πik/N) = cos(2πk/N) - i*sin(2πk/N)
See that dividing N by 2 actually doubles the frequency of the test sin and cos.
Conclusion: Whilst splitting the samples into ever smaller halves we are also doubling the frequency of the sin and cos waveforms we are using as tests for the presence of frequencies. So it all works out:)
Is this a good bit of hand waving or not? Or even correct logic?
Y0,...,N−1 ← ditfft2(X, N, s): DFT of (X0, Xs, X2s, ..., X(N-1)s): if N = 1 then Y0 ← X0 trivial size-1 DFT base case else Y0,...,N/2−1 ← ditfft2(X, N/2, 2s) DFT of (X0, X2s, X4s, ...) YN/2,...,N−1 ← ditfft2(X+s, N/2, 2s) DFT of (Xs, Xs+2s, Xs+4s, ...) for k = 0 to N/2−1 combine DFTs of two halves into full DFT: t ← Yk Yk ← t + exp(−2πi k/N) Yk+N/2 Yk+N/2 ← t − exp(−2πi k/N) Yk+N/2 endfor endifHere, ditfft2(X,N,1), computes Y=DFT(X) out-of-place by a radix-2 DIT FFT, where N is an integer power of 2 and s=1 is the stride of the input X array. X+s denotes the array starting with Xs.Of course the FFT calls return a new list rather than a single sum value. Then each element of one half of the list is multiplied by the "exp gobbledygook" (Our sin and cos test frequencies) prior to adding it to the output.
Graham
The ultimate "simple" explanation might benefit from including the links to these more complicated side references, in a similar order as found in the thread.
I'm guessing new explorers will have to read the final explanation several times, going a bit deeper each pass. So far I've read through several times; its still watering dry ground, but getting better.
This is really good work, I appreciate it and I'm sure other do and more will.
The "F" in "FFT" is starting to take on a whole other meaning for me:)
OK. I think we have uncovered a dirty little secret here.
What we have is a world full of software engineers, electronics engineers, etc etc, loaded up with appropriate degrees, and diplomas in maths, physics, electronics, computer science etc. Many of them use the Fast Fourier Transform from time to time in all kind of contexts without ever knowing how it works or how to make one. And probably without wanting to admit it, perhaps thinking that everyone else "gets it" and they don't want to appear dumb. The secret is, everyone else does not get it either:)
This is kind of odd given the ubiquity of the FFT.
Sure we all use micro-processors without knowing how to design one, and we may use other maths packages that we could not create. But I think this is a little different. The maths is all out there in the open, source codes implementing FFT can be found in many places. The original Cooley-Tukey paper is one of the most sited ever. The thing is not actually that big.
Let's do an informal survey: Stand up everyone here who fits my description of the typical engineer above. (Heater stands up). Hmm... that's a lot.
The way I see it is that at one end there is the impenetrable maths with it's integrals, exponentials, imaginary numbers, Euler identities, complex roots of unity etc ect. It is possible to understand the FT from that and get through the typical questions they might ask of a physics/electronics undergraduate. But at the same time not have a real feel for what is actually going on and for sure being a long way from the "F" in "FFT". That's me.
On the other hand there is the source code you can find everywhere. This is totally unintelligible unless you know all the little tricks that went into creating it and so does not help in understanding at all.
Then there are all the explanations of how to put the "F" into "FFT" and get from the maths to the code. Most of those seem to be written by people who either don't really want to tell you, haven't really understood the thing just regurgitated what they saw elsewhere, or do really know what they are doing but are so familiar with it all that they can't imagine the state of mind of their dumbfounded audience.
Is understanding the FFT like some kind of mental illness? Those who don't have it can't possibly imagine what it is like to be a sufferer. Those who do have it can't possibly explain what it's like:)
Well, yes, mathematics has often been described as "self analysis", you can't just learn it like a bunch of historical facts. No one can explain it to you, although they might give hints and clues. You have to see it yourself and convince yourself of the truth of things.
OK, I'm putting by head protection back on....there must be a way to explain the "F" in "FFT" this without Euler.
Yep, essential, that's what I'm interested in. Especially as I haven't trodden much of the path yet.
Yep, but not just use, it's those that have a nagging itch to know what is going on.
Side trips, yes. I'm hoping those side trips can be simple thought experiments based on common knowledge that lead step by step to an A-ha moment. Like my attempt at showing how bit reversal comes about naturally from the SUM() game. That was part of the magic of Beans CORDIC paper, start with some very familiar ideas, explore them a bit, then put them into the final picture.
If you have any sticking points in Beans document you should speak up. Perhaps Bean could smooth the way in a later version. Often one forgets where one got stuck when one gets to the end of the journey.
I think at this point it is time for a first draft, a target if you will, that is going to take a few broadsides maybe even get sunk.
I would suggest this obvious layout:
A statement of what we are doing and why
Then hit em with Beans graphs and description.
Include some pseudo code (and perhaps some spin in an appendix to pay homage to the prop)
Next we state explicitly the problem ie the magnitude of the computation in a practical sense.
Now we hit em with Heaters description of the butterfly
put the two together and expand the guts (pseudo code?)
Now introduce recursion and show the connection to the butterfly
a summary complete with pseudo code and Spin and if someone is crazy enough Assembly
Ok guys let the shooting begin.
Ray
Just read through the thread and I now feel we need to expand the lead in to include some digital sampling and introduce the digital sampling concept along with enough trig to show how to express the sample points on a graph. this can be further expanded when we get to the "guts" part.
I've done it, reached enlightenment, tasted nirvana, come to the end of a long road. My guru status is assured (In my mind anyway). I am not a dummy, I am a free man. The world is mine. Is that Dr Evil you hear laughing? No, it's me:)
Well, OK, what actually happened is that I managed to:
a) Understand enough of the maths of the DFT to convince myself I know how it works. No, feel how it works.
b) Understand enough of the Cooley-Tukey idea to know why splitting the DFT up into halves repeatedly works, and how to do it.
c) Understand how the pseudo code on wikipedia for the recursive splitting works and that it does actually do what is said of it.
d) Now the biggy, actually implement that pseudo code as a real program in C++.
And it works, it does just as Cooley and Tukey and Gauss said it would! The chain is complete.
Calming down a bit, this may be the end of a long road, it's taken me over twenty years to get this far. But it may also be the start of another long road. However I'm savouring the view from here with a large can of beer.
Ahead is:
1) How to rewrite this in C with out using the short hand provided by the "complex" class of C++.
2) This is recursive and requires separate input and output buffers. How to get rid of the recursion and do all the work "in-place" to save memory.
3) Creating a Spin FFT and then perhaps a PASM FFT.
And worst of all, how to create a paper describing how this all works, I still don't see any escape from the heavy maths for the un-mathematical audience. Which is what this is all about. Perhaps "Fast Fourier Transform For Advanced Dummies" where all the hairy maths is illuminated in minute detail but keeping the scope as narrow as possible
Attached is my effort and a plot showing firstly 1024 samples of a sine wave followed by 512 frequency bins containing a peek at frequency = 16.
By the way someone has already posted an FFT in PASM here a while back, forget who it was, just wondering why he has not piped up on this thread.
I just saw your post and while I'm not quite at your place yet I can see your campfire and hear the fireworks.
The text document I just posted does not pretend to get too deep into the math . I think my target audience is an intelligent engineer, like all of us here, who because of need or pride just has to go a few steps deeper without climbing all the stairs in the ivory tower. My intention is kinda like I treat my car. I like knowing how it works but grease on my hands is icky and not for me.
I made a mistake re: the definition of "butterfly". I have to offer my apologies for any confusion and try to make amends.
Bean posted a question about the "butterfly". I posted an explanation, a revelation for me at the time, of how the magical bit-reversal stage works in the common FFT implementations. I used the SUM() game as an example and referred to it as "the butterfly".
This is not common meaning of "butterfly" in FFT circles.
Turns out that whilst the samples are shuffled around by the bit-reversal stage of a common FFT, producing "butterfly" like patterns, as the actual work progresses the data is being constantly re-shuffled such that it is in the correct order when done. It is this re-shuffling, during the work stage, that is commonly referred to as "the butterfly".
What can I say, sorry, I was was misled in that both shufflings are called butterfly but the later use is far more common in everything I have read since. For example see the flow chart in chapter 12 of the DSP Guide: http://www.dspguide.com/ch12/2.htm
This was brought to my attention by RickB.
There will be a slight delay while this dummy reattaches the missing parts.
I had not realized that the decomposition had to be undone so this is another deeper aha moment for me. I suppose that if I had gotten some grease on my hands by actually trying to implement the algorithm I would have realized I was missing something. But that is what this process is all about.
so yippie - i think.
Post edit - time passes -
So as a result of my misunderstanding I have realized two important things.
first I realized that when we finish the reordering and have the scrambled array that we have shifted from the time domain to the frequency domain. Recall from Beans diagrams that if we have an Nhz waveform and multiply it by an Nhz waveform the result is N or in other words nothing has to be done to treat it as in the frequency domain.
The second thing I realize with your help of course, is that the reverse reordering also takes place in the frequency domain such that the final array is back in the correct sequence and is the frequency spectrum, and that the butterfly (where the guts live) takes place in the frequency domain.
Your reference to chapter 12 of the Scientists and engineers guide is what did the trick for me.
My head hurts but it feels good. I'll modify the first draft text shortly...
Don't forget all the slaves were killed after they said that