在C++中的FFT和IFFT

14

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

我的想法是将输出结果使用正向FFT转换为频域,并对相位进行线性最佳拟合(首先将其展开),然后从相位信息中减去这个最佳拟合。

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

我曾尝试在MATLAB中实现它,但没有成功,因此我转而使用C++。正向FFT工作正常,我从Numerical recipes in C++中获取了基本的配方,并使用一个函数针对复杂输入进行修改:

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;
}

使用从“Numerical recipes C++”中获取的函数ftt_basic()

我的问题是输入形式似乎影响了反FFT的输出。这可能是一个精度问题,但我搜索了一下,似乎没有其他人受到过影响。

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

enter image description here

然而,将正向FFT的功率输出取为real^2+imag^2,并将其复制到数组中,则为:

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

接下来将其用于反向FFT的输入,结果如下图所示:enter image description here

最后,将正向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)。这些值在转换回时域时产生了以下功率输出结果: enter image description here 更细致地观察脉冲结构: enter image description here 第一张图片显示出了美丽、规则的脉冲。我的问题是:是cos和sin函数的精度导致反向FFT的输出变成这样吗?为什么不同的复杂数据输入方法之间存在如此大的差异?而且只有当数据直接反馈到反向FFT中时,时域中的数据才与输入到正向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

再次感谢!


3
出于好奇,您在C ++中能做什么,在MATLAB中做不到? - Hammer
2
@KRS-fun,如果您在评论回复中添加(at)用户名,就会通知该用户您已回复,并且他们有更高的可能性重新访问问题。 - Hammer
3
这对我来说看起来不仅仅是数值不稳定的问题 - 当你写 phase[i]=atan(imag/real) 时,你不想使用 atan2(imag, real) 吗? - Eric
1
@Eric,确实是这样,谢谢你发现了这个问题!我已经使用了atan2(imag,real)函数,现在输出结果与我使用complex(real,imag)格式和complex(Amplitudecos,Amplitudesin)作为反向FFT的输入时产生的结果相同。然而,我仍然很难理解为什么似乎存在两个脉冲周期,并且为什么峰值幅度比第一个图形要小得多? - KRS-fun
1
我同意这在数学上是正确的,但我想知道你的fft和reverse_fft函数的函数签名和返回值定义。fft返回的是能量还是系数? reverse_fft是要将能量还是原始系数作为输入?如果您可以展示调用fft()的完整代码集,进行变换并调用reverse_fft(),那么我们更有可能发现问题。 - Eric
显示剩余11条评论
1个回答

2

可能存在的几个错误:

   Data[0] == 0.0;

也许应该是=,而不是==

fft_basic(Data,fftSize,InverseTransform);

我无法访问此函数的源代码;它是否真的期望您提供的布局在奇数位置具有实部,在偶数位置具有虚部?

    //Swap the fft halfes

正如我在另一个问题中告诉过你的那样,如果你交换了它们,在进行逆变换之前需要将它们交换回来。你是否执行了这个交换操作?

你的数据是否符合fft_basic函数的期望?可能它期望fftSize是2的幂。

你是否对FFT结果进行了归一化处理?离散FFT变换需要归一化。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接