fft1.hpp

00001 // http://www.mit.edu/~emin/source_code/fft/index.html
00002 
00003 #define SIN_2PI_16 0.38268343236508978
00004 #define SIN_4PI_16 0.707106781186547460
00005 #define SIN_6PI_16 0.923879532511286740
00006 #define C_P_S_2PI_16 1.30656296487637660
00007 #define C_M_S_2PI_16 0.54119610014619690
00008 #define C_P_S_6PI_16 1.3065629648763766
00009 #define C_M_S_6PI_16 -0.54119610014619690
00010 
00011 void fft1(float input[16],float output[16] ) {
00012   float temp, out0, out1, out2, out3, out4, out5, out6, out7, out8;
00013   float out9,out10,out11,out12,out13,out14,out15;
00014 
00015   out0=input[0]+input[8]; /* output[0 through 7] is the data that we */
00016   out1=input[1]+input[9]; /* take the 8 point real FFT of. */
00017   out2=input[2]+input[10];
00018   out3=input[3]+input[11];
00019   out4=input[4]+input[12];
00020   out5=input[5]+input[13];
00021   out6=input[6]+input[14];
00022   out7=input[7]+input[15];
00023   out8=input[0]-input[8];   /* inputs 8,9,10,11 are */
00024   out9=input[1]-input[9];   /* the Real part of the */
00025   out10=input[2]-input[10]; /* 4 point Complex FFT inputs.*/
00026   out11=input[3]-input[11]; 
00027   out12=input[12]-input[4]; /* outputs 12,13,14,15 are */
00028   out13=input[13]-input[5]; /* the Imaginary pars of  */
00029   out14=input[14]-input[6]; /* the 4 point Complex FFT inputs.*/
00030   out15=input[15]-input[7];
00031   temp=(out13-out9)*(SIN_2PI_16); 
00032   out9=out9*(C_P_S_2PI_16)+temp; 
00033   out13=out13*(C_M_S_2PI_16)+temp;
00034   out14*=(SIN_4PI_16);
00035   out10*=(SIN_4PI_16);
00036   out14=out14-out10;
00037   out10=out14+out10+out10;
00038   temp=(out15-out11)*(SIN_6PI_16);
00039   out11=out11*(C_P_S_6PI_16)+temp;
00040   out15=out15*(C_M_S_6PI_16)+temp;
00041   out8+=out10;
00042   out10=out8-out10-out10;
00043   out12+=out14;
00044   out14=out12-out14-out14;
00045   out9+=out11;
00046   out11=out9-out11-out11;
00047   out13+=out15;
00048   out15=out13-out15-out15;
00049   output[1]=out8+out9;
00050   output[7]=out8-out9;
00051   output[9]=out12+out13;
00052   output[15]=out13-out12;  
00053   output[5]=out10+out15;    /* implicit multiplies by */
00054   output[13]=out14-out11;   /* a twiddle factor of -j */                            
00055   output[3]=out10-out15;    /* implicit multiplies by */
00056   output[11]=-out14-out11;  /* a twiddle factor of -j */    
00057   out0=out0+out4;
00058   out4=out0-out4-out4;
00059   out1=out1+out5;
00060   out5=out1-out5-out5;
00061   out2+=out6;
00062   out6=out2-out6-out6;
00063   out3+=out7;
00064   out7=out3-out7-out7;  
00065   output[0]=out0+out2;
00066   output[4]=out0-out2;
00067   out1+=out3;
00068   output[12]=out3+out3-out1;
00069   output[0]+=out1;          /* Real Part of X[0] */
00070   output[8]=output[0]-out1-out1; /*Real Part of X[4] */
00071   out5*=SIN_4PI_16;
00072   out7*=SIN_4PI_16;
00073   out5=out5-out7;
00074   out7=out5+out7+out7;
00075   output[14]=out6-out7;     /* Imag Part of X[5] */
00076   output[2]=out5+out4;      /* Real Part of X[7] */
00077   output[6]=out4-out5;      /*Real Part of X[5] */
00078   output[10]=-out7-out6; /* Imag Part of X[7] */
00079 }