0

I've been attempting to code the FFT lately and have not been able to get correct results on sample problems. As I have coded recursively, debugging has been difficult. This is also my first time working with complex vectors, so it is likely that my issue may lie in my declaration of my valarrays. Below is my code (the reason I am using a template is because I plan on measuring single versus double precision to verify stability of the algorithm).

#include <iostream>
#include <complex>
#include "dftfft.h"
#include <valarray>

int main()
{
    std::valarray<std::complex<double>> V{ {-1,0}, {2,0}, {3,0}, {0,0} };
    std::valarray<std::complex<double>> VFFT(FFT<double>(V));
    for (size_t i = 0; i < VFFT.size(); i++)
    {
        std::cout << VFFT[i] << std::endl;
    }
}

In "dftfft.h":

template <typename T>
std::valarray<std::complex<T>> FFT(std::valarray<std::complex<T>> P)
{
    double pi = 2 * acos(0.0);
    size_t n = P.size();
    if (n == 1)
    {
        return P;
    }
    std::valarray<std::complex<T>> P_even = P[std::slice(0, n/2, 2)];
    std::valarray<std::complex<T>> P_odd = P[std::slice(1, n/2, 2)];
    std::valarray<std::complex<T>> y_even = FFT(P_even);
    std::valarray<std::complex<T>> y_odd = FFT(P_odd);
    std::valarray<std::complex<T>> y(n);
    for (size_t i = 0; i < n / 2; i++)
    {
        std::complex<T> omega = std::polar(1.0, 2.0 * pi / n);
        y[i] = P_even[i] + pow(omega, i) * P_odd[i];
        y[i + n / 2] = P_even[i] - pow(omega, i) * P_odd[i];
    }
    return y;
}

For this particular vector V the output is:

(1,0) (3,0) (-3,0) (3,0)

I know this to be incorrect, but cannot figure out why. The correct answer is: (4,0) (-4,-2) (0,0) (-4,2)

Any input would be highly appreciated.

Heggo
  • 1
  • I can't immediately see what's wrong but this looks very similar to the [Rosetta Code C++ FFT implementation](https://rosettacode.org/wiki/Fast_Fourier_transform#C.2B.2B). Comparing your code to that might give you some insights – alter_igel Feb 15 '21 at 02:04
  • `omega` is a constant, it should be computed outside the loop. And you compute `pow(omega, i)` twice. Instead, do `x=omega` outside the loop, and `x*=omega` inside the loop, this should be much cheaper than computing the power. – Cris Luengo Feb 15 '21 at 02:32
  • The Rosetta Code version linked in the comment above uses `-2.0 * pi / n`, you have `2.0 * pi / n`. – Cris Luengo Feb 15 '21 at 02:37

1 Answers1

0

Solved the issue. In my combining loop I needed to use y_even/odd instead of P_even/odd.

Heggo
  • 1