FFT and IFFT in C++

Question

I'm using C++/C to perform forwards and reverse FFT on some data which is supposed to be the pulsed output of a laser.

The idea is to take the output, use a forward FFT to convert to the frequency domain, apply a linear best fit to the phase ( first unwrapping it) and then subtracting this best fit from the phase information.

The resulting phase and amplitude are then converted back to the time domain, with the ultimate aim being the compression of the pulses through phase compensation.

I've attempted to do this in MATLAB unsuccesfully, and have turned to C++ as a result. The forwards FFT is working fine, I took the basic recipe from Numerical recipes in C++, and used a function to modify it for complex inputs as following:

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

with the function ftt_basic() taken from 'Numerical recipes C++'.

My issue is that the form of input seems to affect the output of the Reverse FFT. This could be a precision issue, but I've looked around and it doesn't seem to have affected anyone else before.

Feeding the output of the forwards FFT directly back into the reverse FFT yields pulses identical to the intput:

However taking the power output taken asreal^2+imag^2 of the forwards FFT and copying it to an array such that:

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

and then using this as the input for the reverse FFT yields the following:

And finally, taking the output of the forwards FFT and copying such that:

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

where the Amplitude[i]=(real^2+imag^2)^0.5 and phase[i]=atan(imag/real). Yields the following power output when converting back to the time domain:

with a closer look at the pulse structure:

when the first picture yielded nice, regular pulses.

My question is, is it the precision of the cos and sin functions which cause the output of the reverse fft to become like this? Why is it that there is such a massive a difference between the different methods of inputting the complex data, and why is it that only when the data is directly fed back into the reverse FFT that the data in the time domain is identical to the original input into the forwads FFT?

Thank you.

*Edit here is the implemenntation of the functions :

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

Thanks again!

A

what is phase[i] ? a float or a double? – Examen

out of curiosity, what is it that you can do in c++ you couldn't do in MATLAB? – Lavation

phase[i] is an array of doubles. – A

I tried using the fft in matlab, but applying the phase changes didnt seem to affect the output pulses at all, except for shifting the actual pulse envelope rather than compressing it, theres another question a I posted with the problems I was having in matlab. So I tried doing it in C++ instead+ – A

#11925781 Here is the question concerning MATLAB, if you care to have a look! Thank you!! – A

@KRS-fun, if you add (at)username as seen in this comment, to a comment response, username will be notified that you have responded and will have a higher chance of revisiting the question. – Lavation

On the last graphics, x-axis step size is about 10e-9 but y-axis values are about 60. If its about the same calculation involving sin and cos functions and they are using doubles, precision shouldnt be a problem – Examen

The order of floating point calculations changes the error – Examen

This looks bigger than just numerical instability issues to me - when you write phase[i]=atan(imag/real), don't you want to use atan2(imag, real)? – Dependent

@Dependent I did indeed, thanks for spotting that! I've used atan2(imag,real) and now the output yields the same as when I use the format complex(real,imag) and complex(Amplitudecos,Amplitudesin) as the input for the reverse FFT. However I'm still struggling to understand why there seems to be a two pulse-periodicity, and why the peak amplitude is much smaller than for the first graph? – A

I think you are assuming you can do a transformation that isn't actually valid. I'm worried about the line "taking the power output taken as real^2+imag^2" - what is the actual transformation you're doing from the output of the FFT to the input of the IFFT? – Dependent

@Dependent The line you are referring to is the definition of power from a wave as the complex waveform multiplied by its complex conjugate. I think this is correct, at least for EM waves. The actual transformation to the output of the FFT is a subtraction of the line of best fit to the unwrapped phase, or removal of wavelength dependence of the phase, which should in theory compress the input pulses and increase peak powers. – A

I agree that this is mathematically correct, I'm wondering about the function signatures and definition of return values of your fft and reverse_fft functions. Does fft return the power or just the coefficients? does the reverse_fft want to take the power or the raw coefficients as input? If you could show the full set of code that calls fft(), transforms, and calls reverse_fft(), I think it's more likely we'll be able to spot the problem. – Dependent

@Dependent I've updated my post with the code, thanks again! – A

When you are loading the Data array from DataIn, you're storing the first complex values at Data[1] and Data[2]. Data[0] is left undefined yet you are passing Data to the fft_basic function. Is it expecting to find the data at the second element? If not, your input data is skewed. – Arium

fftw.org ... just sayin' ... – Pricefixing

Answer 1

Several likely errors:

   Data[0] == 0.0;

Probably it should be =, not ==?

fft_basic(Data, fftSize, InverseTransform);

I have no access to the source of this function; does it really expect the layout you're providing with real part at odd places and imaginary at even ones?

    //Swap the fft halfes

As I've told you in the other question, if you do swap them, you need to swap them back before the inverse transform. Do you perform this swap?

Does your data match the expectations of fft_basic function? Probably it expects that fftSize is a power of two.

Do you normalize the FFT results? The discrete FFT transform requires normalization.

Pathognomy