1
$\begingroup$

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

$\endgroup$
  • 8
    $\begingroup$ Any reason why you're trying to roll your own FFT code instead of using something proven like FFTW - or if it looks too big and complicated to you, FFTreal or kissFFT? Also, you'd better try to get your idea working in matlab first, where debugging is usually faster, visualization capabilities better, and where many signal processing primitives are available allowing very compact code... $\endgroup$ – pichenettes Aug 15 '12 at 12:44
  • $\begingroup$ Hi, thanks for the swift response! I'm using this code because it was included in the model I am using, Matlab was actually my first recourse but I ran into certain issues with that as well, I've already posted a question on this in stackoverflow: stackoverflow.com/questions/11925281/… $\endgroup$ – KRS-fun Aug 15 '12 at 12:52
  • 1
    $\begingroup$ It's not clear what your algorithm is and exactly what it tries to accomplish. As you found out, if you strip the phase information from the DFT outputs and try to IDFT the squared magnitude only, what you see in the time domain won't likely resemble the original signal. Your second attempt, where you first convert to amplitude/phase representation and then try to IDFT should work, but it's possible you made a mistake. Since that portion of your code isn't shown, it's hard to tell. As pichenettes recommended, this is a process best hammered out in a MATLAB/Octave-like environment. $\endgroup$ – Jason R Aug 15 '12 at 12:57
4
$\begingroup$

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++
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.