0

我的朋友们，我现在正在编写 FFT 算法，我有以下程序：

#include<iostream>
#include<vector>
#include<complex>
#include<cmath>
#define _USE_MATH_DEFINES
#include<math.h>
#include<stdlib.h>

using namespace std;

const complex<double> I(0,1);
const int pi = M_PI;

//This function will check if a number is a power of certain number 
bool checkpower(float base,float num){
    float a = ceil(log(num)/log(base))-floor(log(num)/log(base));
    if (a==0){
        return 1;
    }
    else{
    return 0;
    }
}


//Fast Fourier Transform for DFT
vector<complex<double>> FFT_DFT(vector<complex<double>> samples){
    int N = samples.size();
    cout << N << endl;
    if(!checkpower(2,N)){
        cout << "Please make sure the sample size is of 2^n!" << endl;
        exit (EXIT_FAILURE);
    }
    else{
        vector<complex<double>> F(N);
        if(N==1){
            F.push_back(samples[0]);
        }
        else{
            int M = N/2;
            cout << M << endl;
            vector<complex<double>> O; //vector to store the odd elements
            vector<complex<double>> E; //vector to store the even elements
            //Reoder the samples
            for(int l=0;l<M;l++){
                E.push_back(samples[2*l]);
                O.push_back(samples[2*l+1]);
            }
            vector<complex<double>> ODFT;
            cout << "Start recursive for odd" << endl;
            ODFT = FFT_DFT(O);
            vector<complex<double>> EDFT;
            cout << "Start recursive for even" << endl;
            EDFT = FFT_DFT(E);
            for(int k=0;k<M;k++){
                cout << real(EDFT[k]) << " + "<< imag(EDFT[k]) << "I" << endl;
                cout << real(ODFT[k]) << " + "<< imag(ODFT[k]) << "I" << endl;
                F[k] = EDFT[k]+exp(-2.0*pi*k*I/(double)N)*ODFT[k];
                F[k+M] = EDFT[k]-exp(-2.0*pi*k*I/(double)N)*ODFT[k]; 
                cout << real(F[k]) << " + "<< imag(F[k]) << "I" << endl;
                cout << real(F[k+M]) << " + "<< imag(F[k+M]) << "I" << endl;
            }
        }
        return F;
    }
}


int main(){
    vector<complex<double>> samples;
    samples.push_back(8.0);
    samples.push_back(4.0);
    samples.push_back(8.0);
    samples.push_back(0.0);

    vector<complex<double>> dft = FFT_DFT(samples);
    vector<complex<double>>::iterator item;

    for(item=dft.begin();item!=dft.end();item++){
        cout << real(*item) << " + "<< imag(*item) << "I" << endl;
    }


    return 0;
}

我使用 Visual Studio 2010 Professional 作为编译器。我不知道我的递归算法出了什么问题，所以我在 VS 中使用了调试模式，并逐行检查了过程，但它似乎总是给我所有值的 0。我用普通的 FFT 算法对其进行了测试，效果很好。那么，谁能帮我看看我的程序？我已经调试了大约 4 个小时，但仍然找不到错误。也许，我做了一些非常愚蠢的事情而我没有注意到。

（只是为了您的注意，在 FFT 函数中，为了看看到底发生了什么，我还添加了许多 cout 行）

谢谢你的协助！

4

2 回答 2

1

代替

F.push_back(samples[0]);

它将一个元素添加到F- 已经有一个元素，零 - 与

F[0] = samples[0];

您将使用一个向量结束递归，该向量的第一个元素是零，第二个元素是样本。
由于您稍后只使用第一个元素，因此一切都变为零。

于 2013-07-30T11:38:05.657 回答
0

我不能立即说出你的问题是什么，但你至少有一个严重的错误。这里：

float a = ceil(log(num)/log(base))-floor(log(num)/log(base));
if (a==0){
    return 1;
}

除了有效的比较是 to 之外0.0f，您不应该进行一些浮点计算并期望该值是一个特定的数字。在您期望的地方0.0，它实际上可能是这样的0.00000000132。因为这就是浮点数的工作方式。不幸的是，他们总是会这样工作。

至于您的 FFT 计算，我会更仔细地查看您分配F[k]和的行F[k+M]。在那些中，你乘以一个值ODFT[k]，我怀疑ODFT[k]可能0存在。

于 2013-07-30T10:55:46.143 回答