FFT离散傅氏变换的快速算法

int n, mmax, m, j, istep, i; double wtemp, wr, wpr, wpi, wi, theta; double tempr, tempi; n = nn << 1; j = 1; for (i = 1; i < n; i += 2) { if (j > i) { tempr = data[j]; data[j] = data[i]; data[i] = tempr; tempr = data[j+1]; data[j+1] = data[i+1]; data[i+1] = tempr; } m = n >> 1; while (m >= 2 && j > m) { j -= m; m >>= 1; } j += m; } mmax = 2; while (n > mmax) { istep = 2*mmax; theta = TWOPI/(isign*mmax); wtemp = sin*theta); wpr = *wtemp*wtemp; wpi = sin(theta); wr = ; wi = ; for (m = 1; m < mmax; m += 2) { for (i = m; i <= n; i += istep) { j =i + mmax; tempr = wr*data[j] - wi*data[j+1]; tempi = wr*data[j+1] + wi*data[j];
相关文档：
• • • • • • • • • •
更多相关文档请访问：
} // 将时域点写入 X1 memcpy(X1, TD, sizeof(complex) * count); // 采用蝶形算法进行快速 Fourier 变换 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] * W[i * (1<
FFT(离散傅氏变换的快速算法)
FFT(离散傅氏变换的快速算法) 目录 1 算法简介 2DFT 算法 3 源码表示 4MATLAB 中 FFT 的使用方法 1 算法简介编辑 FFT(Fast Fourier Transformation),即为快速傅氏变换,是离散傅氏变换的快 速算法,它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法 进行改进获得的。它对傅氏变换的理论并没有新的 FFT 算法图(Bufferfly 算法) 发现,但是对于在计算机系统或者说数字系统中应用离散傅立叶变换,可以说是 进了一大步。 设 x(n)为 N 项的复数序列,由 DFT 变换,任一 X(m)的计算都需要 N 次复数乘 法和 N-1 次复数加法,而一次复数乘法等于四次实数乘法和两次实数加法,一次复数 加法等于两次实数加法,即使把一次复数乘法和一次复数加法定义成一次“运 算”(四次实数乘法和四次实数加法),那么求出 N 项复数序列的 X (m),即 N 点 DFT 变换大约就需要 N^2 次运算。当 N=1024 点甚至更多的时候,需要 N2=1048576 次运算,在 FFT 中,利用 WN 的周期性和对称性,把一个 N 项序列(设 N=2k,k 为正整 数),分为两个 N/2 项的子序列,每个 N/2 点 DFT 变换需要(N/2)2 次运算,再用 N 次 运算把两个 N/2 点的 DFT 变换组合成一个 N 点的 DFT 变换。这样变换以后,总的 运算次数就变成 N+2*(N/2)^2=N+(N^2)/2。继续上面的例子, N=1024 时,总的 运算次数就变成了 525312 次,节省了大约 50%的运算量。而如果我们将这种“一 分为二” 的思想不断进行下去,直到分成两两一组的 DFT 运算单元,那么 N 点的 DFT 变 换就只需要 Nlog2N 次的运算,N 在 1024 点时,运算量仅有 10240 次,是先前的直 接算法的 1%,点数越多,运算量的节约就越大,这就是 FFT 的优越性。 2DFT 算法编辑 For length N input vector x, the DFT is a length N vector X, with elements
data[j] = data[i] - tempr; data[j+1] = data[i+1] - tempi; data[i] += tempr; data[i+1] += tempi; } wr = (wtemp = wr)*wpr - wi*wpi + wr; wi = wi*wpr + wtemp*wpi + wi; } mmax = istep; } } 在 C++环境下的源码 bool FFT(complex * TD, complex * FD, int r) { //一维快速 Fourier 变换。 //complex * TD ——指向时域数组的指针; complex * FD ——指向频域数 组的指针; r ——2 的幂数,即迭代次数 LONG count; // Fourier 变换点数 int i,j,k; // 循环变量 int bfsize,p; // 中间变量 double angle; // 角度 complex *W,*X1,*X2,*X; count = 1 << r; // 计算 Fourier 变换点数为 1 左移 r 位 W = new complex[count / 2]; X1 = new complex[count]; X2 = new complex[count]; // 分配运算所需存储器 // 计算加权系数(旋转因子 w 的 i 次幂表) for(i = 0; i < count / 2; i++) { angle = -i * PI * 2 / count; W[ i ] = complex (cos(angle), sin(angle));
N X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N. n=1 The inverse DFT (computed by IFFT) is given by N x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N. k=1 3 源码表示编辑 在 C 环境下的源码 源码(1): 注:亲测,这个版本无法运行,作者删改了重要内容[1]请参考源码(2) (see pages 507-508 of Numerical Recipes in C) Inputs: data[] : array of complex* data points of size 2*NFFT+1. data[0] is unused, * the n'th complex number x(n), for 0 <= n <= length(x)-1, is stored as: data[2*n+1] = real(x(n)) data[2*n+2] = imag(x(n)) if length(Nx) < NFFT, the remainder of the array must be padded with zeros nn : FFT order NFFT. Thisห้องสมุดไป่ตู้MUST be a power of 2 and >= length(x). isign: if set to 1, computes the forward FFT if set to -1, computes Inverse FFT - in this case the output values have to be manually normalized by multiplying with 1/NFFT. Outputs: data[] : The FFT or IFFT results are stored in data, overwriting the input. */ void four1(double data[], int nn, int isign) {