Thread: 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:

Attachment 9674

Here is the graph I should get as output:

Attachment 9672

And here is what I really am getting as output:

Attachment 9673

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?

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;
}

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.

genius! thanks these simple mistakes get you every time.

it works

Originally Posted by DavidP

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.

Might I ask what library it is you are using to do the graphing in, in C++?

Originally Posted by brewbuck

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.

Originally Posted by IceDane

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.

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;
    }
}

Originally Posted by DavidP

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.