-
- class CComplex
- {
- public:
- short R;
- short I;
- CComplex() {}
- CComplex(short r, short i) { R = r; I = i; }
- void Set(short r, short i) { R = r; I = i; }
- CComplex(CComplex& complex) { R = complex.R; I = complex.I; }
- inline void operator=(CComplex& complex) { R = complex.R; I = complex.I; }
- void operator+=(CComplex& complex)
- {
- R += complex.R;
- I += complex.I;
- }
- void operator-=(CComplex& complex)
- {
- R -= complex.R;
- I -= complex.I;
- }
- inline void Add(CComplex& c1, CComplex& c2)
- {
- R = c1.R + c2.R;
- I = c1.I + c2.I;
- }
- inline void Sub(CComplex& c1, CComplex& c2)
- {
- R = c1.R - c2.R;
- I = c1.I - c2.I;
- }
- inline void Add(CComplex& c1, CComplex& c2, int shift)
- {
- R = (c1.R + c2.R) >> shift;
- I = (c1.I + c2.I) >> shift;
- }
- inline void Sub(CComplex& c1, CComplex& c2, int shift)
- {
- R = (c1.R - c2.R) >> shift;
- I = (c1.I - c2.I) >> shift;
- }
- inline void Mul(CComplex& c1, CComplex& c2, int shift)
- {
- R = (c1.R * c2.R + c1.I * c2.I) >> shift;
- I = (c1.I * c2.R - c1.R * c2.I) >> shift;
- }
- inline void MulConj(CComplex& c1, CComplex& c2, int shift)
- {
- R = (c1.R * c2.R - c1.I * c2.I) >> shift;
- I = (c1.I * c2.R + c1.R * c2.I) >> shift;
- }
- inline void Shift(int n) { R >>= n; I >>= n; }
- };
在嵌入式，构造函数中的 BitReversal, W 数组应该预先离线计算
- class CFFT
- {
- protected:
- static const int N = 32;
- CComplex W[N];
- int BitReversal[N];
- public:
- CComplex Y[N];
- protected:
- void ReverseBit()
- {
- int rev = 0;
- int t = N/2;
- int mask;
- for (int i = 0; i < N-1; i++) {
- BitReversal = rev;
- mask = t;
- while (rev >= mask) {
- rev -= mask;
- mask >>= 1;
- }
- rev += mask;
- }
- BitReversal[N-1] = N-1;
- }
- public:
- CFFT()
- {
- for (int i=0; i<N; i++) {
- W.R = (short) ( ::cos(PI * 2*i / N) * 16384 + 0.5);
- W.I = (short) (-::sin(PI * 2*i / N) * 16384 + 0.5);
- }
- ReverseBit();
- }
- void Calculate(short* x)
- {
- short t;
- int i, j, n, index1, index2, indexW, di;
- CComplex result;
- for(i=0; i<N; i++) {
- n = BitReversal;
- if(n > i) {
- t = x;
- x = x[n];
- x[n] = t;
- }
- }
- for (i=0; i<N; i+=2) {
- t = x[i+1];
- Y.Set(x + t, 0);
- Y[i+1].Set(x - t, 0);
- }
- for (di=N>>2, n=2; n<N; n<<=1, di >>= 1) {
- for (j=0; j<N; j+=(n<<1)) {
- indexW = 0;
- for (i=0; i<n; i++) {
- index1 = i + j;
- index2 = index1 + n;
- result.Mul(Y[index2], W[indexW], 14);
- Y[index2].Sub(Y[index1], result, 1);
- Y[index1].Add(Y[index1], result, 1);
- indexW += di;
- }
- }
- }
- }
- void Inverse()
- {
- int i, j, n, index1, index2, indexW, di;
- CComplex result;
- for(i=0; i<N; i++) {
- n = BitReversal;
- if(n > i) {
- result = Y;
- Y = Y[n];
- Y[n] = result;
- }
- }
- for (i=0; i<N; i+=2) {
- result = Y[i+1];
- Y[i+1].Sub(Y, result, 1);
- Y.Add(Y, result, 1);
- }
- for (di=N>>2, n=2; n<N; n<<=1, di >>= 1) {
- for (j=0; j<N; j+=(n<<1)) {
- indexW = 0;
- for (i=0; i<n; i++) {
- index1 = i + j;
- index2 = index1 + n;
- result.MulConj(Y[index2], W[indexW], 14);
- Y[index2].Sub(Y[index1], result);
- Y[index1].Add(Y[index1], result);
- indexW += di;
- }
- }
- }
- }
- };
