数字信号处理实验讲义.docx

1、数字信号处理实验讲义王实验一 离散信号的matlab实现一、 实验目的1、 熟悉matlab软件，学会matlab语言的编写2、 使用matlab软件产生一些常见的离散信号3、 掌握用matlab软件作信号的相关分析二、 实验环境计算机操作系统、matlab软件三、 实验内容1、用matlab程序产生下列离散信号或连续信号，并画出其波形。 a单位抽样序列的产生参考程序：N=100;x=zeros(1,N); 产生一个1行N列值全为0的矩阵，如看成数组x（1）-x（100）都为0x(1)=1;n=0:N-1;stem(n,x);产生序列参考程序：N=100;x=zeros(1,N);k=20;x

2、(k+1)=1;xn=0:N-1;stem(xn,x);b单位阶跃序列的产生参考程序：N=32;x=ones(1,N);产生一个1行N列值全为1的矩阵n=0:N-1;stem(n,x);产生序列参考程序：N=32;k=20;x1=zeros(1,k);x2=ones(1,N-k);x=x1,x2;xn=0:N-1;stem(xn,x);c模拟信号，以t=0.01n（n=0：N-1）进行采样后的离散信号。参考程序： N=128；n=0:N-1; t=0.01*n; x=2*sin(4*pi*t)+5*cos(8*pi*t); figure(1); subplot(211); stem(t,x);

3、 subplot(212); stem(n,x); d产生一个sinc(t)=sint/t抽样函数参考程序：n=200;step=4*pi/n;t=-2*pi:step:2*pi;y=sinc(t);plot(t,y,t,zeros(size(t);%同时画出y(t）和横轴grid on;plot(t,y,t,zeros(size(t),zeros(size(y),y);%同时画出y(t）和横轴、纵轴 e方波信号square(t) square(t,duty) 产生周期是2pi，幅度为正负1的方波，duty占空比,高电平跟整个周期的比值参考程序：t=0:0.01:2*pi;y=square(t

4、,3);plot(t,y); 试产生一个周期为1，高低电平分别为半个周期的方波信号2、相关分析去除噪声x(n)=sin(2*pi*n)+u(n) 噪声为高斯分布白噪声，使用相关分析去除噪声，噪声1功率为1，噪声2功率为0.1%rxy=xcorr(x,y);%rx=xcorr(x,Mlag,flag) Mlag表示rx的单边长度，总长度为2Mlag+1, flag-biased rx(m)/N -unbiased rx(m)/(N-abs(m)参考程序：N=500;p1=1;p2=0.1;f=1/8;Mlag=60;u=randn(1,N);u2=u*sqrt(p2);n=0:N-1;s=sin

5、(2*pi*f*n);x1=u(1:N)+s;rx1=xcorr(x1,Mlag,biased);subplot(211);plot(-Mlag:Mlag,rx1);x2=u2(1:N)+s;rx2=xcorr(x2,Mlag,biased);subplot(212);plot(-Mlag:Mlag,rx2);3.书后习题1.15实验二、离散信号的傅里叶变换一、 实验目的1、 进一步熟悉matlab软件的使用，熟悉matlab的编程语言2、 用matlab语言编写程序进行离散信号的傅里叶分析二、 实验原理 设离散序列，长度为N，其DTFT定义为： 在实际计算中无法取到无限长序列，通常通过无限长

6、序列加窗作有限长序列的DTFT。三、 实验内容试求序列，其中n=0，1，N-1，N=12的DTFT。并画出其幅频特性曲线和相频特性曲线。 若N=24,36,120,其幅频特性曲线和相频特性曲线如何变化，为什么？参考程序：N=12;n=0:1:N-1;xn=cos(n*pi/6);w=0:0.01:2*pi;M=length(w);xejw=zeros(1,M);for i=0:1:N-1 xejw=xejw+cos(i*pi/6)*exp(-j*w*i);end;xr=abs(xejw);xphase=angle(xejw);subplot(211);plot(w,xr);subplot(21

7、2);plot(w,xphase);实验三、FFT频谱分析及应用一.实验目的1、通过实验，加深对FFT的理解，熟悉FFT子程序2、熟悉应用FFT对典型信号进行频谱分析的方法二、实验原理与方法 在各种信号序列中，有限长序列占有重要地位，对有限长序列，可以利用离散傅里叶变换（DFT）进行分析。DFT不但可以很好地反映序列的频谱特性，而且易于用快速算法（FFT）在计算机上实现。 设离散序列，长度为N，其DFT定义为：k=0，1，2,N-1有限长序列的DFT是其Z变换在单位圆上的等距采样，因此可以用于序列的谱分析。FFT是DFT的一种快速算法，它是对变换式进行一次次分解，使其成为若干小点数的组合，从而

8、减少运算量。常用的FFT是以2为基数的，其长度。它的效率高，程序简单，使用方便。当要变换的序列长度不等于2的整数次方时，为了使用以2为基数的FFT，可以用末尾补零的方法，使其长度延长至2的整数次方。在MATLAB中，可以用函数FFT来实现。格式：y=FFT(x) y=FFT(x,n)说明：y=FFT(x)为利用FFT算法计算矢量的离散傅里叶变换，当x为矩阵时，y为矩阵x每一列的FFT。当x的长度为2的整数次方时，则FFT函数采用基2的FFT算法，否则采用稍慢的混合基算法。y=FFT(x,n)采用n点FFT。当x的长度小于n时，FFT函数在x的尾部补零，以构成n点数据；当x的长度大于n时，FFT

9、会截断序列x。三、实验内容1、模拟信号，以进行采样。求（1）N=40点FFT的幅度频谱，从图中，能否观察出信号的2个频率分量？ （2）提高采样点数，如N=128,再求该信号的幅度频谱，此时幅度频谱发生了什么变化？信号的2个模拟频率和数字频率各为多少？FFT频谱分析结果与理论上是否一致？2、研究高密度频谱与高分辨率频谱 设有连续信号以采样频率对信号采样，分析下列三种情况的幅频特性。（1） 采样数据长度N=16点，做N=16的FFT，并画出幅频特性（2） 采集数据长度N=16点，补零到256点，做256点的FFT，并画出幅频特性（3） 采样数据长度N=256点，做N=256的FFT，并画出幅频特性

10、。参考程序：N=40;n=0:N-1;t=0.01*n;x=2*sin(4*pi)+5*cos(8*pi*t);k=0:1:N/2;w=2*pi/N*k;X=fft(x,N);magx=abs(X(1:1:N/2+1);subplot(211);stem(n,x);title(signal x(n);subplot(212);plot(w/pi,magx);title(fft N=40);xlabel(频率（单位：pi）);ylabel(| X|);grid;取N=128时,程序参考以上程序.实验内容2的程序参考以上程序.实验四、IIR数字滤波器的设计一实验目的1、掌握脉冲响应不变法和双线性变

11、换法设计IIR数字滤波器的具体方法和原理,熟悉双线性变换法和脉冲响应不变法设计低通,带通IIR数字滤波器的计算机编程.2、观察双线性变换法和脉冲响应不变法设计的数字滤波器的频域特性,了解双线性变换法和脉冲响应不变法的特点和区别.3、熟悉Butterworth滤波器和Chebyshev滤波器的频率特性.二实验原理与方法 本实验利用模拟滤波器设计IIR数字滤波器,这是IIR数字滤波器设计最常用的方法.利用模拟滤波器设计,需要将模拟域的转换为数字域的,最常用的转换方法为脉冲响应不变法和双线性变换法.1、脉冲响应不变法 用数字滤波器的单位脉冲响应序列模仿模拟滤波器的冲激响应,让正好等于的采样值,即:

12、其中T为采样间隔,如果以及分别表示的拉氏变换及的z变换,则: 在MATLAB中，可用函数impinvar实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为： b,a=impinvar(c,d,fs) b,a=impinvar(c,d) 其中，c、d分别为模拟滤波器的分子和分母多项式系数向量，fs为采样频率（Hz），缺省值fs=1Hz，b、a分别为数字滤波器分子和分母多项式系数向量。例：已知，T=0.1，利用脉冲响应不变法求H(z)MATLAB程序：c=2 1;d=1 2 3;T=0.1;b,a=impinvar(c,d,1/T); 执行结果：b=0.2 -0.1881 a=1 -1.7

13、916 0.8187 数字滤波器为： 2、 双线性变换法 S平面与z平面之间满足下列映射关系： S平面的虚轴映射在z平面的单位圆上，s平面的左半平面完全映射到z平面的单位圆内。双线性变换不存在频率混叠问题。 在MATLAB中，提供了一个叫做bilinear的函数实现这种映射，调用格式为： b,a=bilinear(c,d,fs) 双线性变换是一种非线性变换，即，这种非线性引起的幅频特性畸变可通过预畸得到校正。3、 一般设计步骤（1） 给定技术指标转换为模拟低通原型设计性能指标。（2） 估计满足性能指标的模拟低通原型阶数和截止频率。利用MATLAB中buttord、cheb1ord、cheb2o

14、rd、ellipord等函数，调用格式如：n,wn=buttord(wp,ws,rp,rs,s)（3）设计模拟低通原型 利用MATLAB中buttap、cheb1ap、cheb2ap、ellipap等函数，调用格式如： z,p,k=buttap(n) 采用上述函数所得到原型滤波器的传递函数为零点、极点、增益形式，需要和函数 c,d=zp2tf(z,p,k)配合使用，以转化为多项式形式。（4）由模拟低通原型经频率变换获得模拟低通、高通、带通或带阻滤波器。 利用MATLAB中lp2lp、lp2hp、lp2bp、lp2bs等函数，调用格式如： c1,d1=lp2lp(c,d,wn)（5）利用脉冲响应

15、不变法或者双线性不变法，实现模拟滤波器到数字滤波器的映射。说明：MATLAB信号处理工具箱还提供了模拟滤波器设计的完全工具函数：butter、cheby1、cheby2、ellip、besself。用户只需一次调用就可自动完成以上设计步骤中的3-4步，调用格式如： c,d=butter(n,wn,ftype,s)三、实验内容1、已知fp=0.3kHz，rp=1.2db，fs=0.2kHz,rs=20db,T=1ms,利用双线性变换法设计一个chebyshevI型数字高通滤波器，观察通带损耗和阻带衰减是否满足要求。 按照要求写出matlab 程序，并附上所设计的滤波器传递函数H(z)及相应的幅频

16、特性曲线。2、要求fp=0.2kHz，rp=1db，fs=0.3kHz,rs=25db，T=1ms，分别用脉冲响应不变法和双线性变换法设计一个butterworth数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。程序1：rp=1.5;rs=40;T=0.001;fp=0.3;fs=0.2;wp=2*pi*fp;ws=2*pi*fs;%wp,ws为数字频率wp1=(2/T)*tan(wp/2);ws1=(2/T)*tan(ws/2);%预畸为模拟频率n,wn=cheb1ord(wp1,ws1,rp,rs,s); %求阶数n，3db截

17、止频率wnz,p,k=cheb1ap(n,rp)； %模拟低通原型c,d=zp2tf(z,p,k); %转化为多项式形式c1,d1=lp2hp(c,d,wn); %模拟低通到模拟高通b,a=bilinear(c1,d1,1/T); %双线性变换法设计出数字滤波器db,mag,pha,grd,w=freqz_m(b,a)； %作数字滤波器的幅频响应和相频响应subplot(211);plot(w/pi,mag);title(幅频特性);xlabel(w(/pi);ylabel(|H(jw)|);subplot(212);plot(w/pi,pha/pi);title(相频特性);xlabel(w

18、(/pi);ylabel(pha(/pi);freqz_m函数如下：functiondb,mag,pha,grd,w=freqz_m(b,a)%滤波器的幅值响应(相对、绝对)%Usage:db,mag,pha,grd,w=freqz_m(b,a) %w 采样频率； b系统函数H(z)的分子项(对FIR，b=h)%a 系统函数H(z)的分母项(对FIR，a=1)H,w=freqz(b,a); %500点的复频响应mag=abs(H); %绝对幅值响应db=20*log10(mag/max(mag); %相对幅值响应pha=angle(H); %相位响应grd=grpdelay(b,a,w); %

19、群延迟响应程序2：fp=0.2;fs=0.3;wp=2*pi*fp;ws=2*pi*fs;%数字频率rp=1;rs=25;T=0.001;wp1=tan(wp/2);ws1=tan(ws/2);%预变形为模拟频率n,wn=buttord(wp1,ws1,rp,rs,s);%设计模拟滤波器阶数c,d=butter(n,wn,s);bz,az=bilinear(c,d,1/2); %双线性变换法 由s到zh,w=freqz(bz,az,1024,1000);plot(w,abs(h);grid on;程序3：fp=200;fs=300;rp=1;rs=25;T=0.001;wp=2*pi*fp;w

20、s=2*pi*fs;% 模拟频率n,wn=buttord(wp,ws,rp,rs,s);%设计模拟滤波器的阶数z,p,k=buttap(n);bp,ap=zp2tf(z,p,k);bs,as=lp2lp(bp,ap,wp);bz,az=impinvar(bs,as,1000);%bz,az=butter(n,wn/pi);%设计数字滤波器db,mag,pha,grd,w=freqz_m(b,a);%作出数字滤波器的频率特性图subplot(211);plot(w/pi,mag);title(幅频特性);xlabel(w(/pi);ylabel(|H(jw)|);subplot(212);plo

21、t(w/pi,pha/pi);title(相频特性);xlabel(w(/pi);ylabel(pha(/pi);实验五、FIR数字滤波器的设计一、实验目的：1、掌握用窗函数法设计FIR滤波器的原理和方法，熟悉相应的计算机编程。2、了解用频率采样法设计FIR滤波器的原理和实现。3、熟悉线性相位FIR滤波器的幅频特性和相频特性。4、了解不同窗函数对滤波器性能的影响。二、实验原理与方法1、窗函数法：窗函数法设计线性相位FIR滤波器步骤：（1）确定数字滤波器的性能要求：临界频率wk，滤波器单位脉冲响应长度N。（2）根据性能要求，合理选择h(n)的奇、偶对称性，从而确定理想频率响应的幅频特性和相位特性

22、。（3）求理想单位脉冲响应。在实际计算中，可对按M等间隔采样，并对其求IDFT得，用代替。（4）选择合适的窗函数w(n),由求所设计的FIR滤波器单位脉冲响应。（5）求，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。 窗函数的傅氏变换的主瓣决定了的过渡带宽。的旁瓣大小和多少决定了在通带和阻带范围内的波动幅度。2.频率采样法 频率采样法是从频域出发，将给定的加以等间隔采样：然后以此作为实际数字滤波器的频率特性的采样值，即令由H(k)通过IDFT可得有限长序列h(n)代入到z变换中可得：用频率采样法设计线性相位FIR滤波器的一般步骤为：（1）由设计

23、要求选择滤波器的种类（2）根据线性相位的约束条件，确定和，进而得到（3）将代入内插公式得到所设计滤波器的频率响应。关于第3步，在MATLAB中可由函数h=real(ifft(H,N)和db,mag,pha,grd,w=freqz_m(h,l)实现。三、实验内容：1、用窗函数法设计一线性相位FIR低通滤波器，设计指标为：（1）选择一个合适的窗函数，取N=15，确定脉冲响应，并给出所设计的滤波器的频率响应图，分析是否满足设计要求。（2）若取N=45，重复这一设计，观察幅频和相位特性变化，分析长度N变化的影响。（3）保持N=45不变，改变窗函数，（如由hamming窗变为blackman窗），观察并

24、记录窗函数对滤波器幅频特性的影响，比较两种窗的特点。2、用kaiser窗设计一个数字带通滤波器，设计指标为：低阻带： 低通带： 高通带： 高阻带： 3、 用频率采样法设计一个低通滤波器：问：（1）采样点数N=33，过渡带设置1个采样点，H(k)=0.5，最小阻带衰减为多少，是否满足设计要求？（2）采样点数N=34，过渡带设置2个采样点，,最小阻带衰减为多少，是否满足设计要求？参考程序1：N=15;M=128;b1=fir1(N,0.15,boxcar(N+1);b2=fir1(N,0.15,hamming(N+1);h1=freqz(b1,1,M);h2=freqz(b2,1,M);t=0:N

25、;subplot(221);stem(t,b2,.);grid;f=0:0.5/M:0.5-0.5/M;subplot(222);plot(f,abs(h1),b-,f,abs(h2),g-);grid;db1=20*log10(abs(h1)/1);db2=20*log10(abs(h2)/1);subplot(224);subplot(224);plot(f,db1,b-,f,db2,g-);参考程序2：N=15;M=128;b2=fir1(N,0.35,0.65,kaiser(N+1);h2=freqz(b2,1,M);t=0:N;subplot(221);stem(t,b2,.);gr

26、id;f=0:0.5/M:0.5-0.5/M;subplot(222);plot(f*2,abs(h2),g-);grid;db2=20*log10(abs(h2)/1);subplot(224);subplot(224);plot(f*2,db2,g-);参考程序3：%过渡带设置1个采样点N=33;alpha=(N-1)/2;k=0:N-1;wk=(2*pi/N)*k;Hk=ones(1,4),0.5,zeros(1,24),0.5,ones(1,3);angH=-alpha*(2*pi)/N*k;H=Hk.*exp(j*angH);h=real(ifft(H,N);hejw,w=freqz(h,N);subplot(211);plot(w/pi,abs(hejw);%过渡带设置2个采样点N=34;alpha=(N-1)/2;k=0:N-1;wk=(2*pi/N)*k;Hk=ones(1,4),0.5925,0.1099,zeros(1,23),-0.1099,-0.5925,-ones(1,3);angH=-alpha*(2*pi)/N*k;H=Hk.*exp(j*angH);h=real(ifft(H,N);hejw,w=freqz(h,N);subplot(212);plot(w/pi,abs(hejw);