FastFourierTransform (FFT)

代码来自互联网开源，仅作为收集整理使用。

FastFourierTransform.h

 1 #pragma once
 2 #include <stdio.h>
 3 #include <math.h>
 4 
 5 #ifndef INCLUDE_FASTTOURIERTRANSFORM
 6 #define INCLUDE_FASTTOURIERTRANSFORM
 7 
 8 /************************************************************************/
 9 /* CFastFourierTransform                                                */
10 /************************************************************************/
11 #define PI_2 6.283185F
12 #define PI   3.1415925F
13 class CFastFourierTransform
14 {
15 private:
16     float* xre;
17     float* xim;
18     float* mag;
19     float* fftSin;
20     float* fftCos;
21     int* fftBr;
22     int ss;
23     int ss2;
24     int nu;
25     int nu1;
26 
27     int BitRev(int j, int nu);
28     void PrepareFFTTables();
29 public:
30     CFastFourierTransform(int pSampleSize);
31     ~CFastFourierTransform(void);
32 
33     float* Calculate(float* pSample, size_t pSampleSize);
34 };
35 
36 #endif

FastFourierTransform.cpp

  1 #include "FastFourierTransform.h"
  2 
  3 /************************************************************************/
  4 /* CFastFourierTransform                                                */
  5 /************************************************************************/
  6 CFastFourierTransform::CFastFourierTransform(int pSampleSize)
  7 {
  8     xre = NULL;
  9     xim = NULL;
 10     mag = NULL;
 11     fftSin = NULL;
 12     fftCos = NULL;
 13     fftBr = NULL;
 14 
 15     ss = pSampleSize;
 16     ss2 = ss >> 1;
 17     nu = (int) (log((float)ss) / log((float)2));
 18     nu1 = nu - 1;
 19 
 20     xre = new float[ss]; // real part
 21     xim = new float[ss]; // image part
 22     mag = new float[ss2];
 23 
 24     PrepareFFTTables();
 25 }
 26 
 27 CFastFourierTransform::~CFastFourierTransform(void)
 28 {
 29     if(xre != NULL)
 30         delete [] xre;
 31 
 32     if(xim != NULL)
 33         delete [] xim;
 34 
 35     if(mag != NULL)
 36         delete [] mag;
 37 
 38     if(fftSin != NULL)
 39         delete [] fftSin;
 40 
 41     if(fftCos != NULL)
 42         delete [] fftCos;
 43 
 44     if(fftBr != NULL)
 45         delete [] fftBr;
 46 
 47     xre = NULL;
 48     xim = NULL;
 49     mag = NULL;
 50     fftSin = NULL;
 51     fftCos = NULL;
 52     fftBr = NULL;
 53 }
 54 
 55 void CFastFourierTransform::PrepareFFTTables()
 56 {
 57     int n2 = ss2;
 58     int nu1 = nu - 1;
 59 
 60     fftSin = new float[nu * n2];
 61     fftCos = new float[nu * n2];
 62 
 63     int k = 0;
 64     int x = 0;
 65     for (int l = 1; l <= nu; l++) {
 66         while (k < ss) {
 67             for (int i = 1; i <= n2; i++) {
 68                 float p = (float)BitRev(k >> nu1, nu);
 69                 float arg = (PI_2 * p) / (float) ss;
 70                 fftSin[x] = (float) sin(arg);
 71                 fftCos[x] = (float) cos(arg);
 72                 k++;
 73                 x++;
 74             }
 75 
 76             k += n2;
 77         }
 78 
 79         k = 0;
 80         nu1--;
 81         n2 >>= 1;
 82     }
 83 
 84     fftBr = new int[ss];
 85     for (k = 0; k < ss; k++)
 86         fftBr[k] = BitRev(k, nu);
 87 }
 88 
 89 int CFastFourierTransform::BitRev(int j, int nu) {
 90     int j1 = j;
 91     int k = 0;
 92     for (int i = 1; i <= nu; i++) {
 93         int j2 = j1 >> 1;
 94         k = ((k << 1) + j1) - (j2 << 1);
 95         j1 = j2;
 96     }
 97 
 98     return k;
 99 }
100 
101 float* CFastFourierTransform::Calculate(float* pSample, size_t pSampleSize) {
102     int n2 = ss2;
103     int nu1 = nu - 1;
104     int wAps = pSampleSize / ss;
105     size_t a = 0;
106 
107     for (size_t b = 0; a < pSampleSize; b++) {
108         xre[b] = pSample[a];
109         xim[b] = 0.0F;
110         a += wAps;
111     }
112 
113     int x = 0;
114     for (int l = 1; l <= nu; l++) {
115         for (int k = 0; k < ss; k += n2) {
116             for (int i = 1; i <= n2; i++) {
117                 float c = fftCos[x];
118                 float s = fftSin[x];
119                 int kn2 = k + n2;
120                 float tr = xre[kn2] * c + xim[kn2] * s;
121                 float ti = xim[kn2] * c - xre[kn2] * s;
122                 xre[kn2] = xre[k] - tr;
123                 xim[kn2] = xim[k] - ti;
124                 xre[k] += tr;
125                 xim[k] += ti;
126                 k++;
127                 x++;
128             }
129         }
130 
131         nu1--;
132         n2 >>= 1;
133     }
134 
135     for (int k = 0; k < ss; k++) {
136         int r = fftBr[k];
137         if (r > k) {
138             float tr = xre[k];
139             float ti = xim[k];
140             xre[k] = xre[r];
141             xim[k] = xim[r];
142             xre[r] = tr;
143             xim[r] = ti;
144         }
145     }
146 
147     mag[0] = (float) sqrt(xre[0] * xre[0] + xim[0] * xim[0]) / (float) ss;
148     for (int i = 1; i < ss2; i++)
149         mag[i] = (2.0F * (float) sqrt(xre[i] * xre[i] + xim[i] * xim[i])) / (float) ss;
150 
151     return mag;
152 }

 

posted @ 2016-04-07 02:20  Akatsuki-  阅读(594)  评论(0)  编辑  收藏  举报