C++ 中的 FFT 和 IFFT答案

【问题标题】：FFT and IFFT in C++C++ 中的 FFT 和 IFFT
【发布时间】：2012-08-11 07:24:54
【问题描述】：

我正在使用 C++/C 对一些应该是激光脉冲输出的数据执行正向和反向 FFT。

这个想法是获取输出，使用前向 FFT 转换到频域，对相位应用线性最佳拟合（首先展开它），然后从相位信息中减去这个最佳拟合。

然后将得到的相位和幅度转换回时域，最终目标是通过相位补偿压缩脉冲。

我曾尝试在 MATLAB 中执行此操作，但未成功，因此转而使用 C++。前向 FFT 工作正常，我从 C++ 中的 Numerical recipes 中获取了基本配方，并使用了一个函数来修改它以适应复杂的输入，如下所示：

void fft(Complex* DataIn, Complex* DataOut, int fftSize, int InverseTransform, int fftShift)
{

        double* Data  = new double[2*fftSize+3];
        Data[0] == 0.0;


        for(int i=0; i<fftSize; i++)
        {
                Data[i*2+1]  = real(DataIn[i]);
                Data[i*2+2]  = imag(DataIn[i]);
        }

        fft_basic(Data, fftSize, InverseTransform);

        for(int i=0; i<fftSize; i++)
        {
                DataOut[i] = Complex(Data[2*i+1], Data[2*i+2]);
        }

        //Swap the fft halfes
        if(fftShift==1)
        {
                Complex* temp = new Complex[fftSize];
                for(int i=0; i<fftSize/2; i++)
                {
                        temp[i+fftSize/2] = DataOut[i];
                }
                for(int i=fftSize/2; i<fftSize; i++)
                {
                        temp[i-fftSize/2] = DataOut[i];
                }
                for(int i=0; i<fftSize; i++)
                {
                        DataOut[i] = temp[i];
                }
                delete[] temp;
        }
        delete[] Data;
}

使用函数 ftt_basic() 取自“Numerical recipes C++”。

我的问题是输入的形式似乎会影响反向 FFT 的输出。这可能是一个精度问题，但我环顾四周，之前似乎没有影响到其他任何人。

将正向 FFT 的输出直接反馈到反向 FFT 会产生与输入相同的脉冲：

但是，将前向 FFT 的功率输出设为real^2+imag^2 并将其复制到一个数组中：

Reverse_fft_input[i]=complex(real(forwardsoutput[i]),imag(forwardsoutput[i]));

然后将其用作反向 FFT 的输入，得到以下结果：

最后，获取正向 FFT 的输出并复制如下：

Reverse_fft_input[i]=complex( Amplitude[i]*cos(phase[i]), Amplitude[i]*sin(phase[i]));

其中 Amplitude[i]=(real^2+imag^2)^0.5 和 phase[i]=atan(imag/real)。转换回时域时产生以下功率输出：

仔细观察脉冲结构：

当第一张图片产生良好、规律的脉冲时。

我的问题是，是 cos 和 sin 函数的精度导致反向 fft 的输出变成这样吗？为什么复杂数据的不同输入方式存在如此巨大的差异，为什么只有直接反馈到反向FFT时，时域数据才与原始数据相同？输入到 forwads FFT 中？

谢谢。

*这里编辑是函数的实现：

void TTWLM::SpectralAnalysis()

{

    Complex FieldSpectrum[MAX_FFT];
    double  PowerFFT[MAX_FFT];
    double  dlambda;
    double  phaseinfo[MAX_FFT]; // Added 07/08/2012 for Inverse FFT
    double  fftamplitude[MAX_FFT]; // Added 07/08/2012 for Inverse FFT after correction
    double  phasecorrect[MAX_FFT]; // Added 07/08/2012 for Inverse FFT after correction
    double  lambdaarray[MAX_FFT]; // Added 07/08/2012 for Inverse FFT after correction
    Complex CompressedFFT[MAX_FFT];
    Complex correctedoutput[MAX_FFT];


    //Calc the wavelength step size
    dlambda = lambda*lambda/CONST_C/DT/fftSize;
    //Calculate the spectrum
    fft(fftFieldData, FieldSpectrum, fftSize, FORWARD, SHIFT); // Forward fft of the output data 'fftFieldData' into frequency domain

    //Get power spectrum
    for(int i=0; i<fftSize; i++)
    {

        PowerFFT[i]                  = norm(FieldSpectrum[i]);
        phaseinfo[i]                 = atan2(imag(FieldSpectrum[i]),real(FieldSpectrum[i]));
        fftamplitude[i] = sqrt(PowerFFT[i]);      // Added 07/08/2012 for Inverse FFT after correction


    }

    // Added 07/08/2012 for Inverse FFT after correction, this loop subtracts line of best fit from the phase

        for(int i=0; i<fftSize; i++)
    {
        lambdaarray[i]=dlambda*(i-fftSize/2)*1e-2;
        phasecorrect[i]=phaseinfo[i]-((1.902e+10*lambdaarray[i])+29619); // Correction from best fit in MATLAB (DONE MANUALLY) with phase unwrapping
        CompressedFFT[i]=(fftamplitude[i]*cos(phaseinfo[i]),fftamplitude[i]*sin(phaseinfo[i]));
            }

       fft(CompressedFFT, correctedoutput, fftSize, REVERSE, SHIFT); // Reverse fft of corrected phase back to time domain, final output is correctedoutput

再次感谢！

【问题讨论】：

出于好奇，有什么可以在 c++ 中做而在 MATLAB 中做不到的？
@KRS-fun，如果您将此评论中的用户名添加到评论回复中，用户名将被通知您已回复，并且将有更高的机会重新访问该问题。
在我看来，这不仅仅是数值不稳定性问题——当你写 phase[i]=atan(imag/real) 时，你不想使用 atan2(imag, real) 吗？跨度>
@Eric 我确实做到了，感谢您发现这一点！我使用了 atan2(imag,real)，现在输出的结果与我使用 complex(real,imag) 和 complex(Amplitudecos,Amplitudesin) 格式作为反向 FFT。但是我仍然很难理解为什么似乎有两个脉冲周期性，以及为什么峰值幅度比第一张图小得多？
我同意这在数学上是正确的，我想知道你的 fft 和 reverse_fft 函数的函数签名和返回值的定义。 fft 是返回功率还是仅返回系数？ reverse_fft 是否要将功率或原始系数作为输入？如果您能展示调用 fft()、转换和调用 reverse_fft() 的完整代码集，我认为我们更有可能发现问题。

标签： c++ transform fft ifft

【解决方案1】：

几个可能的错误：

   Data[0] == 0.0;

应该是=，而不是==？

fft_basic(数据, fftSize, InverseTransform);

我无权访问该函数的源代码；它真的期望您在奇数位置提供实部而在偶数位置提供虚部的布局吗？

    //Swap the fft halfes

正如我在另一个问题中告诉您的那样，如果您确实交换了它们，则需要在逆变换之前将它们交换回来。您是否执行此交换？

您的数据是否符合 fft_basic 函数的预期？可能它期望 fftSize 是 2 的幂。

您是否对 FFT 结果进行标准化？离散 FFT 变换需要归一化。

【讨论】：