吉林大学数字信号处理实验报告.docx

1、吉林大学数字信号处理实验报告数字信号处理课程设计实验报告基础实验实验一 离散时间系统及离散卷积一、 实验目的（ 1）熟悉 MATLAB 软件的使用方法。（2）熟悉系统函数的零极点分布、单位脉冲响应和系统频率响应等概念。（ 3）利用 MATLAB 绘制系统函数的零极点分布图、系统频率响应和单位脉冲 响应。（ 4）熟悉离散卷积的概念，并利用 MATLAB 计算离散卷积。二、实验容1、离散时间系统的单位脉冲响应（1）选择一个离散时间系统；（2）用笔进行差分方程的递推计算；（ 3）编制差分方程的递推计算程序；（4）在计算机上实现递推运算； （5）将程序计算结果与笔算的计算结果进行比较，验证程序运行的正

2、确性；2.离散系统的幅频、相频的分析方法（1）给定一个系统的差分方程或单位取样响应； （2）用笔计算几个特殊的幅频、相频的值，画出示意曲线图；（ 3）编制离散系统的幅频、相频的分析程序； （4）在计算机上进行离散系统的幅频、相频特性，并画出曲线； （5）通过比较，验证程序的正确性；3.离散卷积的计算（1）选择两个有限长序列，用笔计算其线性卷积；（ 2）编制有限长序列线性卷积程序；（ 3）利用计算程序对（ 1）选择的有限长序列进行卷积运算； （4）比较结果验证程序的正确性。三、实验要求（ 1）编制实验程序，并给编制程序加注释； （2）按照实验容项要求完成笔算结果；（3）验证编制程序的正确性，记录

3、实验结果。（4）至少要求一个除参考实例以外的实例，在实验报告中，要描述清楚实例中的系统，并对实验结果进行解释说明。四、实验程序及其结果1、离散时间系统的脉冲响应% y(n)-0.9y(n-1)+ 0.5y(n-2)=x(n)+0.5x(n-1) 设 x(n)=2n; 0=n=15 clear all;b=1,0.5;a=1,-0.9,0.5;n=0:15; x=impseq(0,0,15)h=filter(b,a,x); % 系统冲击响应 subplot(2,1,1);stem(n,h);title( 系统冲激响应 ); xlabel(n);ylabel(h);n=0:15; x=2.*n;

4、nx=0:15;nh=0:15;y=conv_m(x,nx,h,nh);ny=length(y); n=0:ny-1;subplot(2,1,2);stem(n,y);title( 系统对 x(n) 响应 ); xlabel(n);ylabel(y);B=roots(b);A=roots(a);figure; zplane(B,A);2、离散系统的幅频、相频的分析方法% 差分方程为 y(n)-1.76y(n-1)+1.1829y(n-2)-0.2781y(n-3)% =0.0181x(n)+0.0543x(n-1)+0.0543x(n-2)+0.0181x(n-3) b=0.0181,0.05

5、43,0.0543,0.0181;a=1.000,-1.76,1.1829,-0.2781;m=0:length(b)-1;l=0:length(a)-1;K=500;k=0:1:K;w=pi*k/K;H=(b*exp(-j*m*w)./(a*exp(-j*l*w); % 计算频率相应 magH=abs(H); % magH 为幅度angH=angle(H); %a ngH 为相位subplot(2,1,1); plot(w/pi,magH);grid;xlabel(以pi 为单位的频率 ); ylabel(幅度 );title( 幅度响应 );subplot(2,1,2); plot(w/p

6、i,angH); grid; xlabel(以pi 为单位的频率 ); ylabel(相位); title( 相位响应 );3、离散卷积的计算% x=1,4,3,5,3,6,5 , -4=n=2% h=3,2,4,1,5,3, -2=n=0).*(n=4); % 三角波序列 subplot(3,1,1);stem(n,x1);title( 三角波序列 );xlabel(n);N=8; k=0:N-1;X1=fft(x1,N); magX1=abs(X1);phaX1=angle(X1); subplot(3,1,2);stem(k,magX1);xlabel(k);ylabel( 三角波 DF

7、T 的幅度 ); subplot(3,1,3);stem(k,phaX1);xlabel(k);ylabel( 三角波 DFT 的相位 );(2)反三角波clear all; n=0:1:7;x2=(4-n).*(n=0).*(n=4); % 反三角波序列subplot(3,1,1);stem(n,x2);title( 反三角波序列 );xlabel(n);N=8; k=0:N-1;X2=fft(x2,N); magX2=abs(X2);phaX2=angle(X2); subplot(3,1,2);stem(k,magX2); xlabel(k);ylabel( 幅度 );subplot(3

8、,1,3);stem(k,phaX2); xlabel(k);ylabel( 相位 );结果分析： 由图知：三角波和反三角波序列的波形不同，当 N=8 时，正反三角波的幅频特性相同，因为两者的时域只差一个相位3、余弦信号的 FFT 分析clea all;F=50; N=32;T=0.000625; % (1) F=50,N=32,T=0.000625 n=1:N;x=cos(2*pi*F*n*T);figure(1); subplot(2,1,1); plot(n,x);ylabel(x(n); xlabel(n);title(1) F=50,N=32,T=0.000625);X=fft(x)

9、; magX=abs(X); subplot(2,1,2);plot(n,X); ylabel(FFT|X|); xlabel(f(pi);F=50; N=32;T=0.005; % (2) F=50,N=32,T=0.005 n=1:N;x=cos(2*pi*F*n*T); figure(2);subplot(2,1,1);plot(n,x);ylabel(x(n);xlabel(n);title(2) F=50,N=32,T=0.005);X=fft(x); magX=abs(X);subplot(2,1,2); plot(n,X); ylabel(FFT|X|); xlabel(f(pi

10、);F=50,N=32,T=0.0046875; % (3) F=50,N=32,T=0.0046875n=1:N;x=cos(2*pi*F*n*T);figure(3);subplot(2,1,1); plot(n,x); ylabel(x(n);xlabel(n);title(3) F=50,N=32,T=0.0046875);X=fft(x);magX=abs(X);subplot(2,1,2);plot(n,X); ylabel(FFT|X|); xlabel(f(pi);F=50,N=32,T=0.004; % (4) F=50,N=32,T=0.004n=1:N; x=cos(2*

11、pi*F*n*T);figure(4);subplot(2,1,1);plot(n,x); ylabel(x(n);xlabel(n);title(4) F=50,N=32,T=0.004);X=fft(x); magX=abs(X);subplot(2,1,2);plot(n,X);ylabel(FFT|X|); xlabel(f(pi);F=50,N=64,T=0.000625; % (5) F=50,N=64,T=0.000625 n=1:N;x=cos(2*pi*F*n*T);figure(5); subplot(2,1,1); plot(n,x); ylabel(x(n); xlab

12、el(n);title(5) F=50,N=64,T=0.000625);X=fft(x); magX=abs(X);subplot(2,1,2); plot(n,X); ylabel(FFT|X|); xlabel(f(pi);iFFTIX_910OE900 0=L=N,09=J (90-O90(!d)j90-OFFTUl吕2535FFTO 5O15 20 25 30 35 Sx(n)结果分析； 不同的采样间隔会产生不同的栅栏效应， 相当于透过栅栏 观赏风景 ,只能看到频谱的一部分 ,而其它频率点看不见 ,因此很可能 使一部分有用的频率成分被漏掉， 从而产生不同的频谱图， 减小栅栏 效应可用

13、提高采样间隔也就是频率分辨力的方法来解决。 间隔小， 频 率分辨力高，被 “挡住 ”或丢失的频率成分就会越少。但会增加采样点 数，使计算工作量增加。4、卷积计算clear all;n=0:14;x1=1.*(n=0).*(n=0).*(n=9); % 第二组函数n=0:19;h2=(0.8).*sin(0.5.*n);N=27; y=circonvt(x2,h2,N);disp(y); subplot(3,1,2);i=0:N-1; stem(i,y); ylabel(式( 2)卷积 );n=0:9; x3=(1-0.1.*n); h3=0.1.*n; % 第三组函数N=17; y=circo

14、nvt(x3,h3,N);disp(y); subplot(3,1,3);i=0:N-1; stem(i,y); ylabel(式( 3)卷积 );实验三 IIR 数字滤波器设计一、试验目的1、学习模拟数字变换滤波器的设计方法；2.掌握双线性变换数字滤波器设计方法；3.掌握实现数字滤波器的具体方法。二、实验容1、设计一个巴特沃思数字低通滤波器，设计指标如下： 通带 wws=0.35pi 幅度衰减不小于 15dB；2、编制计算设计的数字滤波器幅度特性和相位特性的程序， 并进行 实验验证。3、编制实现该数字滤波器程序并且实现数字滤波(1)分别让满足所设计的滤波器的通带、过渡带、阻带频率特性 的正弦

15、波通过滤波器，验证滤波器性能；(2)改变正弦抽样时间，验证数字低通滤波器的模拟截止频率实 抽样时间的函数。三、实验要求1、编制实验容要求的程序，并给程序加注释；2、根据实验结果，给出自己设计的数字滤波器的幅度特性和相位特 性；3.用所设计的滤波器对不同频率的正弦波信号进 行滤波，以说明其特性；4、 fp=0.2KHz,Rp=1dB,fs=0.3KHz,As=25dB,T=1ms ；分别用脉冲响应 不变法及双线性变换法设计一 Butterworth 数字低通滤波器，观察所 设计数字滤波器的幅频特性曲线， 记录带宽和衰减量， 检查是否满足 要求。比较这两种方法的优缺点。四、实验程序及结果分析1、巴

16、特沃斯低通滤波器clear all;wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;R=10(-Rp/20);A=10(-As/20);b,a=afd_buttap(wp,ws,Rp,As); C,B,A = sdir2cas(b,a); % 将直接型装换成级联型db,mag,pha,w=freqs_m(b,a,0.5*pi); % 计算系统频率响应ha,x,t=impulse(b,a);subplot(2,1,1);plot(w/pi,mag);title( 幅度响应 );grid;ylabel(幅度);xlabel(以pi 为单位的频率 );subplot(2,1,2);pl

17、ot(w/pi,pha);title( 相位响应 );grid;ylabel(相位 );xlabel( 以 pi 为单位的频率 );2、脉冲响应不变法设计巴特沃斯滤波器% 冲击响应不变法设计数字巴特沃斯低通滤波器 clear all;T=0.0001; wp=200*2*pi*T;Rp=1; ws=300*2*pi*T;As=25; % 数字滤波器指标T=1;omegap=wp/T; omegas=ws/T; %转换为模拟域指标 cs,ds=afd_buttap(omegap,omegas,Rp,As); % 模拟巴特沃斯滤波器的计算 % 冲击响应不变法b,a=imp_invr(cs,ds,T

18、);C,B,A=dir2par(b,a); % 将 IIR 滤波器直接型转换为并联型的函数 % 计算数字滤波器频率响应db,mag,pha,grd,w=freqz_m(b,a); subplot(3,1,1);plot(w/pi,mag);title( 冲击响应不变法设计巴特沃斯滤波器 ); axis(0 0.5 0 1);grid;ylabel(|幅度 |);xlabel(以pi 为单位的频率 ); subplot(3,1,2);plot(w/pi,db);axis(0 0.8 -400 0);grid;ylabel(对数幅度( db) );xlabel( 以pi 为单位的频率 );3、双线

19、性变换法设计巴特沃斯滤波器clear all;T=0.0001;Fs=1/T; wp=200*2*pi*T;Rp=1; ws=300*2*pi*T;As=25; % 数字滤波器指标 % 转换成模拟域指标 omegap=(2/T)*tan(wp/2);omegas=(2/T)*tan(ws/2); %模拟巴特沃斯滤波器的计算 cs,ds=afd_buttap(omegap,omegas,Rp,As);%双线性变换b,a=bilinear(cs,ds,Fs);C,B,A=sdir2cas(b,a);db,mag,pha,grd,w=freqz_m(b,a);subplot(3,1,1);plot(

20、w/pi,mag);ylabel(| 幅度 |);title( 双线性变换法设计巴特沃斯滤波器 ); xlabel(以pi 为单位的频率 );axis(0 0.5 0 1);grid;subplot(3,1,2);plot(w/pi,db);axis(0 0.8 -400 0);grid; ylabel(对数幅度( db) );xlabel( 以pi 为单位的频率 ); subplot(3,1,3);plot(w/pi,pha);axis(0 0.4 -4 4);grid; ylabel(相位 );xlabel( 以 pi 为单位的频率 );4、数字滤波clear all;F=0.1;N1=5

21、12;T1=0.1;n=0:5:N1;x=sin(2*pi*F*n*T1); wp=0.2*pi;ws=0.35*pi;Rp=1;As=15; T=1;Fs=1/T;OmegaP=(2/T)*tan(wp/2); %预畸变OmegaS=(2/T)*tan(ws/2); ep=sqrt(10(Rp/10)-1); % 求平方根 Ripple=sqrt(1/(1+ep.2);Attn=1/10(As/20);N=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(OmegaP/OmegaS);OmegaC=OmegaP/(10.(Rp/10)-1).(1/

22、(2*N); cs,ds=u_buttap(N,OmegaC);b,a=bilinear(cs,ds,Fs); %括号为模拟滤波器系数，中括号中的为数字滤波器系数 mag,db,pha,w=freqz_m(b,a);y=filter(b,a,x);subplot(2,1,1);stem(x);axis(0 100 -1 1);title( 输入正弦序列 );subplot(2,1,2);stem(y); axis(0 100 -1 1);title( 滤波输出序列 );结果分析：滤波输出序列前有一段是由于延时造成的实验四 FIR 数字滤波器1、FIR 滤波器的设计clear all;% 用哈明

23、窗函数法来设计 FIR 滤波器 选 N=51 wp=0.2*pi; ws=0.35*pi; N=51;n=0:1:N-1;wc=(ws+wp)/2; % 理想低通的截止频率 hd=ideal_lp(wc,N); % 理想低通的冲激响应 w_ham=(hamming(N);h=hd.*w_ham % FIR 滤波器的冲激响应 db,mag,pha,grd,w=freqz_m(h,1); delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1); % 实际的通带衰减 As=-round(max(db(ws/delta_w+1:1:501); % 实际的阻带衰减

24、 subplot(2,2,1);stem(n,hd);title( 理想冲激响应 );axis(0 N-1 -0.1 0.3);xlabel(n); ylabel(hd(n); subplot(2,2,2);stem(n,h);title( 实际冲激响应 );axis(0 N-1 -0.1 0.3);xlabel(n); ylabel(h(n); subplot(2,2,3);plot(w/pi,pha);title( 相位响应 );grid;xlabel(以pi 为单位的频率 ); ylabel(相位); subplot(2,2,4);plot(w/pi,db);title( 幅度响应( db) );grid; xlabel(以pi 为单位的频率 ); ylabel(H(db);2、分别用矩形窗、汉宁窗、海明窗、三角窗和 Blackman 窗设计一个线性相位的带通滤波器。clear all;w1=0.3*pi; w2=0.5*pi; N=51;n=0:1:N-1;hd=ideal_lp(w2,N)-ideal_lp(w1,N); % 理想低通的冲激响应 w_rec=(rect