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 }
