Using FFT to convert to frequency domain, then IFFT back to time domain C++ I've posted the same question on stackoverflow.com with little success, so I thought I would try here!

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&lt;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&lt;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&lt;fftSize/2; i++)
                {
                        temp[i+fftSize/2] = DataOut[i];
                }
                for(int i=fftSize/2; i&lt;fftSize; i++)
                {
                        temp[i-fftSize/2] = DataOut[i];
                }
                for(int i=0; i&lt;fftSize; i++)
                {
                        DataOut[i] = temp[i];
                }
                delete[] temp;
        }
        delete[] Data;
}


with the function ftt_basic() taken from 'Numerical recipes C++'. The input variable InverseTransform is simply the direction of the FFT: forward or reverse.

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 regular coupled pulses, with the wrong amplitude and periodicity.

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 a really mess output in the time domain, with a peak towards the middle, resembling the power spectrum in the frequency domain.

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? All variables are stored as type 'double'. 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.

( I would have posted pictures but unfortunately I can't do that yet...)
A couple of things:


Your Data[] array seems to start at on offset of "1". That's unusual for C++, where "0" is the first index. Matlab starts at 1. This may be part of the calling convention for fft_basic, but it's a bit strange
Doing the inverse Fourier transform of the Power Spectrum (magnitude squared of the forward transform) gives you the auto correlation. That has nothing to do with phase compensation at all. Wrong algorithm.
You cannot calculate the phase as atan(imag/real) since you'll get half of the quadrants wrong. The correct functions is atan2(imag, real) which takes into account the relative signs of real and imaginary part.
Just adding to the choir: solve it in Matlab first and then port to C++