信号分析与处理实验报告.docx
- 文档编号:24793753
- 上传时间:2023-06-01
- 格式:DOCX
- 页数:22
- 大小:213.26KB
信号分析与处理实验报告.docx
《信号分析与处理实验报告.docx》由会员分享,可在线阅读,更多相关《信号分析与处理实验报告.docx(22页珍藏版)》请在冰豆网上搜索。
信号分析与处理实验报告
华北电力大学
实验报告
|
|
实验名称FFT的软件实现实验(Matlab)
IIR数字滤波器的设计
课程名称信号分析与处理
|
|
专业班级:
电气化1308学生姓名:
袁拉麻加
学号:
2成绩:
指导教师:
杨光实验日期:
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.码位倒置程序流程图
k Y k=k+1 N 按新序号生产输入序列B(过程略) c.蝶形运算程序流程图 五、程序代码与实验结果 a.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); else disp('输入正确,请运行下一步') end %% %码位倒置% fork=0: n-1 forj=1: m%取M位的二进制数% x1(j)=bitget(k,j);%倒取出二进制数% end x1=num2str(x1);%将数字序列转化为字符串% y(k+1)=bin2dec(x1);%二进制序列转化为十进制数% clearx1 end fork=1: n B(k)=A(y(k)+1);%时间抽取序列% end clearA %% %计算% 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; end W=W*W1; end end B%输出X(k)% %% 验证结果: 例6-4 b.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); else disp('输入正确,请运行下一步') end %% %码位倒置% fork=0: n-1 forj=1: m%取M位的二进制数% x1(j)=bitget(k,j);%倒取出二进制数% end x1=num2str(x1);%将数字序列转化为字符串% y(k+1)=bin2dec(x1);%二进制序列转化为十进制数% clearx1 end fork=1: n B(k)=A(y(k)+1);%时间抽取序列% end clearA %% %计算% 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; end W=W*W1; end end B=conj(B);%取共轭% B=B/n%输出x(n)% 验证结果: 六、实验心得与结论 本次实验借助于Matlab软件,我避开了用C平台进行复杂的复数运算,在一定程度上简化了程序,并添加了简单的检错代码,码位倒置我通过查阅资料,使用了一些函数,涉及到十-二进制转换,数字-文本转换,二-文本转换,相对较复杂,蝶运算我参考了书上了流程图,做些许改动就能直接实现。 通过本次实验,我对FFT的过程更加熟悉了,尤其是WN系数的特性,可以通过累乘来实现W系数的更新。 IIR数字滤波器的设计实验 一、实验目的 1.掌握实现模拟信号数字化的编程方法。 2.掌握实现时间抽取快速傅里叶变换(FFT)的编程方法。 3.掌握设计IIR数字滤波器的双线性变换法。 4.加深对信号分析与处理的理解。 二、实验内容 1.生成三个不同频率的正弦信号 、 、 ,将它们叠加成一个信号 。 用满足采样定理的采样频率对模拟信号 进行采样,得到离散时间序列 。 2.编制子程序实现FFT运算,并对 进行FFT运算,观察它的频谱。 3.设计满足已知的设计指标的IIR数字滤波器,并用它对 进行滤波处理。 对滤波后得到的波形再做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)~=M A=[A(: );zeros(2^(fix(M)+1)-N,1)]; end N=length(A);M=log2(N); I=0;J=0;L=1; whileI~=N-1 ifI t=A(J+1); A(J+1)=A(I+1); A(I+1)=t; end K=N/2; whileK<=J J=J-K; K=K/2; end J=J+K; I=I+1; end whileL<=M LE=2^L;R=1; LE1=LE/2; W=1;W1=exp(-j*pi/LE1); whileR<=LE1 P=R; whileP<=N Q=P+LE1;T=A(Q)*W; A(Q)=A(P)-T; A(P)=A(P)+T; P=P+LE; end W=W*W1; R=R+1; end L=L+1; end end 附录二主程序 closeall;clearall s_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));gridon hx1=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))); holdon stem(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;gridon subplot(212) hsx=plot(WS,angle(HS)); holdon stem(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,angle(hs(3))-0.8,... {'\Omega_1=100\pirad/s';[num2str(angle(hs(3))),'rad']}); text(wss(4)-20,angle(hs(4))+0.8,... {'\Omega_2=200\pirad/s';[num2str(angle(hs(4))),'rad']}); text(wss(5)-50,angle(hs(5))-1.3,... {'\Omega_3=400\pirad/s';[num2str(angle(hs(5))),'rad']}); hxsx=xlabel('模拟角频率/rad/s');hysx=ylabel('相值/rad'); htsx=title('相频特性');hasx=gca;gridon [b3,a3]=bilinear(b2,a2,Fs);%双线性变换到数字滤波器 [H,W]=freqz(b3,a3,linspace(0,0.5*pi,512));%求数字滤波器频率响应 [hd,wd]=freqz(b3,a3,[wp,ws,[f1,f2,f3].*2*pi]*Ts);%求通带阻带边缘数字角频率及原始信号数字频率响应 wd=wd/pi; %画数字滤波器频谱并标注 figure (2) subplot(211) h2=plot(W/pi,20*log10(abs(H))); holdon stem(wd,20*log10(abs(hd)),'filled','r: ','linewidth',w_f) text(wd (1)-0.0035,20*log10(abs(hd (1)))-12,... {['\omega_p=',num2str(wd (1)),'\pirad'];[num2str(20*log10(abs(hd (1)))),'dB']}); text(wd (2)-0.003,20*log10(abs(hd (2)))-12,... {['\omega_s=',num2str(wd (2)),'\pirad'];[num2str(20*log10(abs(hd (2)))),'dB']}); text(wd(3)-0.0035,20*log10(abs(hd(3)))-12,... {['\omega_1=',num2str(wd(3)),'\pirad'];[num2str(20*log10(abs(hd(3)))),'dB']}); text(wd(4)-0.003,20*log10(abs(hd(4)))-12,... {['\omega_2=',num2str(wd(4)),'\pirad'];[num2str(20*log10(abs(hd(4)))),'dB']}); text(wd(5)-0.003,20*log10(abs(hd(5)))-12,... {['\omega_3=',num2str(wd(5)),'\pirad'];[num2str(20*log10(abs(hd(5)))),'dB']}); hx2=xlabel('归一化数字角频率/\pirad'); hy1=ylabel('幅值/dB'); ht2=title('幅频特性'); ha2=gca; gridon subplot(212) h3=plot(W/pi,angle(H));grid;title('相频特性'); holdon stem(wd,angle(hd),'filled','r: ','linewidth',w_f) text(wd (1),angle(hd (1))+1,... {['\omega_p=',num2str(wd (1)),'\pirad'];[num2str(angle(hd (1))),'\pirad']}); text(wd (2)-0.002,angle(hd (2))+1,... {['\omega_p=',num2str(wd (2)),'\pirad'];[num2str(angle(hd (2))),'\pirad']}); text(wd(3)-0.003,angle(hd(3))-1,... {['\omega_p=',num2str(wd(3)),'\pirad'];[num2str(angle(hd(3))),'\pirad']}); text(wd(4)-0.003,angle(hd(4))-2,... {['\omega_p=',num2str(wd(4)),'\pirad'];[num2str(angle(hd(4))),'\pirad']}); text(wd(5)-0.003,angle(hd(5))-1.3,... {['\omega_p=',num2str(wd(5)),'\pirad'];[num2str(angle(hd(5))),'\pirad']}); hx3=xlabel('归一化数字角频率/\pirad'); hy2=ylabel('相值/\pirad'); ht3=title('相频特性'); ha3=gca; gridon %滤波 data=filter(b3,a3,Sig); figure(3) h4=plot(t,data,'k',t,Sig,'b',t,s_need,'r'); ht4=title('滤波效果'); xlabel('t/s','fontsize',s_xy); ylabel('幅值','fontsize',s_xy); legend('滤波后波形','滤波前波形','50Hz信号') ha4=gca; gridon P_hou=fftlzr(data); figure(4) subplot(211) h5=stem(2*pi/(N*Ts).*(0: N-1),abs(P_qian)); hx4=xlabel('模拟角频率/rad/s'); ht5=title('滤波前幅度频谱'); axis([0,1.5e3,0,2e3]) text(f1*2*pi-20,abs(P_qian(f1*N*Ts+1))+200,[num2str(f1)'Hz']); text(f2*2*pi-20,abs(P_qian(f2*N*Ts+1))+200,[num2str(f2)'Hz']); text(f3*2*pi-20,abs(P_qian(f3*N*Ts+1))+200,[num2str(f3)'Hz']); ha5=gca;gridon subplot(212) h6=stem(2*pi/(N*Ts).*(0: N-1),abs(P_hou)); hx5=xlabel('模拟角频率/rad/s'); ht6=title('滤波后幅度频谱'); axis([0,1.5e3,0,2e3]) text(f1*2*pi-20,abs(P_hou(f1*N*Ts+1))+200,[num2str(f1)'Hz']); text(f2*2*pi-20,abs(P_hou(f2*N*Ts+1))+200,[num2str(f2)'Hz']); text(f3*2*pi-20,abs(P_hou(f3*N*Ts+1))+200,[num2str(f3)'Hz']); ha6=gca;gridon figure(5) h7=plot(t,Sig,t,s_need,'r',t,s_1,'g',t,s_2,'c'); hx6=xlabel('t/s'); hy3=xlabel('幅值'); ht7=title('生成信号'); set(h7,'linewidth',2) legend('合成信号','50Hz','100Hz','200Hz'); ha7=gca;gridon %图像处理 h=[h1,h2,h3,h4',h5,h6,hsf,hsx]; ht=[ht1,ht2,ht3,ht4,ht5,ht6,ht7,htsf,htsx]; hx=[hx1,hx2,hx3,hx4,hx5,hx6,hxsf,hxsx]; hy=[hy1,hy2,hy3,hysf,hysx];ha=[ha1,ha2,ha3,ha4,ha5,ha6,ha7,hasf,hasx]; set(h([2: 6,9: 10]),'Linewidth',w_f) set(ha,'Linewidth',w_axi,'Fontsize',s_f); set([hxhyht],'Fontsize',s_xy,'FontWeight','bold'); set((1: 6),'outerposition',get(0,'screensize'),'color','w')
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 信号 分析 处理 实验 报告