| /* 函数名称:FFT() * 参数: * complex<double> * TD - 指向时域数组的指针 * complex<double> * FD - 指向频域数组的指针 * r -2的幂数,即迭代次数 * 返回值: 无。 * 说明:该函数用来实现快速傅立叶变换 */ void FFT(complex<double> * TD, complex<double> * FD, int r) { LONG count; // 傅立叶变换点数 int i,j,k; // 循环变量 int bfsize,p; // 中间变量 double angle; // 角度 complex<double> *W,*X1,*X2,*X; count = 1 << r; //傅立叶变换点数 // 分配运算所需存储器 W = new complex<double>[count / 2]; X1 = new complex<double>[count]; X2 = new complex<double>[count]; // 计算加权系数 for(i = 0; i < count / 2; i++) { angle = -i * PI * 2 / count; W[i] = complex<double> (cos(angle), sin(angle)); } // 将时域点写入X1 memcpy(X1, TD, sizeof(complex<double>) * count); // 采用蝶形算法进行快速傅立叶变换 for(k = 0; k < r; k++) { for(j = 0; j < 1 << k; j++) { bfsize = 1 << (r-k); for(i = 0; i < bfsize / 2; i++) { p = j * bfsize; X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2]; X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)]; } } X = X1; X1 = X2; X2 = X; } // 重新排序 for(j = 0; j < count; j++) { p = 0; for(i = 0; i < r; i++) { if (j&(1<<i)) { p+=1<<(r-i-1); } } FD[j]=X1[p]; } // 释放内存 delete W; delete X1; delete X2; } |
| #ifndef FFT_H #define FFT_H #include <complex> using namespace std; extern "C" void __declspec(dllexport) __stdcall FFT(complex<double> * TD, complex<double> * FD, int r); #define PI 3.1415926 #endif fft.cpp的源代码为: /* 文件名:fft.cpp */ #include "fft.h" void __stdcall FFT(complex<double> * TD, complex<double> * FD, int r) { …//读者提供的函数代码 } |