Using FFT to convert to frequency domain, then IFFT back to time domain C++

Question

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<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++'. 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...)

Answer 1

A couple of things:

1. 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
2. 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.
3. 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.
4. Just adding to the choir: solve it in Matlab first and then port to C++

Answer 2

Here I am, reading your original PROBLEM in 2019 (7 yrs later). You start out with TIMEDOMAIN data TD[], a set of equally spaced sampled amplitudes of a realtime laser signal. Let's verify what your data looks like. Does it have dead space followed by pulses (0 to 2V ?). Does it resemble a clock of off-time then on-time pulses ? And, was the data sampled "fast enough" that you get 3 good samples of the on-time duration ? Let's say that each pulse gets 3 good samples of amplitude that are recorded in the TD[array]. Then you take the FFT of the timedomain, to produce FD[], the frequency domain array result. Each location in FD[0] has a purpose. FD[0] is a DC level, FD[1] shows presence of the LOWEST FINDABLE frequency, and FD[N/2] shows presence of the HIGHEST FINDABLE frequency :"the NYQUIST frequency". Each location in FD[] is A + jB real + imag, and the phase at each frequency bin is ARCTAN (B/A). You manually visit each FD array location, and modify the REAL part (with a best-fit algorithm) and modify the IMAG part (with a best-fit algorithm). then all you have done is changed the amplitude of each contributing sinewave, and because each A and B has been changed, the phase INFO at each FD[location] has been modified. Let's say that you correctly changed each FD[].real and FD[].imag for each array location from term [0] to term [N/2] and the mirror image terms from [N/2]+1 to the end of the array. When you do the Inverse FFT of that, you just made the timedomain signal of the modified frequency domain signal. Since you did not ZERO out any of the terms, you just "changed" the size of each contributing sine wave of the original timedomain signal. IF YOU REALLY sampled the original signal with sufficient bandwidth, then the result you get should resemble the original signal, with more rounded corners at each pulse, with lower amplitude pulses. It is not enough to say that you have a timedomain array of the input signal. You also need to ensure that there were more than 3 amplitude samples for each piece of off-time and more than 3 amplitude samples for each piece of on-time of the various laser pulses. Otherwise, you may have just recorded a slow, aliased version of the actual laser signal activity. If you never recorded sufficient samples of the laser signal, the FFT could not give you the desired FREQ information that you were after. - Cesar