《数字信号处理》实验指导.docx

1、数字信号处理实验指导数字信号处理实验指导 康文静一、Matlab与信号处理有关的基本函数1. 单位采样序列zeros函数单位采样序列，或称为单位脉冲序列(编程序常用Delta表示)： 如果在时间轴上延迟了k个单位，得到，即： Matlab函数zeros(1,N)可产生一个包含N个零的行向量。下图为产生delta函数的几种编程方法。2. 单位阶跃序列(Unit step sequence)函数ones 单位阶跃序列： 在MATLAB中可以利用ones( )函数实现， 函数ones(1,N)可产生一个包含N个1的行向量。下图为产生单位阶跃序列的几种编程方法。3. 矩形序列rectangle seq

2、uence矩形序列： 编程： rect=zeros(1,3),ones(1,4),zeros(1,3)rect = 0 0 0 1 1 1 1 0 0 04. 实指数序列stem函数，幂的表示实指数序列： 用matlab实现，并画出相应图形。 n=0:8; x=(0.6).n; stem(n,x);5. 正弦序列函数plot，sin，cos；常量pi用matlab实现如下序列，并画出相应图形。 n=0:0.1:10; x=2*sin(0.6*pi*n)+cos(0.1*pi*n+pi/4); plot(n,x);%对应第一幅图 stem(n,x);%对应第二幅图6. 复指数序列exp函数复指数

3、序列： 在MATLAB中实现： 用matlab实现。 n=1:6; x=exp(j*0.5*n)x = 0.8776 + 0.4794i 0.5403 + 0.8415i 0.0707 + 0.9975i -0.4161 + 0.9093i -0.8011 + 0.5985i -0.9900 + 0.1411i7. 序列的翻转fliplr函数，用matlab实现，并画出相应图形。 n=0:0.1:10; x=2*sin(0.6*pi*n)+cos(0.1*pi*n+pi/4); y=fliplr(x); n=-fliplr(n); plot(n,y); stem(n,y);8. 计算序列的线性

4、卷积函数conv已知，程序如下： x=1,1,2,3; h=2,1,8,6,3 y=conv(x,h)y = 2 3 13 22 28 39 24 99. 序列的傅里叶变换（DTFT）-freqz函数freqz函数计算序列的离散时间傅里叶变换在给定的离散频率点上的抽样值。假设：freqz的调用方式如下：(1)H,w=freqz(b,a,N),其中b和a分别是分子和分母多项式的系数向量。此函数在单位圆上半部0,pi上等间隔地计算N个频率响应。返回该系统的N点频率响应矢量w和N点复数频率响应矢量H，如果N没有说明，缺省值为512。(2)H=freqz(b,a,w),返回矢量w指定的那些频率点上的频

5、率响应，频率范围为0,pi。(3)H,w=freqz(b,a,N，whole),在整个单位圆上等间隔地计算N个频率响应，即频率范围为0,2pi。编程题目：因果系统，利用freqz函数试画出的幅频相应和相频相应。利用函数abs/angle/real/image等计算DTFT的幅度、相位以及实部和虚部。程序如下： b=0.5;a=1,-0.85; H,w=freqz(b,a,200,whole); magH=abs(H(1:101);%计算幅度 phaH=angle(H(1:101);%计算相位 w=w(1:101); %画图 subplot(2,1,1); plot(w/pi,magH); gr

6、id;%加网格 xlabel(frequency Unit:pi); ylabel(Magnitude); title(Magnitude Response); subplot(2,1,2); plot(w/pi,phaH/pi); grid; xlabel(frequency Unit:pi); ylabel(Phase Unit:pi); title(Phase Response);10. Z变换函数tf2zp/ zp2tf/zplane假设：利用函数z,p,k=tf2zp(b,a)可确定分子分母多项式按z的降幂排列的有理z变换式的零点z/极点p/增益常数k。函数b,a= zp2tf (z

7、,p,k)用来实现相反的过程。函数zplane可以用来画出z变换的零极点图，有如下两种调用方式：zplane（zeros，poles）zplane（b，a）编程题目：已知离散系统差分方程为：，求其系统函数，画出零/极点示意图，并判断系统的稳定性。由差分方程可得程序如下： b=1,-0.5,0; a=1,0.75,0.125; z,p,k=tf2zp(b,a); zplane(b,a) title(Pole-Zero Plot)由上图可知，系统的极点全部在单位圆内，因此系统是稳定的。11. 求解差分方程fliter函数filter函数求解给定输入x(n)的差分方程的解，该函数调用形式为：y=fi

8、lter(b,a,x)编程题目：已知离散系统差分方程为：,计算并画出冲击响应h(n)(n=-10,-9,0,1,50);由h(n)确定系统是否稳定。程序如下： b=1;a=1,-1,0.5; %求单位脉冲响应 n=-10:50;x=(n)=0; h=filter(b,a,x) stem(n,h) xlabel(n); ylabel(h(n); title(Impulse Response); axis(-10,50,-1,1.5) sum(abs(h) %12. 计算信号的DFT和IDFT函数fft和ifft调用形式如下： fft(x)函数用来计算M点DFT，其中M是序列x的长度； fft(x

9、,N)计算N点DFT，N是用户指定的长度。若序列x的长度MN，则将序列截短为N点序列，再作N点DFT；若序列x的长度MN，则将序列截短为N点序列，再作N点IDFT；若序列x的长度M n=0:15; x=sin(0.25*pi*n);%序列x(n) dft_16=fft(x);%16点DFT dft_256=fft(x,256);%256点DFT L=0:255; plot(L/256,abs(dft_256); xlabel(normalized frequency); ylabel(Magnitude); title(256点DFT); hold on; plot(n/16,abs(dft_

10、16),*); title(256点DFT和16点DFT（*）);编程题目：已知，为16点的有限长序列，计算其16点和256点DFT。程序如下： n=0:15; x=sin(0.25*pi*n);%序列x(n) dft_16=fft(x);%16点DFT dft_256=fft(x,256);%256点DFT L=0:255; subplot(2,1,1);stem(L/256,abs(dft_256); xlabel(normalized frequency); ylabel(Magnitude); title(256点DFT); subplot(2,1,2);stem(n/16,abs(d

11、ft_16),*); xlabel(normalized frequency); ylabel(Magnitude); title(16点DFT（*）);二、基于Matlab的数字滤波器设计1. IIR滤波器设计函数(1)巴特沃斯模拟滤波器N，Wn=buttord(Wp,Ws,Rp,Rs, s)此函数可获得巴特沃斯滤波器的参数N和Wn，分别是在给定通带边界频率Wp、阻带边界频率Ws、通带最大衰减Rp和阻带最小衰减Rs的条件下，所需要的Butterworth滤波器的最小阶数和对应的3dB带宽，这里Wp和Ws的单位是rad/s。z,p,k=buttap(N)求出N阶巴特沃斯低通滤波器的零点、极点和

12、增益。其中p是长度为N的列向量，分布在单位圆的左半平面，而零点z是空矩阵。B,A=butter(N,Wn,s)获得巴特沃斯低通滤波器系统函数的分子多项式B和分母多项式A的系数。例：设计一个模拟巴特沃斯低通滤波器满足以下指标：通带截止频率为fp =1Hz，通带最大衰减1dB（p=1）；阻带截止频率fs=4Hz，阻带最小衰减30dB（s=30）。 %Butterworth低通滤波器设计 Wp=2*pi*1;Ws=2*pi*4; Rp=3;Rs=30; %设计滤波器 N,Wn=buttord(Wp,Ws,Rp,Rs,s)N = 3Wn = 7.9490 z,p,k=buttap(N)z = p =

13、-0.5000 + 0.8660i -0.5000 - 0.8660i -1.0000 k = 1.0000 B,A=butter(N,Wn,s)B = 0 0 0 502.2695A = 1.0000 15.8980 126.3731 502.2695因此要满足所要求的指标需要3阶巴特沃斯滤波器，其系统函数为： 画图Butterworth低通滤波器f1=linspace(0,Wp,5) f2=linspace(Wp,Ws,15)f3=linspace(Ws,2*pi*10,30) h1=20*log10(abs(freqs(B,A,f1) h2=20*log10(abs(freqs(B,A,

14、f2) h3=20*log10(abs(freqs(B,A,f3) plot(f1 f2 f3/(2*pi),h1,h2,h3) grid xlabel(Frequency in Hz) ylabel(Gain in dB)此外，chebwin, cheb1ord, cheb2ord, cheby1, cheby2, ellip, ellipord等函数可分别设计切比雪夫和椭圆滤波器。(2)脉冲响应不变法的Matlab实现bz,az=impinvar(b,a,Fs)此函数实现用脉冲响应不变法将模拟滤波器转换为数字滤波器。其中b和a分别是模拟滤波器的系统函数H(s)的分子多项式和分母多项式的系数

15、，Fs是脉冲响应不变法的采样频率，单位为Hz。如果Fs没有说明，其缺省值为1Hz。运算结果bz、az分别表示数字滤波器的系统函数H(z)的分子多项式和分母多项式的系数。例：利用impinvar函数将模拟系统函数Ha(s)转变为数字系统函数H(z)，采样周期为Ts=0.5。其中。 b=0,3;a=1,4,3;T=0.5;Fs=1/T; bz,az=impinvar(b,a,Fs)bz = 0 0.2876az = 1.0000 -0.8297 0.1353因此，数字滤波器为(3)双线性变换法的Matlab实现zd,pd,kd=bilinear(z,p,k,Fs)其中z，p，k分别表示模拟滤波器系

16、统函数的零点、极点和增益，Fs是采样频率，zd,pd,kd是对应的数字滤波器系统函数的零点、极点和增益。numd，dend=bilinear(num，den，Fs)其中num，den和Fs分别表示模拟滤波器系统函数分子、分母多项式系数，Fs是采样频率；numd，dend是对应的数字滤波器系统函数的分子、分母多项式系数。例：利用双线性变换法设计低通巴特沃斯数字滤波器，其技术指标为：Wp=0.2*pi,Rp=1dB;Ws=0.4*pi,Rs=15dB,T=1s. Wp=0.2*pi;Rp=1;Ws=0.4*pi;Rs=15;T=1, Fs=1/T omegap=(2/T)*tan(Wp/2);%通

17、带截止频率预畸 omegas=(2/T)*tan(Ws/2);%阻带截止频率预畸 N,Wn=buttord(omegap,omegas,Rp,Rs,s) B,A=butter(N,Wn,s) b,a=bilinear(B,A,Fs)%画图 h,w=freqz(b,a,256);h1=20*log10(abs(h);%增益(dB)plot(w/pi,h1); gridxlabel(Digital Frequency in pi units)ylabel(Gain in dB)axis(0 1 -80 10)2. FIR滤波器设计函数(窗函数法)(1)常用的窗函数（a） 矩形窗调用格式 ：w=bo

18、xcar(n)，根据长度n产生一个矩形窗w。（b） 三角窗调用格式 ：w=triang(n)，根据长度n产生一个三角窗w。（c） 汉宁窗调用格式 ：w=hanning(n)，根据长度n产生一个汉宁窗w。（d） 哈明窗调用格式 ：w=hamming(n)，根据长度n产生一个哈明窗w。（e） 布莱克曼窗调用格式 ：w=blackman(n)，根据长度n产生一个布莱克曼窗w。（f） 凯塞窗调用格式 ：w=kaiser(n，beta)，根据长度n产生一个布莱克曼窗w，beta是控制窗形状的参数。上述窗函数返回的变量w是一个长度为n的列向量，表示窗函数在这n点的取值。(2)窗函数法设计FIR滤波器(a)

19、函数fir1调用格式 ：b=fir1(n,wn,ftype,Window)其中，n为阶数，因此h(n)长度为n+1。wn是截止频率，其取值在01之间，这是以采样频率为基准频率的标称值，故1对应采样频率；若Wn是标量，则用来设计低通滤波器；如果要设计带通滤波器，输入wn=w1,w2为一矢量。ftype是滤波器的类型，低通时省略该参数，高通为high，带阻为stop。Window表示选用的窗函数类型，以列向量形式表示，向量Window的长度必须为n+1，若缺省，则fir1默认使用哈明窗。B对应设计好的滤波器的单位取样响应h(n)。(a)函数fir2调用格式 ：b=fir2(n,f,m,Window

20、)该函数采用窗函数法设计具有任意频率响应的FIR数字滤波器。其中f是频率向量，其值在01(标称值)之间，1对应采样频率，其中第一个点必须是0，最后一个点必须是1，而且频率点必须是递增的。m对应于频率点f处期望的幅频响应，f和m的长度必须相等。例：设计长度为8的线性相位FIR滤波器。其理想幅频特性满足Window=boxcar(8); %长度为8的矩形窗Window b=fir1(7,0.4,Window);%用矩形窗作为脉冲响应的窗函数 freqz(b,1)；%矩形窗对应的频率响应freqz (B,A,.) with no output arguments plots the magnitud

21、e and unwrapped phase of the filter in the current figure window. window=blackman(8); %长度为8的布拉克曼窗Window b=fir1(7,0.4,Window); freqz(b,1)例：设计线性相位带通滤波器，阶数N=15，上下边带截至频率分别为0.3*pi，0.5*pi。 window=blackman(16); b=fir1(15,0.3 0.5,window); freqz(b,1)例： 用矩形窗和哈明窗设计线性相位FIR低通滤波器，wc=0.25*pi ，N=10。 N=10; M=128; b1=fir1(N,0.25,boxcar(N+1); b2=fir1(N,0.25,hamming(N+1); h1=freqz(b1,1,M); h2=freqz(b2,1,M); f=0:0.5/M:0.5-0.5/M; plot(f,abs(h1),-.,f,abs(h2);%画出幅频响应 legend(矩形窗,哈明窗) grid二、建立一个M-file 的方法 (1) 建立一个M-file(2)编写M-file，实现阶跃序列程序。(3)保存