我正在使用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会产生与输入完全相同的脉冲:
然而,将正向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
再次感谢!