Main Page | Namespaces | Classes | Compounds | Files | Compound Members | Related
t2.cpp
00001 #define M_PI 3.14159265359 00002 #define M_PI_F 3.14159265359f 00003 00004 #include <math.h> 00005 00006 #include "fft.h" 00007 00008 00009 void help() {;} 00010 void test() {;} 00011 00012 #define SWAP_F(x,y) {float tmp=x; x=y; y=tmp;} 00013 00014 void revbin_permute_f(float *fr, long n) 00015 { 00016 long m,j; 00017 for (m=1,j=0; m<n-1; m++) 00018 { 00019 long k; 00020 for(k=n>>1; (!((j^=k)&k)); k>>=1); 00021 00022 if (j>m) SWAP_F(fr[m], fr[j]); 00023 } 00024 } 00025 00026 void fft_f(float *fr, float *fi, int ldn, int is) 00027 /* 00028 * radix 2 fft a la cooley tukey 00029 * 00030 */ 00031 { 00032 int n2,ldm,m,mh,j,r; 00033 int t1, t2; 00034 float pi,phi,c,s; 00035 float ur,vr, ui,vi; 00036 00037 n2=1<<ldn; 00038 00039 pi=is*M_PI_F; 00040 00041 revbin_permute_f(fr,n2); 00042 revbin_permute_f(fi,n2); 00043 00044 for (ldm=1; ldm<=ldn; ++ldm) 00045 { 00046 m=(1<<ldm); /* m=2^ldm */ 00047 mh=(m>>1); /* mh=m/2 */ 00048 00049 phi=pi/(float)(mh); 00050 00051 for (j=0; j<mh; ++j) 00052 { 00053 float w=phi*(float)j; 00054 00055 c=(float) cos(w); 00056 s=(float) sin(w); 00057 00058 for (r=0; r<n2; r+=m) 00059 { 00060 /* 00061 * u=f[t1] 00062 * v=f[t2]*exp(+-2*pi*i*j/m) 00063 * f[t1]= (u+v) 00064 * f[t2]= (u-v) 00065 */ 00066 00067 t1=r+j; 00068 t2=t1+mh; 00069 00070 vr=fr[t2]*c-fi[t2]*s; 00071 vi=fr[t2]*s+fi[t2]*c; 00072 00073 ur=fr[t1]; 00074 fr[t1]+=vr; 00075 fr[t2]=ur-vr; 00076 00077 ui=fi[t1]; 00078 fi[t1]+=vi; 00079 fi[t2]=ui-vi; 00080 } 00081 } 00082 } 00083 }
Docs made by Doxygen. Email: Mikael Christensen