资源描述
DSP课程作业用C语言编写FFT程序1,快速傅里叶变换FFT简介 快速傅氏变换(FFT),是离散傅氏变换的快速算法,它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现,但是对于在计算机系统或者说数字系统中应用离散傅立叶变换,可以说是进了一大步。 我们假设 x(n)为N项的复数序列,由DFT变换,任一X(m)的计算都需要N次复数乘法和N-1次复数加法,而一次复数乘法等于四次实数乘法和两次实数加法,一次复数加法等于两次实数加法,即使把一次复数乘法和一次复数加法定义成一次“运算”(四次实数乘法和四次实数加法),那么求出N项复数序列的X(m),即N点DFT变换大约就需要N2次运算。当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+N2/2。继续上面的例子,N=1024时,总的运算次数就变成了525312次,节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去,直到分成两两一组的DFT运算单元,那么N点的DFT变换就只需要Nlog2N次的运算,N在1024点时,运算量仅有10240次,是先前的直接算法的1%,点数越多,运算量的节约就越大,这就是FFT的优越性。2,FFT算法的基本原理 FFT算法的基本思想:利用DFT系数的特性,合并DFT运算中的某些项,吧长序列的DFT短序列的DFT,从而减少其运算量。 FFT算法分类:时间抽选法DIT: Decimation-In-Time;频率抽选法DIF: Decimation-In-Frequency按时间抽选的基-2FFT算法1、算法原理设序列点数 N = 2L,L 为整数。 若不满足,则补零。N为2的整数幂的FFT算法称基-2FFT算法。将序列x(n)按n的奇偶分成两组:则x(n)的DFT: 其中 再利用周期性求X(k)的后半部分: 复数乘法: 当N = 2L时,共有L级蝶形,每级N / 2个蝶形,每个蝶形有1次复数乘法2次复数加法。 2)、运算量 复数乘法:2、运算量 复数加法:比较DFT 3)、算法特点原位计算 倒位序 蝶形运算 对N = 2L点FFT,输入倒位序,输出自然序,第m级运算每个蝶形的两节点距离为 2m1第m级运算: 蝶形运算两节点的第一个节点为k值,表示成L位二进制数,左移L m位,把右边空出的位置补零,结果为r的二进制数。 蝶形运算两节点的第一个节点为k值,表示成L位二进制数,左移L m位,把右边空出的位置补零,结果为r的二进制数。 存储单元 输入序列x(n) : N个存储单元 系数 :N / 2个存储单元 3,快速傅立叶变换的C语言实现方法我们要衡量一个处理器有没有足够的能力来运行FFT算法,根据以上的简单介绍可以得出以下两点:1. 处理器要在一个指令周期能完成乘和累加的工作,因为复数运算要多次查表相乘才能实现。 2. 间接寻址,可以实现增/减1个变址量,方便各种查表方法。FFT要对原始序列进行反序排列,处理器要有反序间接寻址的能力。 #include #include #include #define N 1000 typedef struct double real; double img; complex; void fft(); void ifft(); void initW(); /*初始化变化核*/ void change(); void add(complex ,complex ,complex *); void mul(complex ,complex ,complex *); void sub(complex ,complex ,complex *); void divi(complex ,complex ,complex *); void output(); complex xN, *W;/*输出序列的值*/ int size_x=0;/*输入序列的长度,只限2的N次方*/ double PI; int main() int i,method; system(cls); PI=atan(1)*4;/*pi等于4乘以1.0的正切值*/ printf(Please input the size of x:n); /*输入序列的长度*/ scanf(%d,&size_x); printf(Please input the data in xN:(such as:5 6)n); /*输入序列对应的值*/ for(i=0;isize_x;i+) scanf(%lf %lf,&xi.real,&xi.img); initW(); /*选择FFT或逆FFT运算*/ printf(Use FFT(0) or IFFT(1)?n); scanf(%d,&method); if(method=0) fft(); else ifft(); output(); 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(size_x)/log(2) ;i+) /*一级蝶形运算*/ l=1i; for(j=0;jsize_x;j+= 2*l ) /*一组蝶形运算*/ for(k=0;kl;k+) /*一个蝶形运算*/ mul(xj+k+l,Wsize_x*k/2/l,&product); add(xj+k,product,&up); sub(xj+k,product,&down); xj+k=up; xj+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;jsize_x;j+= 2*l ) /*一组蝶形运算*/ for(k=0;kl;k+) /*一个蝶形运算*/ add(xj+k,xj+k+l,&up); up.real/=2;up.img/=2; sub(xj+k,xj+k+l,&down); down.real/=2;down.img/=2; divi(down,Wsize_x*k/2/l,&down); xj+k=up; xj+k+l=down; change(); /*初始化变化核*/ void initW() int i; W=(complex *)malloc(sizeof(complex) * size_x); for(i=0;isize_x;i+) Wi.real=cos(2*PI/size_x*i); Wi.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;i0 ) j=j1; if(ji) temp=xi; xi=xj; xj=temp; void output() /*输出结果*/ int i; printf(The result are as followsn); for(i=0;i=0.0001) printf(+%.4fjn,xi.img); else if(fabs(xi.img)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);
展开阅读全文