fft.cpp

00001 /*
00002  * This is not my code and I do not take credit for it.
00003  * http://www.mit.edu/~emin/source_code/fft/index.html
00004  */
00005 
00006 #define SIN_2PI_16 0.38268343236508978
00007 #define SIN_4PI_16 0.707106781186547460
00008 #define SIN_6PI_16 0.923879532511286740
00009 #define C_P_S_2PI_16 1.30656296487637660
00010 #define C_M_S_2PI_16 0.54119610014619690
00011 #define C_P_S_6PI_16 1.3065629648763766
00012 #define C_M_S_6PI_16 -0.54119610014619690
00013 
00014 /* INPUT: float input[16], float output[16] */
00015 /* OUTPUT: none */
00016 /* EFFECTS:  Places the 16 point fft of input in output in a strange */
00017 /* order using 10 real multiplies and 79 real adds. */
00018 /* Re{F[0]}= out0 */
00019 /* Im{F[0]}= 0 */
00020 /* Re{F[1]}= out8 */
00021 /* Im{F[1]}= out12 */
00022 /* Re{F[2]}= out4 */
00023 /* Im{F[2]}= -out6 */
00024 /* Re{F[3]}= out11 */
00025 /* Im{F[3]}= -out15 */
00026 /* Re{F[4]}= out2 */
00027 /* Im{F[4]}= -out3 */
00028 /* Re{F[5]}= out10 */
00029 /* Im{F[5]}= out14 */
00030 /* Re{F[6]}= out5 */
00031 /* Im{F[6]}= -out7 */
00032 /* Re{F[7]}= out9 */
00033 /* Im{F[7]}= -out13 */
00034 /* Re{F[8]}= out1 */
00035 /* Im{F[8]}=0 */
00036 /* F[9] through F[15] can be found by using the formula */
00037 /* Re{F[n]}=Re{F[(16-n)mod16]} and Im{F[n]}= -Im{F[(16-n)mod16]} */
00038 
00039 /* Note using temporary variables to store intermediate computations */
00040 /* in the butterflies might speed things up.  When the current version */
00041 /* needs to compute a=a+b, and b=a-b, I do a=a+b followed by b=a-b-b.  */
00042 /* So practically everything is done in place, but the number of adds */
00043 /* can be reduced by doinc c=a+b followed by b=a-b. */
00044 
00045 /* The algorithm behind this program is to find F[2k] and F[4k+1] */
00046 /* seperately.  To find F[2k] we take the 8 point Real FFT of x[n]+x[n+8] */
00047 /* for n from 0 to 7.  To find F[4k+1] we take the 4 point Complex FFT of */
00048 /* exp(-2*pi*j*n/16)*{x[n] - x[n+8] + j(x[n+12]-x[n+4])} for n from 0 to 3.*/
00049 
00050 void fft(float input[16],float output[16] ) {
00051   float temp, out0, out1, out2, out3, out4, out5, out6, out7, out8;
00052   float out9,out10,out11,out12,out13,out14,out15;
00053 
00054   out0=input[0]+input[8]; /* output[0 through 7] is the data that we */
00055   out1=input[1]+input[9]; /* take the 8 point real FFT of. */
00056   out2=input[2]+input[10];
00057   out3=input[3]+input[11];
00058   out4=input[4]+input[12];
00059   out5=input[5]+input[13];
00060   out6=input[6]+input[14];
00061   out7=input[7]+input[15];
00062 
00063 
00064 
00065   out8=input[0]-input[8];   /* inputs 8,9,10,11 are */
00066   out9=input[1]-input[9];   /* the Real part of the */
00067   out10=input[2]-input[10]; /* 4 point Complex FFT inputs.*/
00068   out11=input[3]-input[11]; 
00069   out12=input[12]-input[4]; /* outputs 12,13,14,15 are */
00070   out13=input[13]-input[5]; /* the Imaginary pars of  */
00071   out14=input[14]-input[6]; /* the 4 point Complex FFT inputs.*/
00072   out15=input[15]-input[7];
00073 
00074   /*First we do the "twiddle factor" multiplies for the 4 point CFFT */
00075   /*Note that we use the following handy trick for doing a complex */
00076   /*multiply:  (e+jf)=(a+jb)*(c+jd) */
00077   /*   e=(a-b)*d + a*(c-d)   and    f=(a-b)*d + b*(c+d)  */
00078 
00079   /* C_M_S_2PI/16=cos(2pi/16)-sin(2pi/16) when replaced by macroexpansion */
00080   /* C_P_S_2PI/16=cos(2pi/16)+sin(2pi/16) when replaced by macroexpansion */
00081   /* (SIN_2PI_16)=sin(2pi/16) when replaced by macroexpansion */
00082   temp=(out13-out9)*(SIN_2PI_16); 
00083   out9=out9*(C_P_S_2PI_16)+temp; 
00084   out13=out13*(C_M_S_2PI_16)+temp;
00085   
00086   out14*=(SIN_4PI_16);
00087   out10*=(SIN_4PI_16);
00088   out14=out14-out10;
00089   out10=out14+out10+out10;
00090   
00091   temp=(out15-out11)*(SIN_6PI_16);
00092   out11=out11*(C_P_S_6PI_16)+temp;
00093   out15=out15*(C_M_S_6PI_16)+temp;
00094 
00095   /* The following are the first set of two point butterfiles */
00096   /* for the 4 point CFFT */
00097 
00098   out8+=out10;
00099   out10=out8-out10-out10;
00100 
00101   out12+=out14;
00102   out14=out12-out14-out14;
00103 
00104   out9+=out11;
00105   out11=out9-out11-out11;
00106 
00107   out13+=out15;
00108   out15=out13-out15-out15;
00109 
00110   /*The followin are the final set of two point butterflies */
00111   output[1]=out8+out9;
00112   output[7]=out8-out9;
00113 
00114   output[9]=out12+out13;
00115   output[15]=out13-out12;
00116   
00117   output[5]=out10+out15;        /* implicit multiplies by */
00118   output[13]=out14-out11;        /* a twiddle factor of -j */                            
00119   output[3]=out10-out15;  /* implicit multiplies by */
00120   output[11]=-out14-out11;  /* a twiddle factor of -j */
00121 
00122   
00123   /* What follows is the 8-point FFT of points output[0-7] */
00124   /* This 8-point FFT is basically a Decimation in Frequency FFT */
00125   /* where we take advantage of the fact that the initial data is real*/
00126 
00127   /* First set of 2-point butterflies */
00128     
00129   out0=out0+out4;
00130   out4=out0-out4-out4;
00131   out1=out1+out5;
00132   out5=out1-out5-out5;
00133   out2+=out6;
00134   out6=out2-out6-out6;
00135   out3+=out7;
00136   out7=out3-out7-out7;
00137 
00138   /* Computations to find X[0], X[4], X[6] */
00139   
00140   output[0]=out0+out2;
00141   output[4]=out0-out2;
00142   out1+=out3;
00143   output[12]=out3+out3-out1;
00144 
00145   output[0]+=out1;  /* Real Part of X[0] */
00146   output[8]=output[0]-out1-out1; /*Real Part of X[4] */
00147   /* out2 = Real Part of X[6] */
00148   /* out3 = Imag Part of X[6] */
00149   
00150   /* Computations to find X[5], X[7] */
00151 
00152   out5*=SIN_4PI_16;
00153   out7*=SIN_4PI_16;
00154   out5=out5-out7;
00155   out7=out5+out7+out7;
00156 
00157   output[14]=out6-out7; /* Imag Part of X[5] */
00158   output[2]=out5+out4; /* Real Part of X[7] */
00159   output[6]=out4-out5; /*Real Part of X[5] */
00160   output[10]=-out7-out6; /* Imag Part of X[7] */
00161 }