已经讲完了, fft 没有太多的东西好讲。这样吧, 我放出定点c++ 的 fft 以及逆变换的程序,供大家参考
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;
}
}
}
}
};
|