Fast Fourier Transform(FFT C++)
1 #include <iostream> 2 #include <stdio.h> 3 #include <stdlib.h> 4 #include <cmath> 5 6 #define pi 3.1415926 7 #define Nall 32768 8 #define M 32 9 10 class CComplex 11 { 12 private: 13 float real, imag; 14 int BinaryLength; 15 public: 16 17 CComplex(float r, float i){real=r; imag=i;} 18 CComplex(){real=0; imag=0;} 19 CComplex(float r){real=r; imag=0;} 20 21 float GetReal(){ return real;} 22 float GetImag(){ return imag;} 23 24 void SetValue(){real=0; imag=0;} 25 void SetValue(float r, float i){real=r; imag=i;} 26 27 void Print() 28 { 29 printf("%f", real); 30 if(imag==0) 31 {printf("\n");} 32 else 33 { 34 if (imag>0) 35 {printf("+%fi\n",imag);} 36 else 37 {printf("%fi\n", imag);} 38 } 39 40 } 41 42 operator float(){ return this->real;} 43 44 CComplex operator =(CComplex a) 45 { 46 real=a.real; 47 imag=a.imag; 48 return *this; 49 } 50 51 52 template <class T> operator T() 53 { 54 return real; 55 } 56 57 58 CComplex operator +(CComplex c) 59 { 60 CComplex temp; 61 temp.real =real + c.real; 62 temp.imag =imag + c.imag; 63 return temp; 64 } 65 66 67 template <class T> CComplex operator +(T c) 68 { 69 CComplex temp; 70 temp.real =real + c; 71 temp.imag =imag ; 72 return temp; 73 } 74 75 template <class T> friend CComplex operator +(T a, CComplex b) 76 { 77 CComplex temp; 78 temp.real =a + b.real; 79 temp.imag =b.imag ; 80 return temp; 81 } 82 83 84 85 CComplex operator -(CComplex c) 86 { 87 CComplex temp; 88 temp.real =real - c.real; 89 temp.imag =imag - c.imag; 90 return temp; 91 } 92 93 template <class T> CComplex operator -(T c) 94 { 95 CComplex temp; 96 temp.real =real - c; 97 temp.imag =imag; 98 return temp; 99 } 100 101 template <class T> friend CComplex operator -(T a, CComplex b) 102 { 103 CComplex temp; 104 temp.real =a - b.real; 105 temp.imag =-b.imag ; 106 return temp; 107 } 108 109 110 111 CComplex operator *(CComplex c) 112 { 113 CComplex temp; 114 temp.real=real*c.real-imag*c.imag; 115 temp.imag=real*c.imag+imag*c.real; 116 return temp; 117 } 118 119 template <class T> CComplex operator *(T c) 120 { 121 CComplex temp; 122 temp.real =real * c; 123 temp.imag =imag * c; 124 return temp; 125 } 126 127 template <class T> friend CComplex operator *(T a, CComplex b) 128 { 129 CComplex temp; 130 temp.real =a * b.real; 131 temp.imag =a * b.imag ; 132 return temp; 133 } 134 135 136 CComplex operator /(CComplex c) 137 { 138 CComplex temp; 139 if(c.real==0&&c.imag==0) 140 {printf("Error! Divided by zero!\n"); 141 exit(1);} 142 float A; 143 A=1/(c.real*c.real+c.imag*c.imag); 144 temp.real=(real*c.real+imag*c.imag)*A; 145 temp.imag=(c.real*imag-c.imag*real)*A; 146 return temp; 147 } 148 149 template <class T> CComplex operator /(T c) 150 { 151 CComplex temp; 152 if (c==0) 153 {printf("0ÀÌ ºÐžð°¡ µÇŸúœÀŽÏŽÙ!\n"); 154 exit(1); 155 } 156 157 temp.real =real/c; 158 temp.imag =imag/c; 159 return temp; 160 } 161 162 template <class T> friend CComplex operator /(T a, CComplex b) 163 { 164 CComplex temp; 165 if ((b.real == 0)&&(b.imag==0)) 166 {printf("0 + 0i ÀÌ ºÐžð°¡ µÇŸúœÀŽÏŽÙ!\n"); 167 exit(1); 168 } 169 170 double A; 171 A=1/(b.real*b.real+b.imag*b.imag); 172 temp.real=a*b.real*A; 173 temp.imag=-a*b.imag*A; 174 return temp; 175 176 } 177 178 CComplex Exponential(float phase) 179 { 180 CComplex temp; 181 temp.real=cos(phase); 182 temp.imag=sin(phase); 183 return temp; 184 } 185 186 187 int ReverseBinaryOrder(int n, int N) 188 { 189 int arr[M]; 190 int temp=0; 191 192 BinaryLength =floor(log(float(N))/log(float(2))+0.5); 193 194 for (int j= BinaryLength-1;j>=0;j--) 195 { 196 arr[j]=(n>>j&1); 197 198 } 199 for (int j=0;j<=BinaryLength-1;j++) 200 { temp+=arr[j]*(int)pow((double)2,BinaryLength-1-j);} 201 202 return temp; 203 } 204 205 206 void FFT(CComplex *Out, CComplex *In, int N) 207 { 208 int i, j, k; 209 int J, B, P, L; 210 211 for (i=0; i<N; i++) 212 { 213 Out[i]=In[ReverseBinaryOrder(i, N)]; 214 } 215 216 printf("N = %d, M=%d\n",N, BinaryLength); 217 218 for(L=1;L<=BinaryLength;L++) 219 { 220 B=(int) pow((double) 2, (L-1)); 221 for (J=0;J<B;J++) 222 { 223 P=J*(int) pow((double)2, (BinaryLength-L)); 224 for (k=J;k<N; k+=(int) pow((double)2, L)) 225 { 226 In[k]=Out[k]+Out[k+B]*Exponential(-2*P*pi/N); 227 In[k+B]=Out[k]-Out[k+B]*Exponential(-2*P*pi/N); 228 } 229 } 230 for(i=0;i<N;i++) 231 { 232 Out[i]=In[i]; 233 } 234 } 235 236 } 237 238 template <class T> void FFT(CComplex *Out, T *In, int N) 239 { 240 int i, j, k; 241 int J, B, P, L; 242 CComplex *Buffer = new CComplex[N]; 243 244 for (i=0; i<N; i++) 245 { 246 Out[i]=In[ReverseBinaryOrder(i, N)]; 247 } 248 249 for (i = 0; i<N; i++) 250 { 251 Buffer[i] = In[i]; 252 } 253 254 printf("N = %d, M=%d\n",N, BinaryLength); 255 256 for(L=1;L<=BinaryLength;L++) 257 { 258 B=(int) pow((double) 2, (L-1)); 259 for (J=0;J<B;J++) 260 { 261 P=J*(int) pow((double)2, (BinaryLength-L)); 262 for (k=J;k<N; k+=(int) pow((double)2, L)) 263 { 264 Buffer[k]=Out[k]+Out[k+B]*Exponential(-2*P*pi/N); 265 Buffer[k+B]=Out[k]-Out[k+B]*Exponential(-2*P*pi/N); 266 } 267 } 268 for(i=0;i<N;i++) 269 { 270 Out[i]=Buffer[i]; 271 } 272 } 273 274 delete [] Buffer; 275 } 276 277 278 void IFFT(CComplex *Out, CComplex *In, int N) 279 { 280 int i, j, k; 281 int J, B, P, L; 282 283 for (i=0; i<N; i++) 284 { 285 Out[i]=In[ReverseBinaryOrder(i, N)]; 286 } 287 288 printf("N = %d, M=%d\n",N, BinaryLength); 289 290 for(L=1;L<=BinaryLength;L++) 291 { 292 B=(int) pow((double) 2, (L-1)); 293 for (J=0;J<B;J++) 294 { 295 P=J*(int) pow((double)2, (BinaryLength-L)); 296 for (k=J;k<N; k+=(int) pow((double)2, L)) 297 { 298 In[k]=Out[k]+Out[k+B]*Exponential(2*P*pi/N); 299 In[k+B]=Out[k]-Out[k+B]*Exponential(2*P*pi/N); 300 } 301 } 302 for(i=0;i<N;i++) 303 { 304 Out[i]=In[i]; 305 } 306 } 307 308 for(i=0;i<N;i++) 309 { 310 Out[i]=Out[i]/N; 311 } 312 313 314 } 315 316 317 template <class T> void IFFT(CComplex *Out, T *In, int N) 318 { 319 int i, j, k; 320 int J, B, P, L; 321 CComplex *IBuffer = new CComplex[N]; 322 323 for (i=0; i<N; i++) 324 { 325 Out[i]=In[ReverseBinaryOrder(i, N)]; 326 } 327 328 for (i = 0; i<N; i++) 329 { 330 IBuffer[i] = In[i]; 331 } 332 333 printf("N = %d, M=%d\n",N, BinaryLength); 334 335 for(L=1;L<=BinaryLength;L++) 336 { 337 B=(int) pow((double) 2, (L-1)); 338 for (J=0;J<B;J++) 339 { 340 P=J*(int) pow((double)2, (BinaryLength-L)); 341 for (k=J;k<N; k+=(int) pow((double)2, L)) 342 { 343 IBuffer[k]=Out[k]+Out[k+B]*Exponential(2*P*pi/N); 344 IBuffer[k+B]=Out[k]-Out[k+B]*Exponential(2*P*pi/N); 345 } 346 } 347 for(i=0;i<N;i++) 348 { 349 Out[i]=IBuffer[i]; 350 } 351 } 352 for(i=0;i<N;i++) 353 { 354 Out[i]=Out[i]/N; 355 } 356 357 delete [] IBuffer; 358 359 } 360 361 362 };
1 // FFT_Test.cpp 2 #include "CComplex.h" 3 4 int main() 5 { 6 CComplex *DataIn; 7 CComplex *DataOut; 8 DataIn=new CComplex[Nall]; 9 DataOut=new CComplex[Nall]; 10 CComplex *DoFFT= new CComplex; 11 12 int Data[32768]; 13 float Dataf[32768]; 14 double Datad[32768]; 15 int i; 16 int N1 = 1024; 17 int N2 = 2048; 18 int N3 = 4096; 19 int N4 = 8192; 20 int N5 = 16384; 21 int N6 = 32768; 22 23 ///////////////////////////////////////////////////////////////// 24 ////// N = 1024, FFTC2C 25 //////////////////////////////////////////////////////////////// 26 for (i=0; i<N1; i++) 27 { 28 DataIn[i].SetValue(i, 0); 29 } 30 31 DoFFT->FFT(DataOut, DataIn, N1); 32 33 FILE *fidDataOut= fopen("DataOut.txt","wt"); 34 for (i = 0; i < N1; i++) 35 { 36 if(DataOut[i].GetImag()==0) 37 { 38 fprintf(fidDataOut, "%f\n", DataOut[i].GetReal()); 39 } 40 else 41 { 42 if (DataOut[i].GetImag()>0) 43 { fprintf(fidDataOut, "%f", DataOut[i].GetReal()); 44 fprintf(fidDataOut, "+%fi\n",DataOut[i].GetImag());} 45 else 46 { fprintf(fidDataOut, "%f", DataOut[i].GetReal()); 47 fprintf(fidDataOut, "%fi\n",DataOut[i].GetImag());} 48 } 49 } 50 fclose(fidDataOut); 51 52 ///////////////////////////////////////////////////////////////// 53 ////// N = 2048, FFTC2C 54 //////////////////////////////////////////////////////////////// 55 56 for (i=0; i<N2; i++) 57 { 58 DataIn[i].SetValue(i, 0); 59 } 60 61 DoFFT->FFT(DataOut, DataIn, N2); 62 63 FILE *fidDataOut2= fopen("DataOut2.txt","wt"); 64 for (i = 0; i < N2; i++) 65 { 66 if(DataOut[i].GetImag()==0) 67 { 68 fprintf(fidDataOut2, "%f\n", DataOut[i].GetReal()); 69 } 70 else 71 { 72 if (DataOut[i].GetImag()>0) 73 { fprintf(fidDataOut2, "%f", DataOut[i].GetReal()); 74 fprintf(fidDataOut2, "+%fi\n",DataOut[i].GetImag());} 75 else 76 { fprintf(fidDataOut2, "%f", DataOut[i].GetReal()); 77 fprintf(fidDataOut2, "%fi\n",DataOut[i].GetImag());} 78 } 79 } 80 fclose(fidDataOut2); 81 82 ///////////////////////////////////////////////////////////////// 83 ////// N = 4096, IFFTC2C 84 //////////////////////////////////////////////////////////////// 85 86 for (i=0; i<N3; i++) 87 { 88 DataIn[i].SetValue(i, 0); 89 } 90 91 DoFFT->IFFT(DataOut, DataIn, N3); 92 93 FILE *fidDataOut3= fopen("DataOut3.txt","wt"); 94 for (i = 0; i < N3; i++) 95 { 96 if(DataOut[i].GetImag()==0) 97 { 98 fprintf(fidDataOut3, "%f\n", DataOut[i].GetReal()); 99 } 100 else 101 { 102 if (DataOut[i].GetImag()>0) 103 { fprintf(fidDataOut3, "%f", DataOut[i].GetReal()); 104 fprintf(fidDataOut3, "+%fi\n",DataOut[i].GetImag());} 105 else 106 { fprintf(fidDataOut3, "%f", DataOut[i].GetReal()); 107 fprintf(fidDataOut3, "%fi\n",DataOut[i].GetImag());} 108 } 109 } 110 fclose(fidDataOut3); 111 112 ///////////////////////////////////////////////////////////////// 113 ////// N = 8192, FFT Integer2C 114 //////////////////////////////////////////////////////////////// 115 116 for (i=0; i<N4; i++) 117 { 118 Data[i] = i; 119 } 120 121 DoFFT->FFT(DataOut, Data, N4); 122 123 FILE *fidDataOut4= fopen("DataOut4.txt","wt"); 124 for (i = 0; i < N4; i++) 125 { 126 if(DataOut[i].GetImag()==0) 127 { 128 fprintf(fidDataOut4, "%f\n", DataOut[i].GetReal()); 129 } 130 else 131 { 132 if (DataOut[i].GetImag()>0) 133 { fprintf(fidDataOut4, "%f", DataOut[i].GetReal()); 134 fprintf(fidDataOut4, "+%fi\n",DataOut[i].GetImag());} 135 else 136 { fprintf(fidDataOut4, "%f", DataOut[i].GetReal()); 137 fprintf(fidDataOut4, "%fi\n",DataOut[i].GetImag());} 138 } 139 } 140 fclose(fidDataOut4); 141 142 ///////////////////////////////////////////////////////////////// 143 ////// N = 16384, IFFT f2C 144 //////////////////////////////////////////////////////////////// 145 146 for (i=0; i<N5; i++) 147 { 148 Dataf[i] = i; 149 } 150 151 DoFFT->IFFT(DataOut, Dataf, N5); 152 153 FILE *fidDataOut5= fopen("DataOut5.txt","wt"); 154 for (i = 0; i < N5; i++) 155 { 156 if(DataOut[i].GetImag()==0) 157 { 158 fprintf(fidDataOut5, "%f\n", DataOut[i].GetReal()); 159 } 160 else 161 { 162 if (DataOut[i].GetImag()>0) 163 { fprintf(fidDataOut5, "%f", DataOut[i].GetReal()); 164 fprintf(fidDataOut5, "+%fi\n",DataOut[i].GetImag());} 165 else 166 { fprintf(fidDataOut5, "%f", DataOut[i].GetReal()); 167 fprintf(fidDataOut5, "%fi\n",DataOut[i].GetImag());} 168 } 169 } 170 fclose(fidDataOut5); 171 172 ///////////////////////////////////////////////////////////////// 173 ////// N = 32768, IFFT Doule2C 174 //////////////////////////////////////////////////////////////// 175 176 177 for (i=0; i<N6; i++) 178 { 179 Data[i] = i; 180 } 181 182 DoFFT->IFFT(DataOut, Data, N6); 183 184 FILE *fidDataOut6= fopen("DataOut6.txt","wt"); 185 for (i = 0; i < N6; i++) 186 { 187 if(DataOut[i].GetImag()==0) 188 { 189 fprintf(fidDataOut6, "%f\n", DataOut[i].GetReal()); 190 } 191 else 192 { 193 if (DataOut[i].GetImag()>0) 194 { fprintf(fidDataOut6, "%f", DataOut[i].GetReal()); 195 fprintf(fidDataOut6, "+%fi\n",DataOut[i].GetImag());} 196 else 197 { fprintf(fidDataOut6, "%f", DataOut[i].GetReal()); 198 fprintf(fidDataOut6, "%fi\n",DataOut[i].GetImag());} 199 } 200 } 201 fclose(fidDataOut6); 202 203 204 delete(DoFFT); 205 delete(DataIn); 206 delete(DataOut); 207 208 CComplex c1(2.3, 4.6), c2(3.6, 2.8), c3; 209 210 int a1 = 1; 211 float a2 = 2.5; 212 double a3 = 3.5; 213 214 215 printf("***************c1 c2****************************\n"); 216 c1.Print(); 217 c2.Print(); 218 219 printf("***************c1+-*/c2****************************\n"); 220 221 c3= c1+ c2; 222 c3.Print(); 223 224 c3 = c2 - c1; 225 c3.Print(); 226 227 c3 = c2*c1; 228 c3.Print(); 229 230 c3 = c1/c2; 231 c3.Print(); 232 233 printf("****************** c3 = a1 a2 a3 *************************\n"); 234 c3= c1; 235 c3.Print(); 236 237 c3 = a1; 238 c3.Print(); 239 240 c3 = a2; 241 c3.Print(); 242 243 c3 = a3; 244 c3.Print(); 245 246 printf("***********c3 = c1 + a1 a2 a3 ********************************\n"); 247 248 c3 = c1+a1; 249 c3.Print(); 250 251 c3 = a1+c1; 252 c3.Print(); 253 254 c3 = c1+a2; 255 c3.Print(); 256 257 c3 = a2+c1; 258 c3.Print(); 259 260 261 c3 = c1+a3; 262 c3.Print(); 263 264 c3 = a3 + c1; 265 c3.Print(); 266 267 printf("************c3 = c1 - a1 a2 a3*******************************\n"); 268 269 c3 = c1-a1; 270 c3.Print(); 271 272 c3 = a1-c1; 273 c3.Print(); 274 275 c3 = c1-a2; 276 c3.Print(); 277 278 c3 = a2-c1; 279 c3.Print(); 280 281 c3 = c1-a3; 282 c3.Print(); 283 284 c3 = a3-c1; 285 c3.Print(); 286 287 288 printf("**********c3 = c1 * a1 a2 a3*********************************\n"); 289 c3 = c1*a1; 290 c3.Print(); 291 292 c3 = a1*c1; 293 c3.Print(); 294 295 c3 = c1*a2; 296 c3.Print(); 297 298 c3 = a2*c1; 299 c3.Print(); 300 301 c3 = c1*a3; 302 c3.Print(); 303 304 c3 = a3*c1; 305 c3.Print(); 306 307 printf("***********c3 = c1/ a1 a2 a3********************************\n"); 308 c3 = c1/a1; 309 c3.Print(); 310 311 c3 = a1/c1; 312 c3.Print(); 313 314 c3 = c1/a2; 315 c3.Print(); 316 317 c3 = a2/c1; 318 c3.Print(); 319 320 c3 = c1/a3; 321 c3.Print(); 322 323 c3 = a3/c1; 324 c3.Print(); 325 326 // getchar(); 327 return 0; 328 329 }
./CComplex
N = 1024, M=10
N = 2048, M=11
N = 4096, M=12
N = 8192, M=13
N = 16384, M=14
N = 32768, M=15
***************c1 c2****************************
2.300000+4.600000i
3.600000+2.800000i
***************c1+-*/c2****************************
5.900000+7.400000i
1.300000-1.800000i
-4.599999+23.000000i
1.017308+0.486538i
****************** c3 = a1 a2 a3 *************************
2.300000+4.600000i
1.000000
2.500000
3.500000
***********c3 = c1 + a1 a2 a3 ********************************
3.300000+4.600000i
3.300000+4.600000i
4.800000+4.600000i
4.800000+4.600000i
5.800000+4.600000i
5.800000+4.600000i
************c3 = c1 - a1 a2 a3*******************************
1.300000+4.600000i
-1.300000-4.600000i
-0.200000+4.600000i
0.200000-4.600000i
-1.200000+4.600000i
1.200000-4.600000i
**********c3 = c1 * a1 a2 a3*********************************
2.300000+4.600000i
2.300000+4.600000i
5.750000+11.500000i
5.750000+11.500000i
8.050000+16.100000i
8.050000+16.100000i
***********c3 = c1/ a1 a2 a3********************************
2.300000+4.600000i
0.086957-0.173913i
0.920000+1.840000i
0.217391-0.434783i
0.657143+1.314286i
0.304348-0.608696i