2
$\begingroup$

​ ​Hello, I would like to ask you for help in understanding Fast Fourier Transform.

Most articles about FFT describe a simple DFT example with N=8 number of samples. They divide it on half, to evens and odds. And to that point everything is quite clear for me (according to that I have implemented the algorithm, and I get the expected outcome for the whole range of frequencies - in that case from 1 to 8).

But then they instruct divide it again, analogous to the first divide, and then again. Until you end up with N/8. And there is something about bit reversal. But they don't show it, just say "analogous to the first divide". But it seems in my case "analogous" is still not enough to understand. Sorry for that. So for those next divides I get wrong results. I do something wrong. Don't know what. Please help me.

Please let me show you my way of thinking through a simple code example . First I start with DFT and that is clear for me:

for(int freqBin=1; freqBin <= sampleRate; freqBin++)
    {
        _Sfx[0]  = 0.0f;     // I need to make it zero to perform _Sfx[0] += …
                             // in next coming for loop   

        for (int n=0; n<sampleRate; ++n)
        {
            complex<float> _Wnk_N = exp(-1i * 2.0 * M_PI * n * freqBin / sampleRate);

            _Sfx[0] += inputSignal[n]  * _Wnk_N;
        }  

       output[freqBinK] = _Sfx[0];
    }

Then I make first divide:

  for(int freqBin=1; freqBin <= sampleRate/2; freqBin++)
    {
        for(int i=0; i<2; i++)
        {
            _Sfx[i]  = 0.0f;
        }

        for (int n=0; n<sampleRate/2; ++n)
        {
            complex<float> _Wnk_N = exp(-1i * 2.0 * M_PI * n * freqBin / (sampleRate/2.0));

            _Sfx[0] += inputSignal[2*n   ]   * _Wnk_N;
            _Sfx[1] += inputSignal[2*n +1]   * _Wnk_N;
        }

        complex<float> _Wk_N = exp(-1i * 2.0 * M_PI * freqBin / sampleRate);

        _Sf_I_half    = _Sfx[0]  + _Wk_N * _Sfx[1];

        _Sf_II_half   = _Sfx[0]  - _Wk_N * _Sfx[1];


       output[freqBinK]                 = _Sf_I_half;
       output[freqBinK + sampleRate /2] = _Sf_II_half;
    }

and that is also clear for me - not for sure but I think so, because results are as expected.

But then I make next divide like that:

for(int freqBinK=1; freqBinK <= sampleRate/4; freqBinK++)
    {
        for(int i=0; i<4; i++)
        {
            _Sfx[i]  = 0.0f;
        }

        for (int n=0; n<bufferSize/4; ++n)
        {
            complex<float> _Wnk_N = exp(-1i * 2.0 * M_PI * n * freqBin / (sampleRate/4.0));

            _Sfx[0] += inputSignal[ 2*(2*n)       ]     * _Wnk_N ;
            _Sfx[1] += inputSignal[ 2*(2*n) +1    ]     * _Wnk_N ;
            _Sfx[2] += inputSignal[ 2*(2*n  +1)   ]     * _Wnk_N ;
            _Sfx[3] += inputSignal[ 2*(2*n  +1) +1]     * _Wnk_N ;
        }

        complex<float> _Wk_N   = exp(-1i * 2.0 * M_PI * freqBin / sampleRate);
        complex<float> _Wk_N2 = exp(-1i * 2.0 * M_PI * freqBin / (sampleRate/2));

        _Sf_I_quarter    =              (_Sfx[0]  + _Wk * _Sfx[1]) +
                             _Wk_N2  *  (_Sfx[2]  + _Wk * _Sfx[3]);


        _Sf_II_quarter   =              (_Sfx[0]  - _Wk * _Sfx[1]) +
                             _Wk_N2  *  (_Sfx[2]  - _Wk * _Sfx[3]);


        _Sf_III_quarter  =              (_Sfx[0]  + _Wk * _Sfx[1]) -
                             _Wk_N2  *  (_Sfx[2]  + _Wk * _Sfx[3]);


        _Sf_IV_quarter   =              (_Sfx[0]  - _Wk * _Sfx[1]) -
                             _Wk_N2  *  (_Sfx[2]  - _Wk * _Sfx[3]);



        output[freqBin]                     = _Sf_I_quarter;
        output[freqBin +     sampleRate /4] = _Sf_II_quarter;
        output[freqBin + 2 * sampleRate /4] = _Sf_III_quarter;
        output[freqBin + 3 * sampleRate /4] = _Sf_IV_quarter;
    }

and for that code I get only half range as expected, but above Nyquist (more than N/2) there are some unexpected values. Why? What did I do wrong in that last example of code?

Could you modify it for me and put it as it should be? Please don’t write about c++ programming. When I asked that question on Stackoverflow, most answers were about code, but I am asking about FFT. And please don't show me any algorithm to do complete FFT, I just want to understand just that next step, after divide DFT on half. If I understand it I hope I will be able to create algorithm to compute whole range FFT.

And last thing: of course I know there are lot of solutions for FFT, for free, ready to use. But my goal is not to make FFT, but to understand it.

For any help great thanks in advance. Best regards

EDIT: There was question what exactly outputs I mean when I say "expected" or "unexpected".

So let's say I have input signal like that:

for(int sample=0; sample<8; sample++)
{
    inputSignal[sample] = sinf(1.0 * (float)sample * 2.0* M_Pi / 8.0);
}

When I run that signal through the DFT I get:

output[1] = 4.0;
output[2] = 0.0; 
output[3] = 0.0;
output[4] = 0.0;
output[5] = 0.0;
output[6] = 0.0;
output[7] = 4.0;
output[8] = 0.0;

I need to comment it:

  • All zeros are approximate. Actually to be more precise all zeros are something like 8.52367e-07. But it's not the point, that's why I approximated it.
  • As you can see I start output from 1, not from zero. Ask why? I make my DFT frequency loop also from 1. Because I am not interested in frequency 0 Hz. And those numbers correspond to
    frequency. So for random reader of that question I thought it would
    be more clear that output[1] is for 1 Hz. But actually it's also not the point.

But to the point. Now you can see my outputs. And in my opinion they are expected values. output[1] is more than 0.0, because my input signal has a 1 Hz sinusoid. And output[7] is also more than 0.0, because it's a mirror freq in relative to the Nyquist frequency, which is 4 Hz, and as long as I know it's expected. Do you agree?

In my second example of code (where I divide DFT for evens and odds) gives me exactly the same output values. That's why I suppose I make first "divide and conquer" in proper way.

But my next "divide and conquer" (my third example of code) gives me different values:

output[1] = 4.0;
output[2] = 0.0; 
output[3] = 0.0;
output[4] = 0.0;
output[5] = 2.82843;
output[6] = 0.0;
output[7] = 2.82843;
output[8] = 0.0;

And due to fact they are different values, I think those values are unexpected. And that's why I think I do something wrong. But can't find the answer, what exactly is wrong? Only can see that problem is for frequencies above Nyquist, but why? Please help me.

$\endgroup$
  • $\begingroup$ As a first step, I would recommend you to understand the two flavours of decimation (in time & frequency). As a second step, understand that bit reversal applies to certain cases and therefore you have to decide if you are working on integer powers of two or not. Thirdly, if you want to understand the FFT this edition of Proakis-Manolakis DSP is excellent (Chapter 6). $\endgroup$ – A_A Mar 15 '18 at 11:34
  • $\begingroup$ Hey, great thanks for reply. As I know I try to do radix-2, decimation in time. I also know what is bit reversal, but don’t know how to use that knowledge to implement FFT. Can’t you just modify my third example of code? that one where I divide DFT by 4. I believe there is only small tweaking to do. I just need that to understand it, please. How those bit reversal inputs shoud go? Multiply by what? Add to what? And where is the minus or plus. $\endgroup$ – pajczur Mar 15 '18 at 12:23
  • $\begingroup$ I might attempt a more extensive reply later on but I cannot write the code for you. $\endgroup$ – A_A Mar 15 '18 at 12:33
  • $\begingroup$ @A_A Ok I would be very grateful for any more reply. By the way, why can’t you just modify my code little bit? :) I don’t want you write the code for me. I want somebody show me mistakes in my code which I already written. Even if you modify that small part of code I would never use it in any app, because to use it first I need to implement whole range of FFT, but it’s only second divide of DFT. Usually my apps use about 44,1kHz sample rate, and buffer size about 512. So I need to do more divides. I suppose I need to divide FFT by 512. So if you show me how to divide by 4 it would help me $\endgroup$ – pajczur Mar 15 '18 at 12:44
  • $\begingroup$ @A_A By the way, thanks for any help in advance. $\endgroup$ – pajczur Mar 15 '18 at 12:44
0
$\begingroup$

This is the place where I learned what you are asking: MIT Opencourseware 6.008 Digital Signal Processing . See lecture 18. And if you have the chance also see lectures 19 and 20.

You are trying to implement the decimation in time version of the FFT algorithm. There's also a decimation in frequency explained in lecture 19.

About your bit-reversal doubt. Look at the steps:

  1. original sample indexes 0,1,2,3,4,5,6,7. You divide even and odd
  2. sample indexes {0,2,4,6};{1,3,5,7}. You divide even and odd again in every group
  3. sample indexes {0,4};{2,6};{1,5};{3,7}. You stop here

Compare original indexes 0,1,2,3,4,5,6,7 with final indexes 0,4,2,6,1,5,3,7. Take both of them to bits.

  • original indexes 000,001,010,011,100,101,110,111
  • the final indexes 000,100,010,110,001,101,011,111 .

Now reverse every bit representation in the input to get the bit representations in the output. Do you see how for instance 001 is now 100 and 100 is now 001. Of course if you reverse 000,010,101,111 you get exactly the same positions that is why samples 0,2,5 and 7 of the input remain at their original index positions in the output. This is the famous bit reversal

$\endgroup$
  • $\begingroup$ Hey, OK, but the question is what to do with those bit reversal indexes of samples. Where to put them, add to what and multiply by what. Please see my third example of code. Could you just copy it and modify how it should be please? When I see it then I think I would be able to understand it. $\endgroup$ – pajczur Mar 15 '18 at 6:21
  • $\begingroup$ There are in that code input signals with indexes: [ 2*(2*n)], [ 2*(2*n) +1], [ 2*(2*n +1) ], [ 2*(2*n +1) +1]. And n goes from 0 to 1 (n<sampleRate/4). When you put n sample by sample to those indexes you will get: 0, 4, 1, 5, 2, 6, 3, 7. OK it’s wrong order now. But when I change it to: [ 2*(2*n)], [ 2*(2*n) +2], [ 2*(2*n ) +1 ], [ 2*(2*n +1) +1] then I get order the same like you show. But it doesn’t change my results. Still wrong output I get. So I do something wrong with those bit inverted sample indexes. $\endgroup$ – pajczur Mar 15 '18 at 6:28
  • $\begingroup$ And now I’ve found out you make correction in my question. Thanks for that. English is not my first language. I still learn it :) $\endgroup$ – pajczur Mar 15 '18 at 7:48
  • $\begingroup$ I think problem with my code is somewhere around _Wnk_N, _Wk_N, _Wk_N2, and with quarters _Sf_I_quarter, _Sf_II_quarter, ... $\endgroup$ – pajczur Mar 15 '18 at 7:59
  • $\begingroup$ Don't you even have 17 minutes? See lecture 18 from minute 20 to 37! Everything is there crystal clear. The reason why you moved your question from stackoverflow to dsp is because you wanted to understand. If you want your code reviewed leave it there and they will review it. $\endgroup$ – VMMF Mar 15 '18 at 14:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.