离散傅里叶变换及其快速算法

1、上机实验一：离散傅里叶变换及其快速算法一、设计目的 通过编写程序，深入理解快速傅里叶变换算法（FFT）的含义，完成FFT算法的软件实现。二、设计任务 利用时间抽取算法，编写基2的快速傅立叶变换（FFT）程序，并在FFT程序基础上编写快速傅立叶反变换（IFFT）。三、设计要求1、FFT和IFFT子程序相对独立、具有一般性，并加详细注释；2、验证例5-4，并能得到正确结果；四、设计条件C语言五、编程规则1）程序输入元素的数目为2的整数次幂，即N为2M幂，整个运算需要M级蝶形运算。 2）输入序列按二进制码位倒置排列，输出序列按自然顺序排列。 3）输出数据占用输入数据的存储单元。 4）每一级含N/2个

2、基本蝶形运算。第L级中有N/2L个群，群与群间隔为2L。同一级中各个群的系数W分布相同。处于第L级的群的系数是2L.8）对于第L级的蝶形运算，输入数据的间隔为2L-1。根据上述要求，设计程序源代码如下所示：#include#include#include#define swap(a,b) float T; T=(a);(a)=b;(b)=T;void fft(float A,float B,unsigned M) /蝶形运算程序，A存实部，B存虚部，M是级数unsigned long N,I,J,K,L,LE,LE1,P,Q,R;float Wr,Wi,Wlr,Wli,WTr,WTi,thet

3、a,Tr,Ti;N=1M;/N=1M表示N=2MJ=0;for (I=0;IN-1;I+)/for循环负责码位倒置存储if(J1;/ K=N1表示K=N/2 while (K=2&J=K)/while循环表示须向次高位进一位 J-=K;K=1;/K=1表示K=K/2J+=K;for(L=1;L=M;L+)/for循环为M级FFT运算,外层循环由级数L控制，执行M次LE=1L;/ LE=1L表示2L,是群间隔LE1=LE/2; /每个群的系数W数目 Wr=1.0; Wi=0.0; theta=(-1)*3.1415926536/LE1; Wlr=cos(theta); Wli=sin(theta

4、); for(R=0;RLE1;R+)/中层循环由群系数控制 for(P=R;PN-1;P+=LE)/R是群系数的编号，P、Q是基本蝶形运算两个输入数据在数组中的编号，循环每次完成同一个系数 的蝶形运算 Q=P+LE1;Tr=Wr*AQ-Wi*BQ;Ti=Wr*BQ+Wi*AQ;AQ=AP-Tr;BQ=BP-Ti;AP+=Tr;BP+=Ti;WTr=Wr;WTi=Wi;Wr=WTr*Wlr-WTi*Wli;Wi=WTr*Wli+WTi*Wlr;return;void main()/主程序int i,M,N,lb;cout请输入转换类别(FFT，请输入1； IFFT，请输入0)lb;cout请输

5、入序列长度NN;float *A= new floatN;float *B= new floatN;M=log(N)/log(2);cout请输入序列的实部endl;/输入序列实部for(i=0;iAi;cout请输入序列的虚部endl;/输入序列虚部for(i=0;iBi;cout setiosflags(ios:fixed);/输出格式控制cout您输入的序列为endl;cout setiosflags(ios:fixed);for(i=0;iN;i+)coutAi+jBiendl;coutendl;if(lb=0)for(i=0;iN;i+) Bi=Bi*(-1);fft(A,B,M);

6、for(i=0;iN;i+) Bi=Bi*(-1);for(i=0;iN;i+)Ai=Ai/N;Bi=Bi/N;if(lb=1)fft(A,B,M);cout转换后的序列为endl;/输出序列for(i=0;iN;i+)coutAi+j(Bi)endl;coutendl;delete A;delete B;2.实例验证序列计算FFT得：X(0): 6.00000+j0.00000X(1): 2.00000+j2.00000X(2): -6.00000+j0.00000X(3): 2.00000+j-2.00000这与例题计算结果相同，表示快速傅利叶正变换成功。2.2 将上述做IFFT得：X(0

7、): 1.00000+j0.00000X(1): 2.00000+j-0.00000X(2): -1.00000+j0.00000X(3): 4.00000+j0.00000这与已知相同，表示快速傅利叶反变换成功。实验二 ：IIR数字滤波器的设计及软件实现一、设计目的1 掌握设计IIR的数字滤波器的冲激响应不变法。2 掌握设计IIR数字滤波器的双线性变换法。3 掌握设计IIR数字滤波器的实现方法。二、设计任务按照课本实验五中的设计指标要求，编写IIR数字滤波器的程序。三、设计要求1、IIR数字滤波器程序相对独立，并加详细注释；2、观察不同频率信号序列在滤波前和滤波后的序列幅度变化，理解数字滤波

8、的含义；3、思考不同频率信号滤波差异的原因；四、设计条件 计算机、C语言环境五、实验要求与处理（一）设计IIR数字滤波器的冲击响应不变法的步骤：1.明确实际目标，由于需要设计符合指标的模拟滤波器Ha（S），因此，首先要将在数字域给出的设计指标转换为模拟域设计指标；2.设计符合指标要求得模拟滤波器Ha（S），即确定滤波器的阶数N和截止频率c 。3.将模拟滤波器Ha（S）变换为数字滤波器Ha（Z），首先计算Ha（S）在S平面左半平面N个极点Sk及其留数，然后代入下式得到Ha（Z）式中Ts为采样时间；4.验证数字滤波器的频率特性，首先计算数字滤波器的幅频特性然后选取包括通带截止角频率在内的若干频率点

9、，计算数字滤波器的幅频特性，并与设计指标进行比较。如果不满足要求，就要改进设计，重复上述步骤，直到满足设计要求。（二）设计IIR数字滤波器的双线性变换法的步骤为：1.明确设计指标，双向性变化法需要多设计指标进行预畸变，才能将数字域给出的实际指标变换为模拟域设计指标，预畸变公式为式中，Ts为采样时间，wp和ws分别为数字滤波器的通带截止角频率和阻带截止角频率， p和s分别为模拟滤波器的通带截止角频率和阻带截止角频率。2.设计符合指标要求的模拟滤波器Ha（S），即确定滤波器的阶数N和截止角频率c3.将模拟滤波器Ha（S）变换为数字滤波器H(Z)首先计算Ha（S）在S平面左半平面N个几点Sk并构成H

10、a（S）然后利用下是将Ha（S）转换为H（Z） S=4.验证数字滤波器的频率特性在设计出数字滤波器的系统函数H(Z)后，为实现数字滤波器算法，需要选择一种实现结构形式。常用的结构形式有直接I型，直接型，级联型和并联型等。数字网络中子系统的实现方法可以采用分差方程的递推解法。这里选择级联法。六、IIR数字滤波器的源程序1、IIR数字滤波器的冲击响应不变法的程序；#include#includevoid IIRDF(float A,unsigned long N);void fft(float A,float B,unsigned M);main()float A1024,B1024,C1024,

11、D1024,F1024; /Ai，Di存储输入序列实部，Bi，Fi存储输入序列虚部 unsigned long N,i ,M; printf(input M=); scanf(%d,&M); N=1M; printf(input A:); for(i=0;iN;i+) Ai=4*sin(100*3.1415926*1.25/1000)*i);Di=4*sin(100*3.1415926*1.25/1000)*i); for(i=0;iN;i+) Bi=0;Fi=0;printf(x(n):);for(i=0;iN;i+) printf(x(%d)=%f,i,Ai); /输出采样序列的值 fft

12、(A,B,M); /对输入的离散序列进行FFT，求其频谱 for(i=0;iN;i+) Ci=sqrt(Ai*Ai+Bi*Bi); /Ci为输入离散序列的频谱幅值 printf(x(k):);for(i=0;iN;i+) printf(X%d=%fn,i,Ci); /输出Cifor(i=0;iN;i+) Bi=0; IIRDF(D,N); /Ai已被占用，故用Di存储的输入序列通过低通滤波器，Fi同理 printf(y(n):); /通过低通滤波器后的输出序列for(i=0;iN;i+) printf(y(%d)=%fn,i,Di); fft(D,F,M); /对输出序列进行FFT，求其频谱p

13、rintf(y(k):);for(i=0;iN;i+) Ci=sqrt(Di*Di+Fi*Fi); /Ci为输出离散序列的频谱幅值for(i=0;iN;i+) printf(y(%d)=%fn,i,Ci); /输出Ci printf(n);void IIRDF(float A,unsigned long N)/*级联型低通滤波器*/unsigned long n ; float x3=0,0,0,y13=0,0,0,y23=0,0,0,y3=0,0,0; /y(0)表示y（n），y(1)表示y（n-1），y（2）表示y（n-2） for(n=0;nN;n+) x0=An; y10=1.3143

14、2*y11-0.71489*y12+0.08338*x0+0.16676*x1+0.08338*x2;/*第一级*/ x2=x1; x1=x0;y20=1.0541*y21-0.37534*y22+0.08338*y10+0.16676*y11+0.08338*y12;/*第二级*/ y12=y11; y11=y10; y0=0.94592*y1-0.23422*y2+0.08338*y20+0.16676*y21+0.08338*y22;/*第三级*/ y22=y21; y21=y20; y2=y1; y1=y0; An=y0;#define swap(a,b) T=(a); (a)=(b)

15、; (b)=Tvoid fft(float A,float B,unsigned M)/蝶形运算程序，A存实部，B存虚部，M是级数unsigned long N,I,J,K,L,LE,LE1,P,Q,R; float Wr,Wi,W1r,W1i,WTr,WTi,theta,Tr,Ti; N=1M; J=0;for (I=0;II)float T;swap(AI,AJ);swap(BI,BJ); K=N1;while (K=2&J=K)J-=K; K=1;J+=K;for(L=1;L=M;L+) /for循环为M级FFT运算LE=1L; / LE=1L表示2的L次幂，相当于本级运算的N LE1=

16、LE/2; /LE1=2的L-1次幂，相当于有L-1个Wn Wr=1.0; Wi=0.0; theta=(-1)*3.1415926536/LE1; /theta相当于-2/N， W1r=cos(theta); W1i=sin(theta);for(R=0;RLE1;R+) /for循环为每级有2的L-1次运算，R为系数序号，LE1为系数个数 for(P=R;PN;P+=LE) /for循环为蝶形的群运算 Q=P+LE1; Tr=Wr*AQ-Wi*BQ; Ti=Wr*BQ+Wi*AQ; AQ=AP-Tr; BQ=BP-Ti; AP+=Tr; BP+=Ti; WTr=Wr; WTi=Wi; Wr

17、=WTr*W1r-WTi*W1i; Wi=WTr*W1i+WTi*W1r; return;2、IIR数字滤波器的双线性变换法的程序：#includemath.h #define swap(a,b) temp=(a);(a)=(b);(b)=temp void fft(float A,float B,unsigned M);/该函数为快速傅里叶变换的函数void IIRDF(float A,unsigned long N) /IIR滤波器程序 int n; float x2=0,0,y13=0,0,0,y3=0,0,0; for(n=0;nN;n+) x0=An; /把A(n)赋给x（n）y10

18、=0.853973581*y11+0.073013209*x0+0.073013209*x1; /y1(0)表示，y1（1）表示，y（0）表示，y（1）表示，y（2）表示 x1=x0; /,因下个循环的x(n)参考点后移一位y0=1.831936826*y1-0.85480855*y2+0.005717931*y10+0.011435862*y11+0.005717931*y12;/ y12=y11; / y11=y10; / y2=y1; / y1=y0; / An=y0; /即位运算，即 main()int i,N,M;float A100,B100,C100;clrscr(); prin

19、tf(Input the number of samples x(n):); /输入采样序列的长度N scanf(%d,&N); for(i=0;iN;i+) printf(ninput x(%d): ,i); /输入采样序列x(n) 并赋给数组A scanf(%f,&Ai); Bi=Ai; /Bi、Ci分别表示x(i)的实部和虚部 Ci=0; /因是实序列，故Ci=0 M=ceil(log10(N)/log10(2); /M是下面做FFT运算的级数，,ceil()是向上取整函数 if(Npow(2,M) /如果序列的长度不是，在其后以0补足 for(i=N;ipow(2,M);i+) Ai=

20、0; Bi=0; Ci=0; IIRDF(A,N); /将x(n)用设计的低通滤波器滤波 fft(A,C,M,1); /将滤波后的序列做FFT fft(B,C,M,1); /将滤波前的序列做FFT printf(nFuDu: liu bo qian liu bo hou zeng yin); for(i=0;iN;i+) Bi=sqrt(Bi*Bi+Ci*Ci); /Bi是滤波前序列的幅值 Ai=sqrt(Ai*Ai+Ci*Ci); /Ai是滤波后序列的幅值 Ci=20*log10(Bi/Ai); /Ci是序列点的幅值的增益（dB） printf(x(%d): %.4f %.4f %.4fn,

21、i,Bi,Ai,Ci);七、实验分析（1）数据验证x1(t)=4sin(100t）x2(t)=4sin(200t）x3(t)=4sin(400t） Ts=1.25ms， T=80ms p=2/Ts*tan（wp/2）=519.87rad/s,s=815.24rad/s对于x1(t)=4sin(100t)，运行结果为input M=6input A:x(n):x(0)=0.000000 x(1)=1.530734 x(2)=2.828427 x(3)=3.695518 x(4)=4.000000 x(5)=3.695518 x(6)=2.828427 x(7)=1.530734 x(8)=0.0

22、00000 x(9)=-1.530733 x(10)=-2.828427 x(11)=-3.695518 x(12)=-4.000000 x(13)=-3.695518 x(14)=-2.828427 x(15)=-1.530734 x(16)=-0.000000 x(17)=1.530733 x(18)=2.828427 x(19)=3.695518 x(20)=4.000000 x(21)=3.695518 x(22)=2.828428 x(23)=1.530734 x(24)=0.000001 x(25)=-1.530733 x(26)=-2.828427 x(27)=-3.695518

23、 x(28)=-4.000000 x(29)=-3.695518 x(30)=-2.828428 x(31)=-1.530735 x(32)=-0.000001 x(33)=1.530733 x(34)=2.828427 x(35)=3.695518 x(36)=4.000000 x(37)=3.695518 x(38)=2.828428 x(39)=1.530735 x(40)=0.000001 x(41)=-1.530733 x(42)=-2.828426 x(43)=-3.695518 x(44)=-4.000000 x(45)=-3.695518 x(46)=-2.828428 x(4

24、7)=-1.530735 x(48)=-0.000001 x(49)=1.530733 x(50)=2.828426 x(51)=3.695518 x(52)=4.000000 x(53)=3.695519 x(54)=2.828428 x(55)=1.530735 x(56)=0.000002 x(57)=-1.530732 x(58)=-2.828426 x(59)=-3.695518 x(60)=-4.000000 x(61)=-3.695519 x(62)=-2.828428 x(63)=-1.530735y(n):y(0)=0.000000 y(1)=0.000887 y(2)=0.

25、009904 y(3)=0.053724 y(4)=0.190219 y(5)=0.498033y(6)=1.034063 y(7)=1.778526 y(8)=2.604401 y(9)=3.296343 y(10)=3.615747 y(11)=3.384228y(12)=2.550329 y(13)=1.214164 y(14)=-0.396109 y(15)=-1.983722 y(16)=-3.256230 y(17)=-3.989481y(18)=-4.067703 y(19)=-3.495075 y(20)=-2.382815 y(21)=-0.920360 y(22)=0.66

26、0368 y(23)=2.120864y(24)=3.248211 y(25)=3.881631 y(26)=3.932208 y(27)=3.394365 y(28)=2.347403 y(29)=0.945805y(30)=-0.601653 y(31)=-2.061926 y(32)=-3.212752 y(33)=-3.877175 y(34)=-3.951705 y(35)=-3.423175y(36)=-2.371328 y(37)=-0.956652 y(38)=0.604502 y(39)=2.073426 y(40)=3.225829 y(41)=3.886142y(42)=

27、3.954141 y(43)=3.419967 y(44)=2.365370 y(45)=0.951114 y(46)=-0.607521 y(47)=-2.073434y(48)=-3.223682 y(49)=-3.883314 y(50)=-3.951959 y(51)=-3.419121 y(52)=-2.365817 y(53)=-0.952307y(54)=0.606272 y(55)=2.072646 y(56)=3.223539 y(57)=3.883689 y(58)=3.952555 y(59)=3.419635y(60)=2.366067 y(61)=0.952269 y

28、(62)=-0.606502 y(63)=-2.072921Press any key to continue对于x2(t)=4sin(200t)，运行结果为input M=6input A:x(n):x(0)=0.000000 x(1)=2.828427 x(2)=4.000000 x(3)=2.828427 x(4)=0.000000 x(5)=-2.828427 x(6)=-4.000000 x(7)=-2.828427 x(8)=-0.000000 x(9)=2.828427 x(10)=4.000000 x(11)=2.828428 x(12)=0.000001 x(13)=-2.8

29、28427 x(14)=-4.000000 x(15)=-2.828428 x(16)=-0.000001 x(17)=2.828427 x(18)=4.000000 x(19)=2.828428 x(20)=0.000001 x(21)=-2.828426 x(22)=-4.000000 x(23)=-2.828428 x(24)=-0.000001 x(25)=2.828426 x(26)=4.000000 x(27)=2.828428 x(28)=0.000002 x(29)=-2.828426 x(30)=-4.000000 x(31)=-2.828428 x(32)=-0.00000

30、2 x(33)=2.828426 x(34)=4.000000 x(35)=2.828429 x(36)=0.000002 x(37)=-2.828426 x(38)=-4.000000 x(39)=-2.828429 x(40)=-0.000002 x(41)=2.828426 x(42)=4.000000 x(43)=2.828429 x(44)=0.000002 x(45)=-2.828425 x(46)=-4.000000 x(47)=-2.828429 x(48)=-0.000003 x(49)=2.828425 x(50)=4.000000 x(51)=2.828429 x(52)

31、=0.000003 x(53)=-2.828425 x(54)=-4.000000 x(55)=-2.828429 x(56)=-0.000003 x(57)=2.828425 x(58)=4.000000 x(59)=2.828429 x(60)=0.000003 x(61)=-2.828425 x(62)=-4.000000 x(63)=-2.828429y(n):y(0)=0.000000 y(1)=0.001640 y(2)=0.017590 y(3)=0.090329 y(4)=0.296509 y(5)=0.699064 y(6)=1.253904 y(7)=1.750245 y(

32、8)=1.872052 y(9)=1.382391 y(10)=0.321861 y(11)=-0.928060 y(12)=-1.795388 y(13)=-1.821651 y(14)=-0.945751 y(15)=0.414600 y(16)=1.556292 y(17)=1.866754 y(18)=1.171774 y(19)=-0.151266 y(20)=-1.371766 y(21)=-1.811925 y(22)=-1.231046 y(23)=0.034568 y(24)=1.260962 y(25)=1.749781 y(26)=1.228568 y(27)=0.006

33、567 y(28)=-1.205161 y(29)=-1.705868 y(30)=-1.210756 y(31)=-0.014551 y(32)=1.181934 y(33)=1.681051 y(34)=1.194744 y(35)=0.011249 y(36)=-1.174826 y(37)=-1.669348 y(38)=-1.184445 y(39)=-0.006079 y(40)=1.174258 y(41)=1.664906 y(42)=1.179012 y(43)=0.002114 y(44)=-1.175586 y(45)=-1.663816 y(46)=-1.176631

34、y(47)=0.000236 y(48)=1.176973 y(49)=1.663959 y(50)=1.175827 y(51)=-0.001395 y(52)=-1.177921 y(53)=-1.664377 y(54)=-1.175698 y(55)=0.001863 y(56)=1.178444 y(57)=1.664730 y(58)=1.175788 y(59)=-0.001997 y(60)=-1.178684 y(61)=-1.664950 y(62)=-1.175906 y(63)=0.002000Press any key to continue对于x3(t)=4sin(

35、400t)，运行结果为input M=6input A:x(n):x(0)=0.000000 x(1)=4.000000 x(2)=0.000000 x(3)=-4.000000 x(4)=-0.000000 x(5)=4.000000 x(6)=0.000001 x(7)=-4.000000 x(8)=-0.000001 x(9)=4.000000 x(10)=0.000001 x(11)=-4.000000 x(12)=-0.000001 x(13)=4.000000 x(14)=0.000002 x(15)=-4.000000 x(16)=-0.000002 x(17)=4.000000

36、 x(18)=0.000002 x(19)=-4.000000 x(20)=-0.000002 x(21)=4.000000 x(22)=0.000002 x(23)=-4.000000 x(24)=-0.000003 x(25)=4.000000 x(26)=0.000003 x(27)=-4.000000 x(28)=-0.000003 x(29)=4.000000 x(30)=0.000003 x(31)=-4.000000 x(32)=-0.000003 x(33)=4.000000 x(34)=0.000004 x(35)=-4.000000 x(36)=-0.000004 x(37

37、)=4.000000 x(38)=0.000004 x(39)=-4.000000 x(40)=-0.000004 x(41)=4.000000 x(42)=0.000005 x(43)=-4.000000 x(44)=-0.000005 x(45)=4.000000 x(46)=0.000005 x(47)=-4.000000 x(48)=-0.000005 x(49)=4.000000 x(50)=0.000005 x(51)=-4.000000 x(52)=-0.000006 x(53)=4.000000 x(54)=0.000006 x(55)=-4.000000 x(56)=-0.0

38、00006 x(57)=4.000000 x(58)=0.000006 x(59)=-4.000000 x(60)=-0.000006 x(61)=4.000000 x(62)=0.000007 x(63)=-4.000000y(n):y(0)=0.000000 y(1)=0.002319 y(2)=0.021597 y(3)=0.092564 y(4)=0.241948 y(5)=0.430788 y(6)=0.552539 y(7)=0.525251 y(8)=0.367739 y(9)=0.160862 y(10)=-0.029860 y(11)=-0.162063y(12)=-0.19

39、7903 y(13)=-0.135839 y(14)=-0.035351 y(15)=0.037470 y(16)=0.069588 y(17)=0.076267y(18)=0.054972 y(19)=0.006252 y(20)=-0.035272 y(21)=-0.039091 y(22)=-0.021809 y(23)=-0.012381y(24)=-0.005018 y(25)=0.013906 y(26)=0.026185 y(27)=0.012807 y(28)=-0.006215 y(29)=-0.005662y(30)=0.001329 y(31)=-0.005867 y(3

40、2)=-0.012988 y(33)=-0.001214 y(34)=0.012017 y(35)=0.005000 y(36)=-0.006347 y(37)=-0.000254 y(38)=0.008531 y(39)=-0.000268 y(40)=-0.010779 y(41)=-0.002312 y(42)=0.008993 y(43)=0.001811 y(44)=-0.008376 y(45)=-0.000641 y(46)=0.009472 y(47)=0.001246 y(48)=-0.009462 y(49)=-0.001664 y(50)=0.008904 y(51)=0.001230 y(52)=-0.009076 y(53)=-0.001146 y(54)=0.009309 y(55)=0.001392 y(56)=-0.009152 y(57)=-0.001362 y(58)=0.009080 y(59)=0.001245