问题描述

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++实现就出了问题。

解决方案

找到了！这一行出了问题，K=m/(Nout/nSamples); C++默认=右侧为整形数据 整形赋值给float类型，float没有起到相应的作用，仍然为整形数据。所以抽取频域数据的时候出现了问题~改为K=1.0*m/(Nout/nSamples); 即可

解决方案二：

自己顶~