信号分析与处理实验报告

华北电力大学实验报告|实验名称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+1Yjm?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)h6=stem(2*pi/(N*Ts).*(0:N-1),abs(P_hou);hx5=xla