快速傅里叶变换FFT
DFT是信号分析与处理中的一种重要变换。但直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。 1.直接计算DFT 长度为N的有限长序列x(n)的DFT为：
2.减少运算量的思路和方法 思路：N点DFT的复乘次数等于N2。把N点DFT分解为几个较短的DFT，可使乘法次数大大减少。另外，旋转因子WmN具有周期性和对称性。 (考虑x(n)为复数序列的一般情况，对某一个k值，直接按上式计算X(k)值需要N次复数乘法、(N-1)次复数加法.) 方法： 分解N为较小值：把序列分解为几个较短的序列，分别计算其DFT值，可使乘法次数大大减少； 利用旋转因子WNk的周期性、对称性进行合并、归类处理，以减少DFT的运算次数。 周期性：
对称性：
3.FFT算法思想 不断地把长序列的DFT分解成几个短序列的DFT，，并利用旋转因子的周期性和对称性来减少DFT的运算次数。
再次分解，对N=8点，可分解三次。
c语言程序:
{double real;double img; }complex;void fft();/*快速傅里叶变换*/void ifft();void initW();void change();sub(output(); size_x = 0;/*输入序列的长度，只限2的N次方*/double PI;int main(){int i, method;system(“cls”);PI = atan(1.0) * 4;printf(“Please input the size of x: \n”);/*输入序列的长度*/scanf(“%d”, &size_x);printf(“Please input the data in x[N](such as:5 6):\n”);/*输入序列对应的值*/for (i = 0; i<size_x; i++)scanf(“%lf %lf”, &x[i].real, &x[i].img);initW();/*选择FFT或逆FFT运算*/printf(“Use FFT(0) or IFFT(1) ?\n”);scanf(“%d”, &method);if (method == 0)fft();elseifft();output();system(“pause”);return 0;}/*进行基-2 FFT运算*/void fft(){int i = 0, j = 0, k = 0, l = 0;complex up, down, product;change();for (i = 0; i < log((float)size_x) / log(2.0); i++){l = 1 << i;for (j = 0; j<size_x; j += 2 * l){for (k = 0; k<l; k++){mul(x[j + k + l], W[size_x*k / 2 / l], &product);add(x[j + k], product, &up);sub(x[j + k], product, &down);x[j + k] = up;x[j + k + l] = down;}}} }void ifft(){int i = 0, j = 0, k = 0;complex up, down;for (i = 0; i < log((float)size_x) / log(2.0); i++) /*蝶形运算*/{int l = size_x;l/= 2;for (j = 0; j < size_x; j += 2 * l){for (k = 0; k < l; k++){add(x[j + k], x[j + k + l], &up);up.real /= 2;up.img /= 2;sub(x[j + k], x[j + k + l], &down);down.real /= 2;down.img /= 2;divi(down, W[size_x*k / 2 / l], &down);x[j + k] = up;x[j + k + l] = down;}}}change();}/*初始化变化核*/void initW(){int i;W = (complex *)malloc(sizeof(complex)* size_x);for (i = 0; i < size_x; i++){W[i].real = cos(2 * PI / size_x*i);W[i].img = -1 * sin(2 * PI / size_x*i);}}/*变址计算，将x(n)码位倒置*/void change(){complex temp;unsigned short i = 0, j = 0, k = 0;double t;for (i = 0; i<size_x; i++){k = i;j = 0;t = (log((float)size_x) / log(2.0));while ((t–)>0){j = j << 1;j |= (k & 1);k = k >> 1;}if (j > i){temp = x[i];x[i] = x[j];x[j] = temp;}}}void output()/*输出结果*/{int i;printf(“The result are as follows: \n”);for (i = 0; i<size_x; i++){printf(“%.4f”, x[i].real);if (x[i].img >= 0.0001)printf(“+%.4fi\n”, x[i].img);else if(fabs(x[i].img)<0.0001)printf(“\n”);elseprintf(“%.4fi\n”, x[i].img);}}void add(complex a, complex b, complex *c){c->real = a.real + b.real;c->img = a.img + b.img;}void mul(complex a, complex b, complex *c){c->real = a.real*b.real- a.img*b.img;c->img = a.real*b.img+ a.img*b.real;}void sub(complex a, complex b, complex *c){c->real = a.real – b.real;c->img = a.img – b.img;}void divi(complex a, complex b, complex *c){c->real = (a.real*b.real + a.img*b.img) / (b.real*b.real + b.img*b.img);c->img = (a.img*b.real – a.real*b.img) / (b.real*b.real + b.img*b.img);}可是我要如何在浅薄的纸上为你画上我所有的命轮？
