程式師世界 >> 編程語言 >> 更多編程語言 >> 編程綜合問答 >> 變換-fftw使用的問題頻域抽取失敗高手請進~

變換-fftw使用的問題頻域抽取失敗高手請進~

編輯：編程綜合問答

fftw使用的問題頻域抽取失敗高手請進~

int FFTw_IFFTw_Fun2(IN dDataArray* pRRIData, IN float fmin, IN float fmax, OUT dDataArray& VLFData)
{
try
{
int nSamples = pRRIData->m_nSamples;
int N = pRRIData->m_nLen;
int Nout = floor(N/2)+1;//實數據的DFT具有 Hermitian對稱性，所以只需要計算n/2+1（向下取整）個輸出就可以了。比如對於r2c，輸入in有n個數據，輸出out有floor(n /2)＋1個數據。

    fftw_complex* out1;
    fftw_complex* out2;
    double* out3;
    out1 =  (fftw_complex *)fftw_malloc(sizeof(double)*Nout * 2);
    out2 =  (fftw_complex *)fftw_malloc(sizeof(double)*Nout * 2);
    out3 =  (double *)fftw_malloc(sizeof(double)*N);

    /////////////////////////fft變換/////////////////////////////////////////////////
    fftw_plan p = fftw_plan_dft_r2c_1d(N, pRRIData->m_pData, out1,  FFTW_ESTIMATE);     
    fftw_execute(p);// 執行變換 ... 
    fftw_destroy_plan(p);

    /////////////////////////頻域抽取/////////////////////////////////////////////////
    float f1 = fmin;
    float f2 = fmax;

    float K = 0.0f;
    for (int m=0; m<Nout-1; m++)
    {
        K=m/(Nout/nSamples);
        if ((K < f1 || K > f2) || (K > (nSamples - f2) && K < (nSamples - f1)))//1/dt為一個頻率周期
        {
            out2[m][0] = 0;
            out2[m][1] = 0;
        }
        else 
        {
            out2[m][0] = out1[m][0];
            out2[m][1] = out1[m][1];
        }
    }

    //////////////////////////////////IFFT變換////////////////////////////////////////
    fftw_plan  c2cP;
    c2cP = fftw_plan_dft_c2r_1d( N, out2, out3,  FFTW_ESTIMATE);//需注意FFTW的逆變換沒有除以N，即數據正變換再反變換後是原始數據的N倍。
    fftw_execute(c2cP);
    fftw_destroy_plan(c2cP);

    VLFData.m_nLen = N;
    VLFData.m_nSamples = 2;
    VLFData.m_pData = new double[VLFData.m_nLen];
    memset(VLFData.m_pData,0,sizeof(double)*VLFData.m_nLen);

    for (int i=0; i< N; i++)
    {
        VLFData.m_pData[i] = out3[i]/N;//取實數
    }

    fftw_free(out1);
    fftw_free(out2);
    fftw_free(out3);
    TRACE(_T("int  FFTw_IFFTw_Fun2 函數返回1~"));
    return 1;
}
catch (...)
{
    TRACE(_T("int  FFTw_IFFTw_Fun2 函數catch到錯誤~"));
    return -1;
}

}

        float fmin = 0.0f;
        float fmax = 0.0033f;
        dataOutput.VLFData.Init();
        FFTw_IFFTw_Fun2(&dataOutput.RRIData, fmin, fmax, dataOutput.VLFData);

出來的波形跟原始數據波形極其相似，我用matlab頻域抽取的話，是正確的，用C++實現就出了問題。