zhaopeizhaopeipei
teresa_zp
采纳率100%
2015-04-14 05:24 阅读 3.2k
已采纳

fftw使用的问题 频域抽取失败 高手请进~

1

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

  • 点赞
  • 写回答
  • 关注问题
  • 收藏
  • 复制链接分享

2条回答 默认 最新

  • 已采纳
    zhaopeizhaopeipei teresa_zp 2015-04-16 06:13

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

    点赞 3 评论 复制链接分享
  • zhaopeizhaopeipei teresa_zp 2015-04-15 00:24

    自己顶~

    点赞 1 评论 复制链接分享

相关推荐