数字信号处理实验指导书.doc

数字信号处理实验指导书Matlab在数字信号处理中的应用聂小燕 编著电子科技大学成都学院2012年3月数字信号处理上机实验本课程上机实验8课时，共4个实验。学时实验内容备注2实验1 离散时间信号与系统 2实验2 FFT及其应用2实验3 IIR滤波器的设计2实验4 FIR滤波器的设计实验一 离散时间信号与系统【实验目的】熟悉用matlab分析离散时间信号与系统【实验内容】1用MATLAB实现函数impseq(n0,n1,n2)，使函数实现。（1）该函数的MATLAB程序为function x,n=impseq(n0,n1,n2)n=n1:n2; x =(n-n0)=0; （2）函数调用举例% 2, -5=n=0;（2）函数调用举例% x(n) = nu(n)-u(n-10)+10*exp(-0.3(n-10)(u(n-10)-u(n-20); 0=n=20n = 0:20;x1 = n.*(stepseq(0,0,20)-stepseq(10,0,20);x2 = 10*exp(-0.3*(n-10).*(stepseq(10,0,20)-stepseq(20,0,20);x = x1+x2;stem(n,x);title(调用单位脉冲序列实例)xlabel(n);ylabel(x(n);axis(0,20,-1,11)3用MATLAB实现常见的离散时间序列。（1）参考程序：n=0:10; x=(0.8).n; stem(n,x,ro)（2）x(n) = exp(2+3j)*n)n=-10:10; x=exp(2+3j)*n);stem(n,x,ro)title(复指数序列)xlabel(n);ylabel(x(n);（3）x(n) = cos(0.04*pi*n) + 0.2*w(n)n = 0:50;x = cos(0.04*pi*n)+0.2*randn(size(n);stem(n,x);title(正弦序列随机序列)xlabel(n);ylabel(x(n);axis(0,50,-1.4,1.4)（4）x(n) = .,5,4,3,2,1,5,4,3,2,1,.; -10=n=9 序列的周期化n=-10:9;x=5,4,3,2,1;xtilde=x * ones(1,4);xtilde=(xtilde(:);stem(n,xtilde);title(有限长序列的周期化)xlabel(n);ylabel(xtilde(n);axis(-10,9,-1,6)4计算两个序列的卷积。教材P39页（7）a=8 -2 -1 2 3b=2 3 -1 -3c=conv(a,b)M=length(c)-1n=0:1:Mstem(n,c)xlabel(n)ylabel(幅度);5计算幅频响应与相频响应.教材P39页（6）b = 1 -sqrt(2) 1; % H(Z)分子多项式系数a = 1 -0.67 0.9; % H(Z)分母多项式系数h,w = freqz(b,a); % 计算频率响应am = 20*log10(abs(h); % 计算幅频响应并绘图subplot(2,1,1);plot(w,am);ph = angle(h); % 计算相频响应并绘图subplot(2,1,2);plot(w,ph);6计算离散系统的单位脉冲响应.教材P39页（8）N=50;a=1 -2;b=1 0.1 -0.06;x=1 zeros(1,N-1);k=0:1:N-1;y=filter(a,b,x);stem(k,y)xlabel(n);ylabel(幅度)21实验二 快速傅立叶变换FFT及其应用【实验目的】1. 熟悉matlab中的有关函数2. 应用FFT对典型信号进行频谱分析，并了解可能出现的问题，以便在实际中正确应用FFT3. 应用FFT实现序列的线性卷积。【实验内容】（1） 高斯序列的时域和幅频特性。改变参数P，q，运行程序exam2_1.m（2） 正弦衰减序列的时域和幅频特性。改变参数a，f，运行程序exam2_2.m（3） 三角波序列的时域和幅频特性。（4） 连续时间信号采样截短后的频谱分析。（5） 用FFT分别计算16点的循环卷积和线性卷积。实验三 无限长单位脉冲响应(IIR)滤波器的设计【实验目的】1、掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理和方法；2、观察双线性变换法和脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法和脉冲响应不变法的特点和区别。【实验内容】使用MATLAB编写程序，实现IIR数字滤波器的设计。涉及脉冲响应不变法和双线性变换法设计IIR数字滤波器的方法、不同设计方法得到的IIR滤波器频域特性异同等知识点。【实验原理与方法和手段】1、脉冲响应不变法所谓脉冲响应不变法就是使数字滤波器的单位脉冲响应序列h(n)等于模拟滤波器的单位冲激响应和(t)的采样值，即：，其中，T为采样周期。在MATLAB中，可用函数impinvar实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为：num,den=impinvar(b,a,fs)num,den=impinvar(b,a)其中，b,a分别为模拟滤波器的分子和分母多项式系数向量；fs为采样频率（Hz），缺省值fs=1Hz；num,den分别为数字滤波器分子和分母多项式系数向量。2、双线性变换法：由于s平面和z平面的单值双线性映射关系为s，其中T为采样周期。因此，若已知模拟滤波器的传递函数，将上式代入即可得到数字滤波器的系统函数H(z)。在双线性变换中，模拟角频率和数字角频率的变换关系为： 可见，和w之间的变换关系为非线性的。在MATLAB中，可用函数bilinear实现从模拟滤波器到数字滤波器的双线性变换映射，调用格式为：num,den=bilinear(b,a,fs)【常用函数介绍】1、Matlab信号处理工具箱中提供了设计巴特沃思模拟滤波器的函数buttord、buttap和butter，格式如下：（1） 用于计算巴特沃思模拟低通滤波器的阶N和3dB截止频率Wc (即本书中的符号)。其中，Wp 和Ws分别是滤波器的通带截止频率和阻止截止频率，单位为rad/s；Rp和Rs分别是通带最大衰减系数和阻带最小衰减系数，单位为dB。（2）用于计算N阶巴特沃思归一化（=1）模拟低通滤波器系统函数的零、极点和增益因子，返回长度为N的向量z和p分别给出N个零点和极点，G是滤波器增益。得到的滤波器系统函数形式如下：其中，和分别是向量z和p的第k个元素。如果要从零、极点得到系统函数的分子和分母多项式系数向量B和A，可以调用结构转换函数（3），结构转换后系统函数的形式为其中，M是向量B的长度，N是向量A的长度，分别是向量B和A的第k个元素。（3） 用于计算巴特沃思模拟滤波器系统函数中分子和分母多项式系数向量B和A，其中N和分别是滤波器的阶和3dB截止频率，返回向量B和A中的元素和分别是上面的表示式中的分母和分子系数。ftype缺省时，设计低通滤波器;ftype=high时，设计高通滤波器；ftype=stop时，设计带阻滤波器，此时为向量，且ftype缺省时，设计带通滤波器，带通的频率区间为。S缺省时，设计数字滤波器。例如：设计一个满足下列指标要求的巴特沃思模拟低通滤波器。指标Matlab程序：运行结果如下：N=6Wc=2.2496B= 0 0 0 0 0 0 129.5917A=1.0000 8.6916; 37.7720; 104.0667;191.1447 ;222.5973 ; 129.5917即：幅频响应曲线如图所示。2、Matlab中设计数字滤波器的函数都是采用双线性变换法，将模拟滤波器转换为数字滤波器。这些函数及其凋用格式如下（巴特沃思数字滤波器）：（1）该格式用于计算巴特沃思数字滤波器的阶N和3dB截止频率的归一化值（关于归一化）。调用参数Wp 和Ws分别是数字滤波器的通带截止频率和阻带截止频率的归一化值（关于归一化），要求和，其中1表示数字频率（对应模拟频率，为采样频率）。和分别是通带最大衰减和阻带最小衰减，单位dB。当时，设计高通滤波器；当和是二元向量时，设计带通（）或带阻（）滤波器，这时返回参数也是二元向量。（2）该格式用于计算巴特沃思模拟滤波器系统函数中分子和分母多项式系数向量B和A。调用参数N 和WC分别是巴特沃思数字滤波器的阶和3dB截止频率归一化值（关于归一化）。当ftype缺省时，设计低通滤波器；当ftype=high时，设计高通滤波器；ftype=stop时，设计带阻滤波器，此时为二元向量，和分别是带阻滤波器的通带3dB下，上截止频率的归一化值（关于归一化）；为向量，且ftype缺省时，设计带通滤波器，带通的频率区间为。注意设计出的带通和带阻数字滤波器是2N阶的，这是因为带通滤波器可表示为一个N阶低通滤波器与一个N阶高通滤波器的级联。由函数的返回向量B和A可写出数字滤波器的系统函数为：其中，分别是向量B和A的第k个元素。例如：用双线性变换法设计巴特沃思数字低通滤波器，指标要求通带截止频率，通带最大衰减；阻带截止频率，阻带最小衰减解 设计步骤如下：（1）给出数字滤波器的指标。（2）计算巴特沃思数字滤波器的阶N和3dB截止频率。；（3）用双线性变换法设计巴特沃思数字低通滤波器。；运行结果如下：=0.0335 0.1006 0.1006 0.0335=1.000 -1.4245 0.8827 -1.1900设计出的数字滤波器系统函数为程序如下：例：利用脉冲响应不变法，把系统函数为的模拟滤波器变换成等价的数字滤波器，采样间隔手算： Matlab:B=1,1;A=1,5,6;T=0.1;Fs=1/T;Bz,Az=impinvar(B,A,Fs);%用脉冲响应不变法将模拟滤波器变换成数字滤波器运行结果：Bz=0.1000 -0.0897Az=1.000 -1.5595 0.60651、 设采样周期T=250s（采样频率fs =4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB边界频率为fc =1kHz。 B,A=butter(3,2*pi*1000,s); num1,den1=impinvar(B,A,4000); h1,w=freqz(num1,den1); B,A=butter(3,2/0.00025,s); num2,den2=bilinear(B,A,4000); h2,w=freqz(num2,den2); f=w/pi*2000; plot(f,abs(h1),-.,f,abs(h2),-); grid; xlabel(频率/Hz ) ylabel(幅值/dB) 程序中第一个butter的边界频率21000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图3.1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。同时，也看到双线性变换法，在z=-1即=或f=2000Hz处有一个三阶传输零点，这个三阶零点正是模拟滤波器在=处的三阶传输零点通过映射形成的。 2、 设计一数字高通滤波器，它的通带为400500Hz，通带内容许有0.5dB的波动，阻带内衰减在小于317Hz的频带内至少为19dB，采样频率为1,000Hz。 wc=2*1000*tan(2*pi*400/(2*1000); wt=2*1000*tan(2*pi*317/(2*1000); N,wn=cheb1ord(wc,wt,0.5,19,s); B,A=cheby1(N,0.5,wn,high,s); num,den=bilinear(B,A,1000); h,w=freqz(num,den); f=w/pi*500; plot(f,20*log10(abs(h); axis(0,500,-80,10); grid; xlabel() ylabel(幅度/dB) 图3.2给出了MATLAB计算的结果，可以看到模拟滤波器在=处的三阶零点通过高通变换后出现在=0（z=1）处，这正是高通滤波器所希望得到的。 3、 设计一巴特沃兹带通滤波器，其dB边界频率分别为f2=110kHz和f1=90kHz，在阻带f3 = 120kHz处的最小衰减大于dB，采样频率fs=400kHz。w1=2*400*tan(2*pi*90/(2*400); w2=2*400*tan(2*pi*110/(2*400); wr=2*400*tan(2*pi*120/(2*400); N,wn=buttord(w1 w2,0 wr,3,10,s); B,A=butter(N,wn,s); num,den=bilinear(B,A,400); h,w=freqz(num,den); f=w/pi*200; plot(f,20*log10(abs(h); axis(40,160,-30,10); grid; xlabel(频率/kHz) ylabel(幅度/dB) 图3.3给出了MATLAB计算的结果，可以看出数字滤波器将无穷远点的二阶零点映射为z=1的二阶零点，数字带通滤波器的极点数是模拟低通滤波器的极点数的两倍。 4、 一数字滤波器采样频率fs = 1kHz，要求滤除100Hz的干扰，其dB的边界频率为95Hz和105Hz，原型归一化低通滤波器为 w1=95/500; w2=105/500; B,A=butter(1,w1, w2,stop); h,w=freqz(B,A); f=w/pi*500; plot(f,20*log10(abs(h); axis(50,150,-30,10); grid; xlabel(频率/Hz) ylabel(幅度/dB) 图3.4为MATLAB的计算结果 【实验报告要求】简述实验目的和实验原理，用笔算求出设计的滤波器，并用计算机计算结果验证，总结实验中的主要结论，写出收获和体会。实验四 FIR数字滤波器的设计【实验目的】（1）掌握用窗函数法，频率采样法及优化设计法设计FIR滤波器的原理及方法，熟悉相应的MATLAB编程； （2）熟悉线性相位FIR滤波器的幅频特性和相频特性； （3）了解各种不同窗函数对滤波器性能的影响。【实验原理与方法】(一) 窗口法 窗函数法设计线性相位FIR滤波器步骤： 确定数字滤波器的性能要求：临界频率k，滤波器单位脉冲响应长度N； 根据性能要求，合理选择单位脉冲响应h(n)的奇偶对称性，从而确定理想频率响应Hd(ej)的幅频特性和相频特性； 求理想单位脉冲响应hd(n)，在实际计算中，可对Hd(ej)按M(M远大于N)点等距离采样，并对其求IDFT得hM(n)，用hM(n)代替hd(n)； 选择适当的窗函数w(n)，根据h(n)= hd(n)w(n)求所需设计的FIR滤波器单位脉冲响应； 求H(ej)，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。窗函数的傅式变换W(ej)的主瓣决定了H(ej)过渡带宽。W(ej)的旁瓣大小和多少决定了H(ej)在通带和阻带范围内波动幅度，常用的几种窗函数有：矩形窗w(n)=RN(n) Hann