正在学数字信号处理,感觉上学期信号与系统学得不扎实,因为当时只是死记公式,这学期数信老师提倡动手实践,觉得自己在编程中对公式理解得更加深刻了。
以下是我写的fft,欢迎指教。
/*时间抽选基2fft及ifft算法c语言实现*/
/*author :junyi sun*/
/*copyright 2004-2005*/
/*mail:ccnusjy@yahoo.com.cn*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define n 1000
/*定义复数类型*/
typedef struct{
double real;
double img;
}complex;
complex x[n], *w; /*输入序列,变换核*/
int size_x=0; /*输入序列的大小,在本程序中仅限2的次幂*/
double pi; /*圆周率*/
int main(){
int i,method;
void fft(); /*快速傅里叶变换*/
void ifft();
void initw(); /*初始化变换核*/
void change(); /*变址*/
void add(complex a,complex b,complex *c); /*复数加法*/
void mul(complex a,complex b,complex *c); /*复数乘法*/
void sub(complex a,complex b,complex *c); /*复数减法*/
void divi(complex a,complex b,complex *c);/*复数除法*/
void output(); /*输出结果*/
system("cls");
pi=atan(1)*4;
printf("please input the size of x:\n");
scanf("%d",&size_x);
printf("please input the data in x[n]:\n");
for(i=0;i<size_x;i++)
scanf("%lf%lf",&x.real,&x.img);
initw();
printf("use fft(0) or ifft(1)?\n");
scanf("%d",&method);
if(method==0)
fft();
else
ifft();
output();
return 0;
}
/*快速傅里叶变换*/
void fft(){
int i=0,j=0,k=0,l=0;
complex up,down;
change();
for(i=0;i< (int)( log(size_x)/log(2) );i++){ /*一级蝶形运算*/
l=( 1<<i );
for(j=0;j<size_x;j+= (1<<l) ){ /*一组蝶形运算*/
for(k=0;k<l;k++){ /*一个蝶形运算*/
mul(x[j+k+l],w[size_x*k/2/l],&up);
add(x[j+k],up,&up);
mul(x[j+k+l],w[size_x*k/2/l],&down);
sub(x[j+k],down,&down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*快速傅里叶逆变换*/
void ifft(){
int i=0,j=0,k=0,l=size_x;
complex up,down;
for(i=0;i< (int)( log(size_x)/log(2) );i++){ /*一级蝶形运算*/
l/=2;
for(j=0;j<size_x;j+= (1<<l) ){ /*一组蝶形运算*/