БПФ и ОБПФ в C++

Я использую C++/C для выполнения прямого и обратного БПФ для некоторых данных, которые должны быть импульсным выходом лазера.

Идея состоит в том, чтобы взять выходные данные, использовать прямое БПФ для преобразования в частотную область, применить линейное наилучшее соответствие к фазе (, сначала разворачивая его ), а затем вычитая это наилучшее соответствие из информации о фазе.

Результирующие фаза и амплитуда затем преобразуются обратно во временную область, конечной целью является сжатие импульсов посредством фазовой компенсации.

Я безуспешно пытался сделать это в MATLAB и в результате обратился к C++. Прямое БПФ работает нормально, я взял базовый рецепт из числовых рецептов на 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;
}

с функцией ftt_basic(), взятой из «Численных рецептов C++».

Моя проблема в том, что форма ввода, по-видимому, влияет на вывод обратного БПФ. Это может быть проблема с точностью, но я огляделся, и, похоже, это никого раньше не затрагивало.

Подача выходного сигнала прямого БПФ непосредственно обратно в обратное БПФ дает импульсы, идентичные входным:

Однако принимая выходную мощность, взятую как real^2+imag^2прямого БПФ, и копируя ее в массив, такой, что:

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

а затем используя это в качестве входных данных для обратного БПФ, получаем следующее:

И, наконец, получение результата прямого БПФ и копирование таким образом, что:

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

где Amplitude[i]= (real^2+imag^2 )^0,5 и фаза[i]=atan (imag/real ). Дает следующую выходную мощность при обратном преобразовании во временную область:

при ближайшем рассмотрении структуры импульса:

когда первое изображение дало хорошие, регулярные импульсы.

Мой вопрос в том, является ли точность функций cos и sin причиной того, что вывод обратного 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

Спасибо еще раз!