Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
fftprint.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include <stdlib.h>
00008 #include <stdio.h>
00009 #include <string.h>
00010 #include <math.h>
00011
00012 #undef PI
00013 #define PI (3.141592653589793238462643)
00014
00015 typedef struct complex { float r , i ; } complex ;
00016
00017 static void csfft_print( int mode , int idim , complex * xc ) ;
00018 static void csfft_trigconsts( int idim ) ;
00019
00020 int main( int argc , char * argv[] )
00021 {
00022 int len , ii ;
00023 complex * cx ;
00024
00025 if( argc < 2 ){printf("Usage: fftprint len\n");exit(0);}
00026
00027 len = strtol( argv[1] , NULL , 10 ) ;
00028 if( len < 4 ){ fprintf(stderr,"Illegal length\n"); exit(1); }
00029
00030 ii = 2 ; do{ ii *= 2 ; } while( ii < len ) ;
00031 if( ii != len ){ fprintf(stderr,"Illegal length\n"); exit(1); }
00032
00033 cx = (complex *) malloc( sizeof(complex) * len ) ;
00034 csfft_print(-1,len,cx) ;
00035 exit(0) ;
00036 }
00037
00038 static complex * csplus = NULL , * csminus = NULL ;
00039 static int nold = -666 ;
00040
00041 static void csfft_trigconsts( int idim )
00042 {
00043 register unsigned int m, n, i0, i1, i2, i3, k;
00044 register complex *r0, *csp;
00045 register float co, si, f0, f1, f2, f3, f4;
00046 double al;
00047
00048 if( idim == nold ) return ;
00049
00050 if( idim > nold ){
00051 if( csplus != NULL ){ free(csplus) ; free(csminus) ; }
00052 csplus = (complex *) malloc( sizeof(complex) * idim ) ;
00053 csminus = (complex *) malloc( sizeof(complex) * idim ) ;
00054 if( csplus == NULL || csminus == NULL ){
00055 fprintf(stderr,"\n*** csfft cannot malloc space! ***\n"); exit(1) ;
00056 }
00057 }
00058 nold = n = idim ;
00059
00060 f1 = 1.0 ;
00061 m = 1; k = 0;
00062 while (n > m) {
00063 i3 = m << 1;
00064 f2 = m; al = f1*PI/f2;
00065 co = cos(al); si = sin(al);
00066 (csplus + k)->r = 1.;
00067 (csplus + k)->i = 0.;
00068 for (i0=0; i0 < m; i0++) {
00069 k++;
00070 csp = csplus + k; r0 = csp - 1;
00071 csp->r = r0->r * co - r0->i * si;
00072 csp->i = r0->i * co + r0->r * si;
00073 }
00074 m = i3;
00075 }
00076
00077 f1 = -1.0 ;
00078 m = 1; k = 0;
00079 while (n > m) {
00080 i3 = m << 1;
00081 f2 = m; al = f1*PI/f2;
00082 co = cos(al); si = sin(al);
00083 (csminus + k)->r = 1.;
00084 (csminus + k)->i = 0.;
00085 for (i0=0; i0 < m; i0++) {
00086 k++;
00087 csp = csminus + k; r0 = csp - 1;
00088 csp->r = r0->r * co - r0->i * si;
00089 csp->i = r0->i * co + r0->r * si;
00090 }
00091 m = i3;
00092 }
00093 return ;
00094 }
00095
00096
00097
00098
00099
00100
00101
00102
00103 #define EPS 1.e-5
00104
00105 void csfft_print( int mode , int idim , complex * xc )
00106 {
00107 register unsigned int m, n, i0, i1, i2, i3, k;
00108 register complex *r0, *r1, *csp;
00109 register float co, si, f0, f1, f2, f3, f4;
00110 double al;
00111
00112 csfft_trigconsts( idim ) ;
00113 csp = csplus ;
00114
00115 n = idim;
00116 i2 = idim >> 1;
00117 i1 = 0;
00118
00119
00120
00121 printf( "/**************************************/\n"
00122 "/* FFT routine unrolled of length %3d */\n"
00123 "\n"
00124 "void fft%d( int mode , complex * xc )\n"
00125 "{\n"
00126 " register complex * csp , * xcx=xc;\n"
00127 " register float f1,f2,f3,f4 ;\n"
00128 "\n"
00129 " /** perhaps initialize **/\n"
00130 "\n"
00131 " if( nold != %d ) csfft_trigconsts( %d ) ;\n"
00132 "\n"
00133 " csp = (mode > 0) ? csplus : csminus ; /* choose const array */\n"
00134 "\n"
00135 " /** data swapping part **/\n" ,
00136 idim,idim,idim,idim ) ;
00137
00138
00139
00140 for (i0=0; i0 < n; i0 ++) {
00141 if ( i1 > i0 ) {
00142
00143 printf("\n") ;
00144 printf(" f1 = xcx[%d].r ; f2 = xcx[%d].i ;\n",i0,i0) ;
00145 printf(" xcx[%d].r = xcx[%d].r ; xcx[%d].i = xcx[%d].i ;\n" , i0,i1,i0,i1) ;
00146 printf(" xcx[%d].r = f1 ; xcx[%d].i = f2 ;\n" , i1,i1) ;
00147
00148 }
00149 m = i2;
00150 while ( m && !(i1 < m) ) {
00151 i1 -= m;
00152 m >>= 1;
00153 }
00154 i1 += m;
00155 }
00156
00157
00158
00159 printf("\n /** butterflying part **/\n") ;
00160
00161 m = 1;
00162 k = 0;
00163 while (n > m) {
00164 i3 = m << 1;
00165 for (i0=0; i0 < m; i0 ++) {
00166 for (i1=i0; i1 < n; i1 += i3) {
00167
00168 printf("\n") ;
00169 if( csp[k].r == 1.0 ){
00170 printf(" f1 = xcx[%d].r ; f3 = xcx[%d].i ; /* cos=1 sin=0 */\n", i1+m,i1+m ) ;
00171
00172 } else if( fabs(csp[k].r) < EPS ){
00173 printf(" f1 = - xcx[%d].i * csp[%d].i ; /* cos=0 twiddles */\n", i1+m,k ) ;
00174 printf(" f3 = xcx[%d].r * csp[%d].i ;\n", i1+m,k ) ;
00175
00176 } else {
00177 printf(" f1 = xcx[%d].r * csp[%d].r - xcx[%d].i * csp[%d].i ; /* twiddles */\n",
00178 i1+m,k , i1+m,k ) ;
00179 printf(" f3 = xcx[%d].r * csp[%d].i + xcx[%d].i * csp[%d].r ;\n",
00180 i1+m,k , i1+m,k ) ;
00181 }
00182
00183 printf(" f2 = xcx[%d].r ; f4 = xcx[%d].i ;\n" , i1,i1 ) ;
00184 printf(" xcx[%d].r = f2-f1 ; xcx[%d].i = f4-f3 ;\n",i1+m,i1+m) ;
00185 printf(" xcx[%d].r = f2+f1 ; xcx[%d].i = f4+f3 ;\n",i1,i1) ;
00186
00187 }
00188 k++;
00189 }
00190 m = i3;
00191 }
00192
00193 printf("\n return ;\n}\n") ;
00194
00195 return ;
00196 }