14

我正在使用 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()取自“数值食谱 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时,时域数据才与原始数据相同?输入到前锋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

再次感谢!

4

1 回答 1

2

几个可能的错误:

   Data[0] == 0.0;

大概应该是===

fft_basic(数据, fftSize, InverseTransform);

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

    //Swap the fft halfes

正如我在另一个问题中告诉您的那样,如果您确实交换了它们,则需要在逆变换之前将它们交换回来。你执行这个交换吗?

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

您是否标准化 FFT 结果?离散 FFT 变换需要归一化。

于 2012-08-21T20:36:30.093 回答