Matlab FFT and home brewed FFT

I'm trying to verify an FFT algorithm I should use for a project VS the same thing on Matlab. The point is that with my own C FFT function I always get the right (the second one) part of the double sided FFT spectrum evaluated in Matlab and not the first one as "expected".

For instance if my third bin is in the form a+i*b the third bin of Matlab's FFT is a-i*b. A and b values are the same but i always get the complex conjugate of Matlab's. I know that in terms of amplitudes and power there's no trouble (cause abs value) but I wonder if in terms of phases I'm going to read always wrong angles.

Im not so skilled in Matlab to know (and I have not found useful infos on the web) if Matlab FFT maybe returns the FFT spectre with negative frequencies first and then positive... or if I have to fix my FFT algorithm... or if it is all ok because phases are the unchanged regardless wich part of FFT we choose as single side spectrum (but i doubt about this last option).

Example:

If S is the sample array with N=512 samples, Y = fft(S) in Matlab return the FFT as (the sign of the imaginary part in the first half of the array are random, just to show the complex conjugate difference for the second part):

1   A1    +  i*B1     (DC, B1 is always zero)
2   A2    +  i*B2
3   A3    -  i*B3
4   A4    +  i*B4
5   A5    +  i*B5
...
253 A253  -  i*B253
254 A254  +  i*B254
255 A255  +  i*B255
256 A256  +  i*B256
257 A257  -  i*B257   (Nyquyst, B257 is always zero)
258 A256  -  i*B256
259 A255  -  i*B255
260 A254  -  i*B254
261 A253  +  i*B253
...
509 A5    -  i*B5
510 A4    -  i*B4
511 A3    +  i*B3
512 A2    -  i*B2

My FFT implementation returns only 256 values (and that's ok) in the the Y array as:

1   1    A1    +  i*B1     (A1 is the DC, B1 is Nyquist, both are pure Real numbers)
2   512  A2    -  i*B2
3   511  A3    +  i*B3
4   510  A4    -  i*B4
5   509  A5    +  i*B5
...
253 261  A253  +  i*B253
254 260  A254  -  i*B254
255 259  A255  -  i*B255
256 258  A256  -  i*B256

Where the first column is the proper index of my Y array and the second is just the reference of the relative row in the Matlab FFT implementation.

As you can see my FFT implementation (DC apart) returns the FFT like the second half of the Matlab's FFT (in reverse order).

To summarize: even if I use fftshift as suggested, it seems that my implementation always return what in the Matlab FFT should be considered the negative part of the spectrum. Where is the error???

This is the code I use:

Note 1: the FFT array is not declared here and it is changed inside the function. Initially it holds the N samples (real values) and at the end it contains the N/2 +1 bins of the single sided FFT spectrum.

Note 2: the N/2+1 bins are stored in N/2 elements only because the DC component is always real (and it is stored in FFT[0]) and also the Nyquyst (and it is stored in FFT[1]), this exception apart all the other even elements K holds a real number and the oven elements K+1 holds the imaginary part.

void Fft::FastFourierTransform( bool inverseFft ) {
   double twr, twi, twpr, twpi, twtemp, ttheta;
   int i, i1, i2, i3, i4, c1, c2;
   double h1r, h1i, h2r, h2i, wrs, wis;
   int nn, ii, jj, n, mmax, m, j, istep, isign;
   double wtemp, wr, wpr, wpi, wi;    
   double theta, tempr, tempi;

   // NS is the number of samples and it must be a power of two
   if( NS == 1 )
      return;

   if( !inverseFft ) {
      ttheta = 2.0 * PI / NS;
      c1 = 0.5;
      c2 = -0.5;
   }
   else {
      ttheta = 2.0 * PI / NS;
      c1 = 0.5;
      c2 = 0.5;
      ttheta = -ttheta;
      twpr = -2.0 * Pow( Sin( 0.5 * ttheta ), 2 );
      twpi = Sin(ttheta);
      twr = 1.0+twpr;
      twi = twpi;
      for( i = 2; i <= NS/4+1; i++ ) {
         i1 = i+i-2;
         i2 = i1+1;
         i3 = NS+1-i2;
         i4 = i3+1;
         wrs = twr;
         wis = twi;
         h1r = c1*(FFT[i1]+FFT[i3]);
         h1i = c1*(FFT[i2]-FFT[i4]);
         h2r = -c2*(FFT[i2]+FFT[i4]);
         h2i = c2*(FFT[i1]-FFT[i3]);
         FFT[i1] = h1r+wrs*h2r-wis*h2i;
         FFT[i2] = h1i+wrs*h2i+wis*h2r;
         FFT[i3] = h1r-wrs*h2r+wis*h2i;
         FFT[i4] = -h1i+wrs*h2i+wis*h2r;
         twtemp = twr;
         twr = twr*twpr-twi*twpi+twr;
         twi = twi*twpr+twtemp*twpi+twi;
      }
      h1r = FFT[0];
      FFT[0] = c1*(h1r+FFT[1]);
      FFT[1] = c1*(h1r-FFT[1]);
   }
   if( inverseFft )
      isign = -1;
   else
      isign = 1;
   n = NS;
   nn = NS/2;
   j = 1;
   for(ii = 1; ii <= nn; ii++) {
      i = 2*ii-1;
      if( j>i ) {
         tempr = FFT[j-1];
         tempi = FFT[j];
         FFT[j-1] = FFT[i-1];
         FFT[j] = FFT[i];
         FFT[i-1] = tempr;
         FFT[i] = tempi;
      }
      m = n/2;
      while( m>=2 && j>m ) {
         j = j-m;
         m = m/2;
      }
      j = j+m;
   }
   mmax = 2;
   while(n>mmax) {
      istep = 2*mmax;
      theta = 2.0 * PI /(isign*mmax);
      wpr = -2.0 * Pow( Sin( 0.5 * theta ), 2 );
      wpi = Sin(theta);
      wr = 1.0;
      wi = 0.0;
      for(ii = 1; ii <= mmax/2; ii++) {
         m = 2*ii-1;
         for(jj = 0; jj <= (n-m)/istep; jj++) {
            i = m+jj*istep;
            j = i+mmax;
            tempr = wr*FFT[j-1]-wi*FFT[j];
            tempi = wr*FFT[j]+wi*FFT[j-1];
            FFT[j-1] = FFT[i-1]-tempr;
            FFT[j] = FFT[i]-tempi;
            FFT[i-1] = FFT[i-1]+tempr;
            FFT[i] = FFT[i]+tempi;
         }
         wtemp = wr;
         wr = wr*wpr-wi*wpi+wr;
         wi = wi*wpr+wtemp*wpi+wi;
      }
      mmax = istep;
   }
   if( inverseFft )
      for(i = 1; i <= 2*nn; i++)
         FFT[i-1] = FFT[i-1]/nn;
   if( !inverseFft ) {
      twpr = -2.0 * Pow( Sin( 0.5 * ttheta ), 2 );
      twpi = Sin(ttheta);
      twr = 1.0+twpr;
      twi = twpi;
      for(i = 2; i <= NS/4+1; i++) {
         i1 = i+i-2;
         i2 = i1+1;
         i3 = NS+1-i2;
         i4 = i3+1;
         wrs = twr;
         wis = twi;
         h1r = c1*(FFT[i1]+FFT[i3]);
         h1i = c1*(FFT[i2]-FFT[i4]);
         h2r = -c2*(FFT[i2]+FFT[i4]);
         h2i = c2*(FFT[i1]-FFT[i3]);
         FFT[i1] = h1r+wrs*h2r-wis*h2i;
         FFT[i2] = h1i+wrs*h2i+wis*h2r;
         FFT[i3] = h1r-wrs*h2r+wis*h2i;
         FFT[i4] = -h1i+wrs*h2i+wis*h2r;
         twtemp = twr;
         twr = twr*twpr-twi*twpi+twr;
         twi = twi*twpr+twtemp*twpi+twi;
      }
      h1r = FFT[0];
      FFT[0] = h1r+FFT[1];  // DC
      FFT[1] = h1r-FFT[1];  // FS/2 (NYQUIST)
   }            
   return;
}

回答1:

In matlab try using fftshift(fft(...)). Matlab doesn't automatically shift the spectrum after the FFT is called which is why they implemented the fftshift() function.

回答2:

It is simply a matlab formatting thing. Basically, matlab arrange Fourier transform in following order

DC, (DC-1), .... (Nyquist-1), -Nyquist, -Nyquist+1, ..., DC-1

Let's say you have a 8 point sequence: [1 2 3 1 4 5 1 3]

In your signal processing class, your professor probably draws the Fourier spectrum based on a Cartesian system ( negative -> positive for x axis); So your DC should be located at 0 (the 4th position in your fft sequence, assuming position index here is 0-based) on your x axis.

In matlab, the DC is the very first element in the fft sequence, so you need to to fftshit() to swap the first half and second half of the fft sequence such that DC will be located at 4th position (position is 0-based indexed)

I am attaching a graph here so you may have a visual:

where a is the original 8-point sequence; FT(a) is the Fourier transform of a.

The matlab code is here:

a = [1 2 3 1 4 5 1 3];
A = fft(a);

N = length(a);
x = -N/2:N/2-1;

figure
subplot(3,1,1), stem(x, a,'o'); title('a'); xlabel('time')
subplot(3,1,2), stem(x, fftshift(abs(A),2),'o'); title('FT(a) in signal processing'); xlabel('frequency')
subplot(3,1,3), stem(x, abs(A),'o'); title('FT(a) in matlab'); xlabel('frequency')