您好,欢迎访问三七文档
当前位置:首页 > 商业/管理/HR > 其它文档 > 信号分析与处理实验报告
华北电力大学实验报告||实验名称FFT的软件实现实验(Matlab)IIR数字滤波器的设计课程名称信号分析与处理||专业班级:电气化1308学生姓名:袁拉麻加学号:201301000827成绩:指导教师:杨光实验日期:2015-12-17快速傅里叶变换实验一、实验目的及要求通过编写程序,深入理解快速傅里叶变换算法(FFT)的含义,完成FFT和IFFT算法的软件实现。二、实验内容利用时间抽取算法,编写基2点的快速傅立叶变换(FFT)程序;并在FFT程序基础上编写快速傅里叶反变换(IFFT)的程序。三:实验要求1、FFT和IFFT子程序相对独立、具有一般性,并加详细注释;2、验证例6-4,并能得到正确结果。3、理解应用离散傅里叶变换(DFT)分析连续时间信号频谱的数学物理基础。四、实验原理:a.算法原理1、程序输入序列的元素数目必须为2的整数次幂,即N=2M,整个运算需要M级蝶形运算;2、输入序列应该按二进制的码位倒置排列,输出序列按自然序列排列;3、每个蝶形运算的输出数据军官占用其他输入数据的存储单元,实现“即位运算”;4、每一级包括N/2个基本蝶形运算,共有M*N/2个基本蝶形运算;5、第L级中有N/2L个群,群与群的间隔为2L。6、处于同一级的各个群的系数W分布相同,第L级的群中有2L-1个系数;7、处于第L级的群的系数是(p=1,2,3,…….,2L-1)而对于第L级的蝶形运算,两个输入数据的间隔为2L-1。b.码位倒置程序流程图开始检测A序列长度nk=0j=1x1(j)=bitget(k,j);j=j+1jm?YNx1=num2str(x1);y(k+1)=bin2dec(x1);clearx1k=k+1c.蝶形运算程序流程图kn?YN按新序号生产输入序列B(过程略)五、程序代码与实验结果接码位倒置程序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结束Yesa.FFT程序:%%clearall;closeall;clc;%输入数据%A=input('输入x(n)序列','s');A=str2num(A);%A=[1,2,-1,4];%测试数据%%%%校验序列,%n=length(A);m=log2(n);if(fix(m)~=m)disp('输入序列长度错误,请重新输入!');A=input('输入x(n)序列','s');A=str2num(A);elsedisp('输入正确,请运行下一步')end%%%码位倒置%fork=0:n-1forj=1:m%取M位的二进制数%x1(j)=bitget(k,j);%倒取出二进制数%endx1=num2str(x1);%将数字序列转化为字符串%y(k+1)=bin2dec(x1);%二进制序列转化为十进制数%clearx1endfork=1:nB(k)=A(y(k)+1);%时间抽取序列%endclearA%%%计算%forL=1:m%分解为M级进行运算%LE=2^L;%第L级群间隔为2^L%LE1=2^(L-1);%第L级中共有2^(L-1)个Wn乘数,进行运算蝶运算的两数序号相隔LE1%W=1;W1=exp(-1i*pi/LE1);forR=1:LE1%针对第R个Wn系数进行一轮蝶运算,共进行LE1次%forP=R:LE:n%每个蝶的大小为LE%Q=P+LE1;T=B(Q)*W;B(Q)=B(P)-T;B(P)=B(P)+T;endW=W*W1;endendB%输出X(k)%%%验证结果:例6-4b.IFFT程序:%%clearall;closeall;clc;%输入数据%A=input('输入X(k)序列','s');A=str2num(A);%A=[6,2+2i,-6,2-2i];%测试数据%%%%校验序列,%n=length(A);m=log2(n);if(fix(m)~=m)disp('输入序列长度错误,请重新输入!');A=input('输入x(n)序列','s');A=str2num(A);elsedisp('输入正确,请运行下一步')end%%%码位倒置%fork=0:n-1forj=1:m%取M位的二进制数%x1(j)=bitget(k,j);%倒取出二进制数%endx1=num2str(x1);%将数字序列转化为字符串%y(k+1)=bin2dec(x1);%二进制序列转化为十进制数%clearx1endfork=1:nB(k)=A(y(k)+1);%时间抽取序列%endclearA%%%计算%forL=1:m%分解为M级进行运算%LE=2^L;%第L级群间隔为2^L%LE1=2^(L-1);%第L级中共有2^(L-1)个Wn乘数,进行运算蝶运算的两数序号相隔LE1%W=1;W1=exp(-1i*pi/LE1);forR=1:LE1%针对第R个Wn系数进行一轮蝶运算,共进行LE1次%forP=R:LE:n%每个蝶的大小为LE%Q=P+LE1;T=B(Q)*W;B(Q)=B(P)-T;B(P)=B(P)+T;endW=W*W1;endendB=conj(B);%取共轭%B=B/n%输出x(n)%验证结果:六、实验心得与结论本次实验借助于Matlab软件,我避开了用C平台进行复杂的复数运算,在一定程度上简化了程序,并添加了简单的检错代码,码位倒置我通过查阅资料,使用了一些函数,涉及到十-二进制转换,数字-文本转换,二-文本转换,相对较复杂,蝶运算我参考了书上了流程图,做些许改动就能直接实现。通过本次实验,我对FFT的过程更加熟悉了,尤其是WN系数的特性,可以通过累乘来实现W系数的更新。IIR数字滤波器的设计实验一、实验目的1.掌握实现模拟信号数字化的编程方法。2.掌握实现时间抽取快速傅里叶变换(FFT)的编程方法。3.掌握设计IIR数字滤波器的双线性变换法。4.加深对信号分析与处理的理解。二、实验内容1.生成三个不同频率的正弦信号1()xt、2()xt、3()xt,将它们叠加成一个信号()xt。用满足采样定理的采样频率对模拟信号()xt进行采样,得到离散时间序列()xn。2.编制子程序实现FFT运算,并对()xn进行FFT运算,观察它的频谱。3.设计满足已知的设计指标的IIR数字滤波器,并用它对()xn进行滤波处理。对滤波后得到的波形再做FFT运算,观察滤波后的信号频谱。三、实验结果1.生成的信号图1原始信号及合成信号的时域波形图2原始信号的幅度谱2.滤波器的性能检验图3模拟滤波器的频率特性图4数字滤波器的频率特性3.滤波后的信号图5滤波前后信号的时域波形对比图6滤波前后信号的幅度谱对比附录附录一:FFT子程序function[A]=fftlzr(A)N=length(A);M=log2(N);j=complex(0,1);iffix(M)~=MA=[A(:);zeros(2^(fix(M)+1)-N,1)];endN=length(A);M=log2(N);I=0;J=0;L=1;whileI~=N-1ifIJt=A(J+1);A(J+1)=A(I+1);A(I+1)=t;endK=N/2;whileK=JJ=J-K;K=K/2;endJ=J+K;I=I+1;endwhileL=MLE=2^L;R=1;LE1=LE/2;W=1;W1=exp(-j*pi/LE1);whileR=LE1P=R;whileP=NQ=P+LE1;T=A(Q)*W;A(Q)=A(P)-T;A(P)=A(P)+T;P=P+LE;endW=W*W1;R=R+1;endL=L+1;endend附录二主程序closeall;clearalls_f=15;s_xy=15;s_t=20;s_mar=10;w_f=2.5;w_axi=2;%信号产生T0=8e-2;Ts=1/800;N=T0/Ts;t=0:Ts:(N-1)*Ts;f1=50;f2=100;f3=200;fp=80;fs=120;Sig=4*sin(f1*2*pi*t)+4*sin(f2*2*pi*t)+4*sin(f3*2*pi*t);s_need=4*sin(f1*2*pi*t);s_1=4*sin(f2*2*pi*t);s_2=4*sin(f3*2*pi*t);P_qian=fftlzr(Sig);figure(1)h1=stem(2*pi/(N*Ts).*(0:N-1),abs(P_qian));gridonhx1=xlabel('模拟角频率/rad/s');ht1=title('滤波前幅度频谱');ha1=gca;axistight%滤波器设计Fs=1/Ts;wp=2/Ts*tan(fp*2*pi*Ts/2);ws=2/Ts*tan(fs*2*pi*Ts/2);rp=1;rs=15;[Nl,Wn]=buttord(wp,ws,rp,rs,'s');%求阶数及3dB截止频率[z,p,k]=buttap(Nl);%求极零点增益[b,a]=zp2tf(z,p,k);%系统传递函数[b2,a2]=lp2lp(b,a,Wn);%转换为所需低通滤波器[HS,WS]=freqs(b2,a2,linspace(0,2000,512));%求0到2000rad/s的频率响应[hs,wss]=freqs(b2,a2,[wp,ws,[f1,f2,f3].*2*pi]);%求通带阻带边缘角频率及原始信号频率响应%画模拟滤波器频谱并标注figure(6)subplot(211)hsf=plot(WS,20*log10(abs(HS)));holdonstem(wss,20*log10(abs(hs)),'filled','r:','linewidth',w_f)text(wss(1)-100,20*log10(abs(hs(1)))-10,...{'\Omega_p=160\pirad/s';[num2str(20*log10(abs(hs(1)))),'dB']});text(wss(2)-65,20*log10(abs(hs(2)))-12,...{'\Omega_s=240\pirad/s';[num2str(20*log10(abs(hs(2)))),'dB']});text(wss(3)-100,20*log10(abs(hs(3)))-10,...{'\Omega_1=100\pirad/s';[num2str(20*log10(abs(hs(3)))),'dB']});text(wss(4)-65,20*log10(abs(hs(4)))-12,...{'\Omega_2=200\pirad/s';[num2str(20*log10(abs(hs(4)))),'dB']});text(wss(5)-70,20*log10(abs(hs(5)))-10,...{'\Omega_3=400\pirad/s';[num2str(20*log10(abs(hs(5)))),'dB']});hxsf=xlabel('模拟角频率/rad/s');hysf=ylabel('幅值/dB');htsf=title('幅频特性');hasf=gca;gridonsubplot(212)hsx=plot(WS,angle(HS));holdonstem(wss,angle(hs),'filled','r:','linewidth',w_f)text(wss(1),angle(hs(1))+0.7,...{'\Omega_p=160\pirad/s';[num2str(angle(hs(1))),'rad']});text(wss(2),angle(hs(2))+0.7,...{'\Omega_s=240\pirad/s';[num2str(angle(hs(2))),'rad']});text(wss(3)-150,an
本文标题:信号分析与处理实验报告
链接地址:https://www.777doc.com/doc-2691186 .html