title: FFT and IFFT in C++
body: <p>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.</p>
  
  <p>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.</p>
  
  <p>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.</p>
  
  <p>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:</p>
  
  <pre><code>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;
  }
  </code></pre>
  
  <p>with the function <code>ftt_basic()</code> taken from 'Numerical recipes C++'.</p>
  
  <p>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.</p>
  
  <p>Feeding the output of the forwards FFT directly back into the reverse FFT yields pulses identical to the intput:</p>
  
  <p><img src="https://i.stack.imgur.com/T2RMX.png" alt="enter image description here"></p>
  
  <p>However taking the power output taken as<code>real^2+imag^2</code> of the forwards FFT and copying it to an array such that:</p>
  
  <pre><code>Reverse_fft_input[i]=complex(real(forwardsoutput[i]),imag(forwardsoutput[i]));
  </code></pre>
  
  <p>and then using this as the input for the reverse FFT yields the following:
  <img src="https://i.stack.imgur.com/cZBVl.png" alt="enter image description here"></p>
  
  <p>And finally, taking the output of the forwards FFT and copying such that:</p>
  
  <pre><code>Reverse_fft_input[i]=complex( Amplitude[i]*cos(phase[i]), Amplitude[i]*sin(phase[i]));
  </code></pre>
  
  <p>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:</p>
  
  <p><img src="https://i.stack.imgur.com/ezRt0.png" alt="enter image description here"></p>
  
  <p>with a closer look at the pulse structure:</p>
  
  <p><img src="https://i.stack.imgur.com/7fAV2.png" alt="enter image description here"></p>
  
  <p>when the first picture yielded nice, regular pulses. </p>
  
  <p>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?</p>
  
  <p>Thank you.</p>
  
  <p>*Edit here is the implemenntation of the functions :</p>
  
  <pre><code>void TTWLM::SpectralAnalysis()
  </code></pre>
  
  <p>{</p>
  
  <pre><code>    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&lt;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&lt;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
  </code></pre>
  
  <p>Thanks again!</p>