FFT и IFFT в С++

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

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

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

Я попытался сделать это в MATLAB неубедительно, и в результате обратился к С++. Форвард FFT работает нормально, я взял основной рецепт из Numericical 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(), взятой из "Численных рецептов С++".

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

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

enter image description here

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

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

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

И, наконец, взяв вывод форварда БПФ и скопировав его так:

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

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

enter image description here

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

enter image description here

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

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

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

Ответы

Ответ 1

Несколько возможных ошибок:

   Data[0] == 0.0;

Вероятно, это должно быть =, а не ==?

fft_basic (данные, fftSize, InverseTransform);

У меня нет доступа к источнику этой функции; действительно ли вы ожидаете, что макет, который вы предоставляете с реальной частью в нечетных местах, и мнимый на четных?

    //Swap the fft halfes

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

Соответствуют ли ваши данные ожиданиям функции fft_basic? Вероятно, он ожидает, что fftSize является мощностью двух.

Вы нормализуете результаты FFT? Для дискретного преобразования БПФ требуется нормализация.