fft.rar_FFT算法_MFCC_fft_fft.h

/************************************************************************** * 文件名：FFT&IFFT.cpp * * 函数： * * FFT() - 快速傅立叶变换 * IFFT() - 快速傅立叶反变换 * FFT2() - 图像的快速傅立叶变换 * IFFT2() - 图像的快速傅立叶反变换 * *************************************************************************/ #include "stdafx.h" #include "FFT&IFFT.h" #include <math.h> #include <complex> using namespace std; // 常数π #define PI 3.1415926535 /************************************************************************* * * 函数名称： * FFT() * * 参数: * complex<double> * TD - 指向时域数组的指针 * complex<double> * FD - 指向频域数组的指针 * r - 2的幂数，即迭代次数，2的r次方=N * * 说明: * 该函数用来实现快速付立叶变换。 * ************************************************************************/ 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; //1左移r位，即2的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++) //FFT的级数 { for(j = 0; j < 1 << k; j++) //组数 { bfsize = 1 << (r-k); //本级的跨度*2 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; } /************************************************************************* * * 函数名称： * IFFT() * * 参数: * complex<double> * FD - 指向频域值的指针 * complex<double> * TD - 指向时域值的指针 * r - 2的幂数，即迭代次数，2的r次方=N * * 说明: * 该函数用来实现快速付立叶反变换。 * ************************************************************************/ void IFFT(complex<double> * FD, complex<double> * TD, int r) { // 付立叶变换点数 long count; // 循环变量 int i; complex<double> *X; // 计算付立叶变换点数 count = 1 << r; //1左移r位，即2的r次方 // 分配运算所需存储器 X = new complex<double>[count]; // 将频域点写入X memcpy(X, FD, sizeof(complex<double>) * count); // 求共轭 for(i = 0; i < count; i++) { X[i] = complex<double> (X[i].real(), -X[i].imag()); } // 调用快速付立叶变换 FFT(X, TD, r); // 求时域点的共轭 for(i = 0; i < count; i++) { TD[i] = complex<double> (TD[i].real() / count, -TD[i].imag() / count); } // 释放内存 delete X; } /************************************************************************* * * 函数名称： * FFT2() * * 参数: * complex<double> * TD - 指向时域数组的指针 * complex<double> * FD - 指向频域数组的指针 * width, height - 2维变换的宽度和高度 * * 说明: * 该函数用来实现对图像的快速付立叶变换。 * ************************************************************************/ void FFT2(complex<double> * TD, complex<double> * FD, long width, long height) { // 求宽度、高度的2的幂(这里暂时要求输入图像的宽度和高度都是2的整数次方) int wr=0; // 2的幂数，即FFT迭代次数，2的r次方=width wr = log(width)/log(2); int hr=0; // 2的幂数，即FFT迭代次数，2的r次方=height hr = log(height)/log(2); // 循环变量 long i, j; for(i = 0; i < height; i++) { // 对x方向进行快速付立叶变换 FFT(&TD[width * i], &FD[width * i], wr); //每一行的FFT } // 保存变换结果,行列互换便于y方向变换 for(i = 0; i < height; i++) { for(j = 0; j < width; j++) { TD[i + height * j] = FD[j + width * i]; } } //临时变量 complex<double> *FDtemp; FDtemp = new complex<double>[width * height]; //保存频域变换结果 for(i = 0; i < width; i++) { // 对y方向进行快速付立叶变换 FFT(&TD[i * height], &FDtemp[i * height], hr); //每一列的FFT } //将FDtemp再行列互换一次才是最终结果 for(i = 0; i < height; i++) { for(j = 0; j < width; j++) { FD[i * width + j] = FDtemp[j * height + i]; } } //释放内存 delete FDtemp; } /************************************************************************* * * 函数名称： * IFFT2() * * 参数: * complex<double> * FD - 指向频域值的指针 * complex<double> * TD - 指向时域值的指针 * width, height - 2维变换的宽度和高度 * * 说明: * 该函数用来实现对图像的快速付立叶反变换。 * ************************************************************************/ void IFFT2(complex<double> * FD, complex<double> * TD, long width, long height) { // 求宽度、高度的2的幂(这里暂时要求输入图像的宽度和高度都是2的整数次方) int wr=0; // 2的幂数，即FFT迭代次数，2的r次方=width wr = log(width)/log(2); int hr=0; // 2的幂数，即FFT迭代次数，2的r次方=height hr = log(height)/log(2); // 循环变量 long i, j; for(i = 0; i < height; i++) { // 对x方向进行快速付立叶反变换 IFFT(&FD[i * width], &TD[i * width], wr); //每一行的IFFT } // 保存变换结果,行列互换便于y方向变换 for(i = 0; i < height; i++) { for(j = 0; j < width; j++) { FD[i + height * j] = TD[j + width * i]; } } //临时变量 complex<double> *TDtemp; TDtemp = new complex<double>[width * height]; //保存时域变换结果 for(i = 0; i < width; i++) { // 对y方向进行快速付立叶反变换 IFFT(&FD[height * i], &TDtemp[height * i], hr); //每一列的IFFT } //将TDtemp再行列互换一次才是最终结果 for(i = 0; i < height; i++) { for(j = 0; j < width; j++) { TD[i * width + j] = TDtemp[j * height + i]; } } //释放内存 delete TDtemp; }