信号分析与处理实验报告.docx

华北电力大学实 验 报 告| 实验名称 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+1 y jm? nx1=num2str(x1);y(k+1)=bin2dec(x1);clear x1k=k+1ykn?n按新序号生产输入序列b（过程略）c. 蝶形运算程序流程图五、程序代码与实验结果a.fft程序：%clear all;close all;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%码位倒置%for k=0:n-1 for j=1:m %取m位的二进制数% x1(j)=bitget(k,j); %倒取出二进制数% end x1=num2str(x1); %将数字序列转化为字符串% y(k+1)=bin2dec(x1); %二进制序列转化为十进制数% clear x1endfor k=1:n b(k)=a(y(k)+1); %时间抽取序列%endclear a%计算%for l=1:m %分解为m级进行运算% le=2l; %第l级群间隔为2l% le1=2(l-1); %第l级中共有2(l-1)个wn乘数，进行运算蝶运算的两数序号相隔le1% w=1; w1=exp(-1i*pi/le1); for r=1:le1 %针对第r个wn系数进行一轮蝶运算，共进行le1次% for p=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; endendb %输出x(k)%验证结果：例6-4b.ifft程序：%clear all;close all;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%码位倒置%for k=0:n-1 for j=1:m %取m位的二进制数% x1(j)=bitget(k,j); %倒取出二进制数% end x1=num2str(x1); %将数字序列转化为字符串% y(k+1)=bin2dec(x1); %二进制序列转化为十进制数% clear x1endfor k=1:n b(k)=a(y(k)+1); %时间抽取序列%endclear a%计算%for l=1:m %分解为m级进行运算% le=2l; %第l级群间隔为2l% le1=2(l-1); %第l级中共有2(l-1)个wn乘数，进行运算蝶运算的两数序号相隔le1% w=1; w1=exp(-1i*pi/le1); for r=1:le1 %针对第r个wn系数进行一轮蝶运算，共进行le1次% for p=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; endendb=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);if fix(m)=m a=a(:);zeros(2(fix(m)+1)-n,1);endn=length(a);m=log2(n);i=0;j=0;l=1;while i=n-1 if ij t=a(j+1); a(j+1)=a(i+1); a(i+1)=t; end k=n/2; while k=j j=j-k; k=k/2; end j=j+k; i=i+1;endwhile l=m le=2l; r=1; le1=le/2; w=1;w1=exp(-j*pi/le1); while r=le1 p=r; while p=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;endend附录二 主程序close all;clear alls_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);grid onhx1=xlabel(模拟角频率/ rad/s);ht1=title(滤波前幅度频谱);ha1=gca;axis tight%滤波器设计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);hold onstem(wss,20*log10(abs(hs),filled,r:,linewidth,w_f)text(wss(1)-100,20*log10(abs(hs(1)-10,. omega_p=160pi rad/s;num2str(20*log10(abs(hs(1),db);text(wss(2)-65,20*log10(abs(hs(2)-12,. omega_s=240pi rad/s;num2str(20*log10(abs(hs(2),db);text(wss(3)-100,20*log10(abs(hs(3)-10,. omega_1=100pi rad/s;num2str(20*log10(abs(hs(3),db);text(wss(4)-65,20*log10(abs(hs(4)-12,. omega_2=200pi rad/s;num2str(20*log10(abs(hs(4),db);text(wss(5)-70,20*log10(abs(hs(5)-10,. omega_3=400pi rad/s;num2str(20*log10(abs(hs(5),db);hxsf=xlabel(模拟角频率/ rad/s);hysf=ylabel(幅值/db);htsf=title(幅频特性);hasf=gca;grid onsubplot(212)hsx=plot(ws,angle(hs);hold onstem(wss,angle(hs),filled,r:,linewidth,w_f)text(wss(1),angle(hs(1)+0.7,. omega_p=160pi rad/s;num2str(angle(hs(1),rad);text(wss(2),angle(hs(2)+0.7,. omega_s=240pi rad/s;num2str(angle(hs(2),rad);text(wss(3)-150,angle(hs(3)-0.8,. omega_1=100pi rad/s;num2str(angle(hs(3),rad);text(wss(4)-20,angle(hs(4)+0.8,. omega_2=200pi rad/s;num2str(angle(hs(4),rad);text(wss(5)-50,angle(hs(5)-1.3,. omega_3=400pi rad/s;num2str(angle(hs(5),rad);hxsx=xlabel(模拟角频率/ rad/s);hysx=ylabel(相值/ rad);htsx=title(相频特性);hasx=gca;grid onb3,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);hold onstem(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),pi rad;num2str(20*log10(abs(hd(1),db);text(wd(2)-0.003,20*log10(abs(hd(2)-12,. omega_s=,num2str(wd(2),pi rad;num2str(20*log10(abs(hd(2),db);text(wd(3)-0.0035,20*log10(abs(hd(3)-12,. omega_1=,num2str(wd(3),pi rad;num2str(20*log10(abs(hd(3),db);text(wd(4)-0.003,20*log10(abs(hd(4)-12,. omega_2=,num2str(wd(4),pi rad;num2str(20*log10(abs(hd(4),db);text(wd(5)-0.003,20*log10(abs(hd(5)-12,. omega_3=,num2str(wd(5),pi rad;num2str(20*log10(abs(hd(5),db);hx2=xlabel(归一化数字角频率/ pi rad);hy1=ylabel(幅值/db);ht2=title(幅频特性);ha2=gca;grid onsubplot(212)h3=plot(w/pi,angle(h);grid;title(相频特性);hold onstem(wd,angle(hd),filled,r:,linewidth,w_f)text(wd(1),angle(hd(1)+1,. omega_p=,num2str(wd(1),pi rad;num2str(angle(hd(1),pi rad);text(wd(2)-0.002,angle(hd(2)+1,. omega_p=,num2str(wd(2),pi rad;num2str(angle(hd(2),pi rad);text(wd(3)-0.003,angle(hd(3)-1,. omega_p=,num2str(wd(3),pi rad;num2str(angle(hd(3),pi rad);text(wd(4)-0.003,angle(hd(4)-2,. omega_p=,num2str(wd(4),pi rad;num2str(angle(hd(4),pi rad);text(wd(5)-0.003,angle(hd(5)-1.3,. omega_p=,num2str(wd(5),pi rad;num2str(angle(hd(5),pi rad);hx3=xlabel(归一化数字角频率/ pi rad);hy2=ylabel(相值/ pi rad);ht3=title(相频特性);ha3=gca;grid on%滤波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;grid onp_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;grid onsubplot(212