FFT 算法中的一个错误

A bug in FFT Algorithm

本文关键字：一个 错误 算法 FFT      更新时间：2023-10-16

我的朋友们，我现在正在编写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线(

谢谢你帮助我！

替换

F.push_back(samples[0]);

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

F[0] = samples[0];

您将使用一个向量结束递归，该向量以零作为其第一个元素，样本作为其第二个元素。
由于您后来只使用第一个元素，因此所有内容都变为零。

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

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

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

至于你的FFT计算，我会更仔细地研究你分配F[k]和F[k+M]的行。在这些中，您将值乘以ODFT[k]，我怀疑ODFT[k]可能在那里0。

相关文章：