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

上机实验一：离散傅里叶变换及其快速算法一、设计目的通过编写程序，深入理解快速傅里叶变换算法（FFT）的含义，完成FFT算法的软件实现。二、设计任务利用时间抽取算法，编写基2的快速傅立叶变换（FFT）程序，并在FFT程序基础上编写快速傅立叶反变换（IFFT）。三、设计要求1、FFT和IFFT子程序相对独立、具有一般性，并加详细注释；2、验证例5-4，并能得到正确结果；四、设计条件C语言五、编程规则1）程序输入元素的数目为2的整数次幂，即N为2𝑀幂，整个运算需要M级蝶形运算。2）输入序列按二进制码位倒置排列，输出序列按自然顺序排列。3）输出数据占用输入数据的存储单元。4）每一级含N/2个基本蝶形运算。第L级中有N/2𝐿个群，群与群间隔为2𝐿。同一级中各个群的系数W分布相同。处于第L级的群的系数是2𝐿.8）对于第L级的蝶形运算，输入数据的间隔为2𝐿−1。输入A[I]，MI←0,J←0,N←2MI≥J？T←A[J]A[J]←A[I]A[I]←TK←N/2KJ?J←J←KK←K/2J←J+KI←I+1I=N-1?输出A[I]NoYesNoYesNo接蝶形运算程序Yes接码位倒置程序L=1LE=2LLE1=LE/2W=1W1=exp(-jπ/LE1)R=0P=RQ=P+LE1T=A[Q]*WA[Q]=A[P]-TA[P]=A[Q]+TP=P+LEPN-1?NoW=W*W1YesR=R+1RLE1-1?L=L+1NoYesLM?No结束Yes根据上述要求，设计程序源代码如下所示：#includemath.h#includeiostream.h#includeiomanip.h#defineswap(a,b){floatT;T=(a);(a)=b;(b)=T;}voidfft(floatA[],floatB[],unsignedM)//蝶形运算程序，A存实部，B存虚部，M是级数{unsignedlongN,I,J,K,L,LE,LE1,P,Q,R;floatWr,Wi,Wlr,Wli,WTr,WTi,theta,Tr,Ti;N=1M;//N=1M表示N=2^MJ=0;for(I=0;IN-1;I++)//for循环负责码位倒置存储{if(JI){swap(A[I],A[J]);swap(B[I],B[J]);}K=N1;//K=N1表示K=N/2while(K=2&&J=K)//while循环表示须向次高位进一位{J-=K;K=1;//K=1表示K=K/2}J+=K;}for(L=1;L=M;L++)//for循环为M级FFT运算,外层循环由级数L控制，执行M次{LE=1L;//LE=1L表示2^L,是群间隔LE1=LE/2;//每个群的系数W数目Wr=1.0;Wi=0.0;theta=(-1)*3.1415926536/LE1;Wlr=cos(theta);Wli=sin(theta);for(R=0;RLE1;R++)//中层循环由群系数控制{for(P=R;PN-1;P+=LE)//R是群系数的编号，P、Q是基本蝶形运算两个输入数据在数组中的编号，循环每次完成同一个系数的蝶形运算{Q=P+LE1;Tr=Wr*A[Q]-Wi*B[Q];Ti=Wr*B[Q]+Wi*A[Q];A[Q]=A[P]-Tr;B[Q]=B[P]-Ti;A[P]+=Tr;B[P]+=Ti;}WTr=Wr;WTi=Wi;Wr=WTr*Wlr-WTi*Wli;Wi=WTr*Wli+WTi*Wlr;}}return;}}voidmain()//主程序{inti,M,N,lb;cout请输入转换类别(FFT，请输入1；IFFT，请输入0)endl;//确定转换类别cinlb;cout请输入序列长度Nendl;cinN;float*A=newfloat[N];float*B=newfloat[N];M=log(N)/log(2);cout请输入序列的实部endl;//输入序列实部for(i=0;iN;i++){cinA[i];}cout请输入序列的虚部endl;//输入序列虚部for(i=0;iN;i++){cinB[i];}coutsetiosflags(ios::fixed);//输出格式控制cout您输入的序列为endl;coutsetiosflags(ios::fixed);for(i=0;iN;i++)coutA[i]+jB[i]endl;coutendl;if(lb==0){for(i=0;iN;i++){B[i]=B[i]*(-1);}fft(A,B,M);for(i=0;iN;i++){B[i]=B[i]*(-1);}for(i=0;iN;i++){A[i]=A[i]/N;B[i]=B[i]/N;}}if(lb==1)fft(A,B,M);cout转换后的序列为endl;//输出序列for(i=0;iN;i++)coutA[i]+j(B[i])endl;coutendl;delete[]A;delete[]B;}2.实例验证序列(){1,2,1,4}xn计算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将上述()Xk做IFFT得：X(0):1.00000+j0.00000X(1):2.00000+j-0.00000X(2):-1.00000+j0.00000X(3):4.00000+j0.00000这与已知()xn相同，表示快速傅利叶反变换成功。实验二：IIR数字滤波器的设计及软件实现一、设计目的1掌握设计IIR的数字滤波器的冲激响应不变法。2掌握设计IIR数字滤波器的双线性变换法。3掌握设计IIR数字滤波器的实现方法。二、设计任务按照课本实验五中的设计指标要求，编写IIR数字滤波器的程序。三、设计要求1、IIR数字滤波器程序相对独立，并加详细注释；2、观察不同频率信号序列在滤波前和滤波后的序列幅度变化，理解数字滤波的含义；3、思考不同频率信号滤波差异的原因；四、设计条件计算机、C语言环境五、实验要求与处理（一）设计IIR数字滤波器的冲击响应不变法的步骤：1.明确实际目标，由于需要设计符合指标的模拟滤波器Ha（S），因此，首先要将在数字域给出的设计指标转换为模拟域设计指标；2.设计符合指标要求得模拟滤波器Ha（S），即确定滤波器的阶数N和截止频率Ωc。3.将模拟滤波器Ha（S）变换为数字滤波器Ha（Z），首先计算Ha（S）在S平面左半平面N个极点Sk及其留数，然后代入下式得到Ha（Z）NkzTsZsk11e-1Ak?)(H式中Ts为采样时间；4.验证数字滤波器的频率特性，首先计算数字滤波器的幅频特性)()(zHeeHjzj然后选取包括通带截止角频率在内的若干频率点，计算数字滤波器的幅频特性，并与设计指标进行比较。如果不满足要求，就要改进设计，重复上述步骤，直到满足设计要求。（二）设计IIR数字滤波器的双线性变换法的步骤为：1.明确设计指标，双向性变化法需要多设计指标进行预畸变，才能将数字域给出的实际指标变换为模拟域设计指标，预畸变公式为)2tan(2),2tan(2ssdSwTwT式中，Ts为采样时间，wp和ws分别为数字滤波器的通带截止角频率和阻带截止角频率，Ωp和Ωs分别为模拟滤波器的通带截止角频率和阻带截止角频率。2.设计符合指标要求的模拟滤波器Ha（S），即确定滤波器的阶数N和截止角频率Ωc3.将模拟滤波器Ha（S）变换为数字滤波器H(Z)首先计算Ha（S）在S平面左半平面N个几点Sk并构成Ha（S）然后利用下是将Ha（S）转换为H（Z）S=1112ZHZTs4.验证数字滤波器的频率特性在设计出数字滤波器的系统函数H(Z)后，为实现数字滤波器算法，需要选择一种实现结构形式。常用的结构形式有直接I型，直接Ⅱ型，级联型和并联型等。数字网络中子系统的实现方法可以采用分差方程的递推解法。这里选择级联法。六、IIR数字滤波器的源程序1、IIR数字滤波器的冲击响应不变法的程序；#includemath.h#includestdio.hvoidIIRDF(floatA[],unsignedlongN);voidfft(floatA[],floatB[],unsignedM);main(){floatA[1024],B[1024],C[1024],D[1024],F[1024];//A[i]，D[i]存储输入序列实部，B[i]，F[i]存储输入序列虚部unsignedlongN,i,M;printf(inputM=);scanf(%d,&M);N=1M;printf(inputA[]:);for(i=0;iN;i++){A[i]=4*sin((100*3.1415926*1.25/1000)*i);D[i]=4*sin((100*3.1415926*1.25/1000)*i);}for(i=0;iN;i++){B[i]=0;F[i]=0;}printf(x(n):);for(i=0;iN;i++)printf(x(%d)=%f,i,A[i]);//输出采样序列的值fft(A,B,M);//对输入的离散序列进行FFT，求其频谱for(i=0;iN;i++)C[i]=sqrt(A[i]*A[i]+B[i]*B[i]);//C[i]为输入离散序列的频谱幅值printf(x(k):);for(i=0;iN;i++)printf(X[%d]=%f\n,i,C[i]);//输出C[i]for(i=0;iN;i++)B[i]=0;IIRDF(D,N);//A[i]已被占用，故用D[i]存储的输入序列通过低通滤波器，F[i]同理printf(y(n):);//通过低通滤波器后的输出序列for(i=0;iN;i++)printf(y(%d)=%f\n,i,D[i]);fft(D,F,M);//对输出序列进行FFT，求其频谱printf(y(k):);for(i=0;iN;i++)C[i]=sqrt(D[i]*D[i]+F[i]*F[i]);//C[i]为输出离散序列的频谱幅值for(i=0;iN;i++)printf(y(%d)=%f\n,i,C[i]);//输出C[i]printf(\n);}voidIIRDF(floatA[],unsignedlongN)/*级联型低通滤波器*/{unsignedlongn;floatx[3]={0,0,0},y1[3]={0,0,0},y2[3]={0,0,0},y[3]={0,0,0};//y(0)表示y（n），y(1)表示y（n-1），y（2）表示y（n-2）for(n=0;nN;n++){x[0]=A[n];y1[0]=1.31432*y1[1]-0.71489*y1[2]+0.08338*x[0]+0.16676*x[1]+0.08338*x[2];/*第一级*/x[2]=x[1];x[1]=x[0];y2[0]=1.0541*y2[1]-0.37534*y2[2]+0.08338*y1[0]+0.16676*y1[1]+0.08338*y1[2];/*第二级*/y1[2]=y1[1];y1[1]=y1[0];y[0]=0.94592*y[1]-0.23422*y[2]+0.08338*y2[0]+0.16676*y2[1]+0.08338*y2[2];/*第三级*/y2[2]=y2[1];y2[1]=y2[0];y[2]=y[1];y[1]=y[0];A[n]=y[0];}}#defineswap(a,b)T=(a);(a)=(b);(b)=Tvoidfft(