C语言编写FFT程序

1、DSP课程作业 用C语言编写FFT程序 1快速傅里叶变换FFT简介 快速傅氏变换（FFT）,是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、 偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并 没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说 是进了 一大步。 我们假设x（n）为N项的复数序列，由DFT变换，任一 X（ m）的计算都需要 N次复 数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复 数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四 次实数乘法和四次实数加法），那么求

2、出N项复数序列的X（ m ,即N点DFT变换大约就 需要NT次运算。当 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+（N/2）2=N+N2/2。 继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%勺运算 量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组

3、的DFT运 算单元，那么N点的DFT变换就只需要 Nlog2N次的运算，N在1024点时，运算量仅有 10240次，是先前的直接算法的1%点数越多，运算量的节约就越大，这就是 FFT的优 越性。 2, FFT算法的基本原理 FFT 算法的基本思想：利用DFT系数的特性，合并DFT运算中的某些项， 吧长序列的 DFT-短序列的DFT,从而减少其运算量。 FFT 算法分类：时 间 抽选法 DIT: Decimation-ln-Time;频率抽选法 DIF: Decimati on-I n-F reque ncy 按时间抽选的基-2FFT算法 1、算法原理 设序列点数N = 2L , L为整数。 若

4、不满足，则补零。 N为2的整数幕的 FFT算法称基-2FFT算法。将序列 x（n）按n 的奇偶分成两组 x 2r r r N d 0,1,., 1 x 2r 1 x2 r 2 则 x(n)的 DFT: N 1 N nk 1 nk N 1 nk X k x n WN x n WN x n WN n 0n 0n 0 2 i x1(r)W；k r 02 N I 2 I r 0 x(2r)叫 x(2r)W 2 (k 0,1,.； 1) Xi (A) -V2(A) NN彳 2I2 1 x 2r wN2rk x 2r 2r 1 k 1 Wn 丁1 r 0 N 1 22 rk k 2 rk X1 rWn W

5、N X2 r Wn r 0 r 0 N 1 N 2 1 X1 r WN2 WNk X2 r W,/2 r 0 r 0 k X1 k WnX2 k (r,k 0,1,.N 1) 2 其中 X!(k) Xi(k) 再利用周期性求X(k)的后半部分： QXi k汛2 k是以N为周期的 X1 k NX1 k X2 k N X2 k 2 2 k NN 又Wn 刁 Wn2W,w, k X(k) Xi(k) WnX2(G Nk X(k ) Xi(k) WzX2(k) 图4时闻抽选法蝶形运算流图符号 n为偶数 n为奇数 V(0) 点 V(2) DI T XiO) 4 v1) x；(! )=a (3) 点 DF

6、T a-(3)=a7) ffl4 2按时间抽选，将仆V点DFT分鹏为 两个22总DFT Vz( 2) 分解后的运算量: 妲数乘法 复数加法 一个N/2点 DFT (V/2)2 N/2(N/2-l) 两个N/2点DF1 N2/2 /V(V/2 1) 个蝶形 1 2 和2个蝶形 /V/2 N 总计 /VJ/2 + A/2 *Af2/2 N(N/2-)-N f/2 运STSWT近半! 2 ）、运算量 当N = 2L时，共有 L级蝶形, 每级 N / 2个蝶形， 每 N N 复数乘法： 叶 L g N 2 2 复数加法： aF NL Nlog2N 比较DFT mF (DFT ) N 2 2 N m f

7、 ( FFT ) N . 2 log 2 N log 2 N 3 )、算法特点 原位计算 蝶形运算两节点的第一个节点为 把右边空出的位置补零，结果为 k值，表示成L位二进制数，左移 L r的二进制数。 m位, Xm(k) Xmi(k) Xmi(j)WN Xm(j) Xmi(k) Xmi(j)WN A i (/) W Xi 一计一,匕(/)=匕恋)_K泊呎 你-I 图47按时间抽选蝶形运算结构 倒位序 蝶形运算 对N = 2L点FFT，输入倒位序，输出自然 序， 第m级运算每个蝶形的两节点距离为2m - 倒位序 自然序 000 0 0 000 100 4 1 001 010 2 2 010 11

8、0 6 3 ()11 001 1 4 100 101 5 5 101 011 3 6 110 111 7 7 111 m 1 r Xm(k) Xmi(k) Xmi(k 2 )Wn Xm(k 2m1) Xmi(k) Xmi(k 2m1)wN .VnUAU因一兀*朋叭 V(A)-Vm.|(A)+.；ll0 ,Yju_?(/) * ffl 4-7按吋间抽选蝶形运算结构 wN的确定 蝶形运算两节点的第一个节点为k值，表示成L位二进制数，左移L -m 位，把右边空出的位置补零，结果为 r的二进制数。 存储单元 输入序列x(n) : N个存储单元 系数WN : N / 2个存储单元 3，快速傅立叶变换的C

9、语言实现方法 我们要衡量一个处理器有没有足够的能力来运行FFT算法，根据以上的简单介绍可以得 出以下两点: 1. 处理器要在一个指令周期能完成乘和累加的工作，因为复数运算要多次查表相乘 才能实现。 2. 间接寻址，可以实现增/减1个变址量，方便各种查表方法。FFT要对原始序列进 行反序排列，处理器要有反序间接寻址的能力。 #i nclude #in clude #in clude #defi ne N 1000 typedef struct double real; double img; complex; void fft(); void ifft(); void initW(); /* 初

10、始化变化核 */ 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=O;/* 输入序列的长度，只限 2的N次方*/ double PI; int main() int i,method; syste

11、m(cls); PI=atan(1)*4;/*pi等于4乘以 1.0 的正切值 */ printf(Please input the size of x:n); /* 输入序列的长度 */ scanf(%d, printf(Please input the data in xN:(such as:5 6)n); /* 输入序列对应的值 */ for(i=0;isize_x;i+) scanf(%lf %lf, initW(); /*选择FFT或逆FFT运算*/ printf(Use FFT(0) or IFFT(1)?n); scanf(%d, if(method=0) fft(); else

12、 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, add(xj+k,product, sub(xj+k,product, xj

13、+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.real/=2;up.img/=2; sub(xj+k,xj+k+l, down.real/=2;down.img/=2; divi(down,Wsize

14、_x*k/2/l, xj+k=up; xj+k+l=down; change(); /* 初始化变化核 */ void initW() int i; size_x); W=(complex *)malloc(sizeof(complex) 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 )

15、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