Having issues with my FFT

This is a discussion on Having issues with my FFT within the C++ Programming forums, part of the General Programming Boards category; I coded up an FFT, and I'm having issues with it. I can't figure out what the problem is. It's ...

  1. #1
    l'Anziano DavidP's Avatar
    Join Date
    Aug 2001
    Location
    Plano, Texas, United States
    Posts
    2,738

    Having issues with my FFT

    I coded up an FFT, and I'm having issues with it. I can't figure out what the problem is. It's not the 1st time I've coded the FFT, and I've compared my new FFT code (in C++) to a previous time I coded the FFT (in C#), and I don't personally see any mistakes in the way I've coded it (although the two implementations have some minor semantic differences).

    Here is the graph of the signal I am sending into the function:

    Name:  signal.PNG
Views: 473
Size:  7.6 KB

    Here is the graph I should get as output:

    Name:  dft_result.PNG
Views: 356
Size:  4.5 KB

    And here is what I really am getting as output:

    Name:  fft_result.PNG
Views: 451
Size:  9.6 KB

    The output off the fft looks oddly like the sync function...which shouldn't be the result for a simple sine wave...

    Here is my code:
    Code:
    std::vector < std::complex<double> > InternalFastFourierTransform (std::vector < std::complex<double> > input)
    {
        static std::complex<double> i = std::complex<double>(0, 1);
        static std::complex<double> NegativeTwoPiI = -2.0 * PI * i;
    
        int N = (int)input.size();
        int halfN = N / 2;
    
        if (N == 1)
        {
            return input;
        }
        else
        {
            //Recurse and get evens/odds
            std::vector< std::complex<double> > evens = InternalFastFourierTransform( GetVectorRange(input, 0, N, 2) );
            std::vector< std::complex<double> > odds = InternalFastFourierTransform( GetVectorRange(input, 1, N, 2) );
    
            //Combine
            std::vector< std::complex<double> > result = std::vector< std::complex<double> > (N);
            for(int k = 0; k < halfN; k++)
            {
                std::complex<double> exponent = std::exp(NegativeTwoPiI * double(k / N));
                result[k] = evens[k] + exponent * odds[k];
                result[k + halfN] = evens[k] - exponent * odds[k];
            }
    
            return result;
        }
    }
    Do you guys see anything wrong? The only thing that I could see is that possibly my function GetVectorRange which takes a slice of a vector from some starting index to some ending index at a stride of 2 could be functiong incorrectly...I haven't thoroughly tested that function...but is that really the problem, or do you think it could be some math mistake based off of the graph output?
    Last edited by DavidP; 04-05-2010 at 10:03 AM. Reason: grammatical correction
    My Website

    "Circular logic is good because it is."

  2. #2
    l'Anziano DavidP's Avatar
    Join Date
    Aug 2001
    Location
    Plano, Texas, United States
    Posts
    2,738
    I might as well include the code for the GetVectorRange function if that is suspicious as being the culprit, although I don't think it is. Here it is anyways:

    Code:
    std::vector< std::complex<double> > GetVectorRange(std::vector< std::complex<double> > input, int startIndex, 
    int endIndex, int stride)
    {
        std::vector< std::complex<double> > output;
        for(int k = startIndex; k < endIndex && k < (int)input.size(); k += stride)
        {
            output.push_back(input[k]);
        }
    
        return output;
    }
    Last edited by DavidP; 04-05-2010 at 10:29 AM.
    My Website

    "Circular logic is good because it is."

  3. #3
    Just a pushpin. bernt's Avatar
    Join Date
    May 2009
    Posts
    426
    Code:
    std::complex<double> exponent = std::exp(NegativeTwoPiI * double(k / N));
    I have completely no idea how the math works, but this part stood out as a possible culprit. I think it ought to be
    Code:
    std::complex<double> exponent = std::exp(NegativeTwoPiI * double(k)/double(N));
    if you want a decimal rather than an int there.
    Consider this post signed

  4. #4
    l'Anziano DavidP's Avatar
    Join Date
    Aug 2001
    Location
    Plano, Texas, United States
    Posts
    2,738
    genius! thanks these simple mistakes get you every time.

    it works
    My Website

    "Circular logic is good because it is."

  5. #5
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,249
    Quote Originally Posted by DavidP View Post
    genius! thanks these simple mistakes get you every time.

    it works
    Good deal -- those complex roots of unity ("twiddle factors") are going to be constant across all transforms of the same size. You might as well precompute them and stash them away somewhere instead of calling complex std::exp() in that inner loop.
    Code:
    //try
    //{
    	if (a) do { f( b); } while(1);
    	else   do { f(!b); } while(1);
    //}

  6. #6
    Ex scientia vera
    Join Date
    Sep 2007
    Posts
    478
    Might I ask what library it is you are using to do the graphing in, in C++?
    "What's up, Doc?"
    "'Up' is a relative concept. It has no intrinsic value."

  7. #7
    l'Anziano DavidP's Avatar
    Join Date
    Aug 2001
    Location
    Plano, Texas, United States
    Posts
    2,738
    Quote Originally Posted by brewbuck View Post
    Good deal -- those complex roots of unity ("twiddle factors") are going to be constant across all transforms of the same size. You might as well precompute them and stash them away somewhere instead of calling complex std::exp() in that inner loop.
    Agreed. I need to study the concept of roots of unity more. When I took an algorithms class about 2 years ago, we discussed the complex roots of unity, and after lots of studying I had a basic, but useable, understanding of them. Now I've lost much of my understanding of the roots of unity...so I need to go back over them.

    Quote Originally Posted by IceDane View Post
    Might I ask what library it is you are using to do the graphing in, in C++?
    I'm actually not outputting the graphs in my C++ program in this case. I output my array data to file, and then I open up a program called Octave (an open source clone of Matlab), and use its plot function to do the plotting of the data for me. Octave uses gnuplot as its plotting utility.
    My Website

    "Circular logic is good because it is."

  8. #8
    l'Anziano DavidP's Avatar
    Join Date
    Aug 2001
    Location
    Plano, Texas, United States
    Posts
    2,738
    Okay, even though the problem posted at the beginning of the thread has already been solved, in follow up to brewbuck's suggestion, I did some research on roots of unity and refreshed myself a little bit.

    I believe this optimizes out the call to std::exp correctly:

    Code:
    std::vector < std::complex<double> > InternalFastFourierTransform (std::vector < std::complex<double> > input)
    {
        static std::complex<double> i = std::complex<double>(0, 1);
        static std::complex<double> NegativeTwoPiI = -2.0 * PI * i;
    
        int N = (int)input.size();
        int halfN = N / 2;
    
        std::complex<double> omega_n = std::exp(NegativeTwoPiI / double(N));
        std::complex<double> omega = std::complex<double> (1, 0);
    
        if (N == 1)
        {
            return input;
        }
        else
        {
            //Recurse and get evens/odds
            std::vector< std::complex<double> > evens = InternalFastFourierTransform( GetVectorRange(input, 0, N, 2) );
            std::vector< std::complex<double> > odds = InternalFastFourierTransform( GetVectorRange(input, 1, N, 2) );
    
            //Combine
            std::vector< std::complex<double> > result = std::vector< std::complex<double> > (N);
            for(int k = 0; k < halfN; k++)
            {
                result[k] = evens[k] + omega * odds[k];
                result[k + halfN] = evens[k] - omega * odds[k];
                omega = omega * omega_n;
            }
    
            return result;
        }
    }
    My Website

    "Circular logic is good because it is."

  9. #9
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,249
    Quote Originally Posted by DavidP View Post
    Agreed. I need to study the concept of roots of unity more. When I took an algorithms class about 2 years ago, we discussed the complex roots of unity, and after lots of studying I had a basic, but useable, understanding of them. Now I've lost much of my understanding of the roots of unity...so I need to go back over them.
    It might be easier to imagine them as phasors, little stopwatches where the second hand rotates around some integral number of revolutions -- the value of k dictates how quickly the phasor hand rotates. Essentially, what it's doing is correlating your signal with a harmonic basis vector, and the number of revolutions of the phasor goes up by integral amounts for each basis -- the first one makes one full revolution over the transform period, the second makes two revolutions, etc.
    Code:
    //try
    //{
    	if (a) do { f( b); } while(1);
    	else   do { f(!b); } while(1);
    //}

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. I hate pointers or Pointer Issues
    By bolivartech in forum C Programming
    Replies: 9
    Last Post: 11-14-2009, 10:48 AM
  2. FFT and convolution
    By DavidP in forum General Discussions
    Replies: 7
    Last Post: 07-03-2009, 09:31 AM
  3. FFT returning values in wrong order
    By DavidP in forum C# Programming
    Replies: 3
    Last Post: 10-25-2007, 01:15 PM
  4. hexdump issues
    By daluu in forum C Programming
    Replies: 2
    Last Post: 03-04-2003, 08:01 PM
  5. FFT discussion, anyone?
    By Sebastiani in forum A Brief History of Cprogramming.com
    Replies: 4
    Last Post: 12-26-2002, 03:06 AM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21