1

我正在尝试为 FFT 编写 Cooley Tukey 算法。现在，该算法运行良好，但仅适用于 2 个数字 - 仅此而已。例如，我使用了在线 FFT 计算，输入相同的数据并得到相同的结果。这是算法的代码：

#include <iostream>
#include <math.h>
#include <vector>
#include <complex>

using namespace std;

#define pi 3.14159

void ComplexFFT(vector<float> &realData, vector<float> &actualData, unsigned long sample_num, unsigned int sample_rate, int sign)
{
    unsigned long n, mmax, m, j, istep, i;
    double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;

    // CHECK TO SEE IF VECTOR IS EMPTY;

    actualData.resize(2*sample_num, 0);

    for(n=0; (n < sample_rate); n++)
    {
        if(n < sample_num)
        {
            actualData[2*n] = realData[n];
        }else{
            actualData[2*n] = 0;
            actualData[2*n+1] = 0;
        }
    }

    // Binary Inversion
    n = sample_rate << 1;
    j = 0;

    for(i=1; (i< n); i+=2)
    {
        if(j > i)
        {
            swap(actualData[j - 1], actualData[i - 1]);
            swap(actualData[j], actualData[i]);
        }
         m = sample_rate;
         while (m >= 2 && j > m) {
          j -= m;
          m >>= 1;
         }
         j += m;
     }
     mmax=2;

     while(n > mmax) {

        istep = mmax << 1;
        theta = sign * (2*pi/mmax);
        wtemp = sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi = sin(theta);
        wr = 1.0;
        wi = 0.0;

        for(m=1; (m < mmax); m+=2) {
            for(i=m; (i <= n); i += istep)
            {
                j = i + mmax;
                tempr = wr*actualData[j-1]-wi*actualData[j];
                tempi = wr*actualData[j]+wi*actualData[j-1];

                actualData[j-1] = actualData[i-1] - tempr;
                actualData[j] = actualData[i]-tempi;
                actualData[i-1] += tempr;
                actualData[i] += tempi;
            }
            wr = (wtemp=wr)*wpr-wi*wpi+wr;
            wi = wi*wpr+wtemp*wpi+wi;
        }
        mmax = istep;
    }

    // determine if the fundamental frequency
    int fundemental_frequency = 0;
    for(i=2; (i <= sample_rate); i+=2)
    {
        if((pow(actualData[i], 2)+pow(actualData[i+1], 2)) > pow(actualData[fundemental_frequency], 2)+pow(actualData[fundemental_frequency+1], 2)) {
            fundemental_frequency = i;
        }

    }
}
int main(int argc, char *argv[]) {

    vector<float> numbers;
    vector<float> realNumbers;

    numbers.push_back(50);
    numbers.push_back(10);
    numbers.push_back(50);
    numbers.push_back(10);
    ComplexFFT(numbers, realNumbers, 4, 8, -1);

    for(int i=0; (i < realNumbers.size()); i++)
    {
        cout << realNumbers[i] << ' ';
    }
}

如果我输入：

50, 10

那么结果将是：

60.00, 0.0
40.00, 0.0

而在计算器中，完全相同。

好像我进入的地方：

60.00, 0.0
40.00, 0.0

计算器返回结果：

120.0, 0.0
0.0, 0.0
80.0, 0.0
0.0, 0.0

我的结果返回：

60, 60
24.6446, -45.3553
90, -10.0001
4.64479, 45.3555

有没有人有什么建议？

PS 该算法是否仅适用于 2 位数值？

谢谢 ：）

4

1 回答 1

0

您的代码中有一些错误：

         for(m=1; (m < mmax); m+=2) {
            for(i=m; (i <= n); i += istep)
            {
                j = i + mmax;      // !!! Here, j will be bigger than actualData.size() - 1 
                tempr = wr*actualData[j-1]-wi*actualData[j];
                tempi = wr*actualData[j]+wi*actualData[j-1];

                actualData[j-1] = actualData[i-1] - tempr;
                actualData[j] = actualData[i]-tempi;
                actualData[i-1] += tempr;
                actualData[i] += tempi;
            }
            wr = (wtemp=wr)*wpr-wi*wpi+wr;
            wi = wi*wpr+wtemp*wpi+wi;
        }
于 2012-12-07T09:42:51.417 回答