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 }