数字信号处理实验_2014

1、实验1 常见离散信号产生和实现一、实验目的1、加深对常用离散信号的理解；2、熟悉使用MATLAB在时域中产生一些基本的离散时间信号。 二、实验原理MATLAB 语言提供了一系列函数用来产生信号，如exp, sin, cos, square, sawtooth，ones, zeros等函数。1. 基本信号序列1) 单位抽样序列 在MATLAB中可以利用zeros()函数实现。x=1 zeros(1, n-1)示范程序:% Program P1_1% Generation of a Unit Sample Sequence clf;% Generate a vector from -10 to 2

2、0n = -10:20;% Generate the unit sample sequenceu = zeros(1,10) 1 zeros(1,20);% Plot the unit sample sequencestem(n,u);xlabel('Time index n');ylabel('Amplitude');title('Unit Sample Sequence');axis(-10 20 0 1.2); 如果在时间轴上延迟了k个单位，得到即： Request 1: 编写一个的函数。2) 单位阶跃序列 在MATLAB中可以利用ones

3、()函数实现。 Request 2: 编写一个的函数。3) 实指数序列MATLAB实现：4) 复指数序列MATLAB实现：5) 随机序列 MATLAB提供了两种随机信号： rand(1,N) 产生0,1上均匀分布的随机矢量。 randn(1,N) 产生均值为0，方差为1的高斯随机序列，即白噪声序列。2. 基本周期波形1) 方波MATLAB工具箱函数square可以产生方波;t=0:0.1*pi:6*pi;y=square(t);plot(t,y);xlabel(时间 t);ylabel(幅度y); axis(0 7*pi -1.5 1.5);2) 正弦波在MATLAB中3) 锯齿波工具箱函数s

4、awtooth函数可以产生锯齿波Fs=10000;t=0:1/Fs:1.5; %抽样长度1.5s, 抽样频率为10kHzx=sawtooth(2*pi*50*t); %信号频率为50Hzplot(t, x);axis()0 0.2 -1 1; %画出0.2秒的波形3. 基本非周期波形工具箱函数chirp能产生一种扫射频率信号，其特点是信号的瞬时频率随时间按照一定规律变化 t=0:1/1000:2 %抽样频率1kHz, 抽样时间2s. x=chirp(t, 0.1, 150) %0时刻为DC信号，1s时频率为150Hz。 specgram(x, 256, 1000, 256, 250);4. s

5、inc信号MATLAB实现： t=linspace(-5,5); x=sinc(t); plot(t,x);5. 序列的操作1) 信号加 x(n)=x1(n)+x2(n)MATLAB实现： x=x1+x2;注意：x1和x2序列应该具有相同的长度，位置对应，才能相加。2) 信号乘 x(n)=x1(n)*x2(n)MATLAB实现： x=x1.*x2; %数组乘法3) 改变比例 y(n)=k*x(n)MATLAB实现： y=k*x;4) 折叠 y(n)=x(-n);MATLAB实现： y=fliplr(x);5) 抽样和 MATLAB实现： y=sum(x(n1:n2);6) 抽样积 MATLAB

6、实现： y=prod(x(n1:n2);7) 信号能量 MATLAB实现： Ex=sum(abs(x).2);8) 信号功率 MATLAB实现：Px=sum(abs(x).2)/ N;三、预习要求1、预先阅读附录A（MATLAB基础介绍）；2、熟悉实验内容，预先编写好程序，准备好预习报告。四、实验内容 编制程序产生上述5种信号，长度可输入确定，函数需要的参数可输入确定，并绘出其图形。五、实验报告要求（1）简述实验目的及实验原理。（2）解释示范程序中MATLAB各个语句的功能（3）按实验步骤附上实验过程中的实验程序和信号序列曲线，并对所得结果进行分析和解释。（4）Request3: 讨论复指数序

7、列的性质。绘出信号，当、时、时的信号实部和虚部图；当时呢？此时信号周期为多少？实验2 离散系统的时域分析一、实验目的1、熟悉并掌握离散系统的差分方程表示法；2、加深对冲激响应和卷积分析方法的理解。 二、实验原理在时域中，离散时间系统对输入信号或者延迟信号进行运算处理，生成具有所需特性的输出信号，具体框图如下:其输入、输出关系可用以下差分方程描述：输入信号分解为冲激信号,记系统单位冲激响应，则系统响应为如下的卷积计算式：当时，hn是有限长度的（），称系统为FIR系统；反之，称系统为IIR系统。三、预习要求1、在MATLAB中，熟悉利用函数实现差分方程的仿真；2、在MATLAB中，熟悉用函数 计算

8、卷积，用求系统冲激响应的过程。 四、实验内容1、用MATLAB求系统响应1） 卷积的实现线性移不变系统可由它的单位脉冲响应来表征。若已知了单位脉冲响应和系统激励就可通过卷积运算来求取系统响应，即程序：%Illustration of Convolutionx=input(Type in the input sequence=);h=input(Type in the impulse response sequence=);y=conv(x,h);N=length(y-1);n=0:1:N;disp(output sequence=); disp(y);stem(n,y);xlabel(Time

9、 index n); ylable(Amplitude);输入为： x=-2 0 1 -1 3 h=1 2 0 -12） 单位脉冲响应的求取线性时不变因果系统可用MATLAB的函数filter来仿真 y=filter(b,a,x); 其中，x和y是长度相等的两个矢量。矢量x表示激励，矢量a，b表示系统函数形式滤波器的分子和分母系数，得到的响应为矢量y。例如计算以下系统的单位脉冲响应y(n)+0.7y(n-1)-0.45y(y-2)-0.6y(y-3)=0.8x(n)-0.44x(n-1)+0.36x(n-2)+0.02x(n-3) 程序实现： % Illustration of Impulse

10、 Response Computation N=input(Desired impuse response length=);b=input(Type in the vector b=);a=input(Type in the vector a=);x=1 zeros(1,N-1);y=filter(b,a,x);k=0:1:N-1;stem(k,y);xlabel(Time index n); ylable(Amplitude);输入数据：N=41b=0.8 -0.44 0.36 0.02a=1 0.7 -0.45 -0.62、以下程序中分别使用conv和filter函数计算h和x的卷积y和

11、y1，运行程序，并分析y和y1是否有差别，为什么要使用xn补零后的x1来产生y1；具体分析当hn有i个值，xn有j个值，使用filter完成卷积功能，需要如何补零？clf;h = 3 2 1 -2 1 0 -4 0 3; %impulse responsex = 1 -2 3 -4 3 2 1; %input sequencey = conv(h,x);n = 0:14;subplot(2,1,1);stem(n,y);xlabel('Time index n'); ylabel('Amplitude');title('Output Obtained b

12、y Convolution'); grid;x1 = x zeros(1,8);y1 = filter(h,1,x1);subplot(2,1,2);stem(n,y1);xlabel('Time index n'); ylabel('Amplitude');title('Output Generated by Filtering'); grid; 3、编制程序求解下列两个系统的单位冲激响应，分别用filter 和 impz实现，并绘出其图形。给出理论计算结果和程序计算结果并讨论。五、实验报告要求（1）简述实验目的及实验原理。（2）按实验

13、步骤附上实验过程及结果。（3）总结实验中的主要结论。实验3 离散系统的变换域分析一、实验目的1、熟悉对离散系统的频率响应分析方法；2、加深对零、极点分布的概念理解。 二、实验原理离散系统的时域方程为其变换域分析方法如下：频域:系统的频率响应为:Z域:系统函数为:分解因式:，其中和称为零、极点。三、预习要求1. 在MATLAB中，熟悉函数tf2zp、zplane、freqz、residuez、zp2sos的使用，其中：z，p，K=tf2zp（num，den）求得有理分式形式的系统函数的零、极点；zplane（z，p）绘制零、极点分布图；h=freqz(num,den,w)求系统的单位频率响应；r

14、，p，k=residuez（num，den）完成部分分式展开计算；sos=zp2sos（z，p，K）完成将高阶系统分解为2阶系统的串联。 2. 阅读扩展练习中的实例，学习频率分析法在MATLAB中的实现； 3. 编程实现系统参数输入，绘出幅度频率响应和相位响应曲线和零、极点分布图。 四、实验内容求系统的零、极点和幅度频率响应和相位响应。五、扩展练习例1: 求下列直接型系统函数的零、极点，并将它转换成二阶节形式解:用MATLAB计算程序如下：num=1 -0.1 -0.3 -0.3 -0.2;den=1 0.1 0.2 0.2 0.5;z,p,k=tf2zp(num,den);m=abs(p);

15、disp('零点');disp(z);disp('极点');disp(p);disp('增益系数');disp(k);sos=zp2sos(z,p,k);disp('二阶节');disp(real(sos);zplane(num,den)输入到“num”和“den”的分别为分子和分母多项式的系数。计算求得零、极点增益系数和二阶节的系数：零点:0.9615-0.5730-0.1443 + 0.5850i-0.1443 - 0.5850i极点:0.5276+0.6997i0.5276-0.6997i-0.5776+0.5635i-0.

16、5776-0.5635i增益系数:1二阶节:1.0000 -0.3885 -0.5509 1.0000 1.15520 0.65111.0000 0.28850 0.36300 1.0000 -1.0552 0.7679系统函数的二阶节形式为：极点图见图3.1:图3.1 系统函数的零、极点图例2: 求差分方程所对应的系统的频率响应。解:差分方程所对应的系统函数为：用MATLAB计算的程序如下：k=256;%采样频率为256Hznum=0.8 -0.44 0.36 0.02;den=1 0.7 -0.45 -0.6;w=0:pi/k:pi;h=freqz(num,den,w);subplot(2

17、,2,1);plot(w/pi,real(h);grid;title('实部');xlabel('omega/pi');ylabel('幅度');subplot(2,2,2);plot(w/pi,imag(h);grid;title('虚部');xlabel('omega/pi');ylabel('Amplitude');subplot(2,2,3);plot(w/pi,abs(h);grid;title('幅度谱');xlabel('omega/pi');ylabe

18、l('幅值');subplot(2,2,4);plot(w/pi,angle(h);grid;title('相位谱');xlabel('omega/pi');ylabel('弧度'); 图3.2六、实验报告要求（1）简述实验目的及实验原理。（2）按实验步骤附上实验过程（3）结合扩展练习，总结实验中的主要结论。实验4 FFT算法的应用一、实验目的1、加深对离散信号的DFT的理解；2、在MATLAB中实现FFT算法。二、实验原理N点序列的DFT和IDFT变换定义式如下：,利用旋转因子具有周期性，可以得到快速算法（FFT）。在MATLA

19、B中，可以用函数和计算N点序列的DFT正、反变换。三、预习要求1、在MATLAB中，熟悉函数fft、ifft的使用；2、阅读扩展练习中的实例，学习在MATLAB中的实现FFT算法的实现；3、利用MATLAB编程完成计算，绘出相应图形。并与理论计算相比较，说明实验结果的原因。 四、实验内容1、已知某序列在单位圆上的N=64等分样点的Z变换为:。用N点IFFT程序计算出。五、扩展练习例1：对连续的单一频率周期信号 按采样频率采样，截取长度N分别选N =20和N =16，观察其DFT结果的幅度谱。解：此时离散序列，即k=8。用MATLAB计算并作图，函数fft用于计算离散傅里叶变换DFT，程序如下：

20、k=8;n1=0:1:19;xa1=sin(2*pi*n1/k);subplot(2,2,1)plot(n1,xa1)xlabel('t/T');ylabel('x(n)');xk1=fft(xa1);xk1=abs(xk1);subplot(2,2,2)stem(n1,xk1)xlabel('k');ylabel('X(k)');n2=0:1:15;xa2=sin(2*pi*n2/k);subplot(2,2,3)plot(n2,xa2)xlabel('t/T');ylabel('x(n)');x

21、k2=fft(xa2);xk2=abs(xk2);subplot(2,2,4)stem(n2,xk2)xlabel('k');ylabel('X(k)');图 4.1 不同的截取长度的正弦信号及其DFT结果计算结果示于图1，(a)和(b)分别是N=20时的截取信号和DFT结果，由于截取了两个半周期，频谱出现泄漏；(c) 和(d) 分别是N=16时的截取信号和DFT结果，由于截取了两个整周期，得到单一谱线的频谱。上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。六、实验报告要求（1）简述实验目的及实验原理。（2）结合实验中所得给定典型序列幅频特性曲

22、线，与理论结果比较，并分析说明误差产生原因以及用FFT作谱分析时有关参数的选择方法。按实验步骤附上实验过程（3）结合扩展练习，总结实验中的主要结论。（4）思考：如果周期信号的周期预先不知道，如何用FFT进行谱分析？实验5 无限冲激响应数字滤波器设计一、实验目的1、掌握双线性变换法及脉冲相应不变法设计IIR数字滤波器的具体设计方法;2、熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。二、实验原理在MATLAB中，可以用下列函数辅助设计IIR数字滤波器：1)利用buttord和cheb1ord可以确定低通原型巴特沃斯和切比雪夫滤波器的阶数和截止频率；如：求阶数N

23、, Wn = cheb1ord (Wp, Ws, Rp, Rs,s)选择项说明：high类别。缺省为low； s模/数，缺省为数Rp即 ap, Rs即as ； WnChebyshev自然频率（3dB频率），数字设计： Wpwp/p Wsws/p 。2)num,den=butter（N,Wn）(巴特沃斯)和num,den=cheby1（N,Wn）,num,den=cheby2（N,Wn）(切比雪夫1型和2型)可以进行滤波器的设计；3）lp2hp，lp2bp，lp2bs可以完成低通滤波器到高通、带通、带阻滤波器的转换；4)使用bilinear可以对模拟滤波器进行双线性变换，求得数字滤波器的传输函数

24、系数；5)利用impinvar可以完成脉冲响应不变法的模拟滤波器到数字滤波器的转换。三、预习要求1. 在MATLAB中，熟悉函数butter、cheby1、cheby2的使用，其中：num,den=butter（N,Wn）巴特沃斯滤波器设计；num,den=cheby1（N,Wn）切比雪夫1型滤波器设计；num,den=cheby2（N,Wn）切比雪夫2型滤波器设计。2. 阅读扩展练习中的实例，学习在MATLAB中进行数字滤波器的设计； 3. 给出IIR数字滤波器参数和滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。 4. 掌握双线性变换法和直接法。双线性变换法：

25、a求模拟边缘频率Wp = 2/T*tan(wp/2); Ws = 2/T*tan(ws/2)b. 求巴特沃斯滤波器阶数 N, Wn = buttord(Wp, Ws, 1, 15, s)c求巴特沃斯滤波器系数 B, A = butter(N, Wn, s)d双线性变换求相应的数字滤波器 num, den = bilinear(B,A,0.5)直接法:a. 求数字滤波器阶数 N, Wn = buttord (wp/(fs/2), ws/(fs/2), 1, 15)b. 求数字滤波器系数 num, den = butter(N, Wn)四、实验内容利用MATLAB编程，用脉冲响应不变法和双线性变换

26、法设计一个数字带通滤波器，指标要求如下：通带边缘频率：，通带峰值起伏：；阻带边缘频率：，最小阻带衰减： 。五、扩展练习例1：设采样周期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(nu

27、m2,den2);f=w/pi*2000; %采样定理plot(f,abs(h1),'-.',f,abs(h2),'-');grid;xlabel('频率/Hz ')ylabel('幅值/dB')程序中第一个butter的边界频率2×1000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特

28、性变差。图5.1给出了MATLAB计算的结果。图 6.1 例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

29、=freqz(num,den);f=w/pi*500;plot(f,20*log10(abs(h);axis(0,500,-80,10);grid;xlabel('')ylabel('幅度/dB') 图5.2给出了MATLAB计算的结果。图 5.2例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);

30、w3=2*400*tan(2*pi*120/(2*400); N,wn=buttord(w1 w2,0.1 w3,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') 直接法：w1=90/(400/2); w2=110/(400/2);

31、w3=120/(400/2); N,wn=buttord(w1 w2,0.1 w3,3,10);B,A=butter(N,wn); h,w=freqz(B,A);f=w/pi*200; plot(f,20*log10(abs(h); axis(40,160,-30,10); grid; xlabel('频率/kHz') ylabel('幅度/dB') 图5.3给出了MATLAB计算的结果。图 5.3 例4: 一数字滤波器采样频率fs=1kHz，要求滤除100Hz的干扰，其dB的边界频率为95Hz和105Hz，原型归一化低通滤波器为：w1=95/(1000/2);

32、 w2=105/(1000/2); 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')图5.4为MATLAB的计算结果。图 5.4 六、实验报告要求（1）简述实验目的及实验原理。（2）按实验步骤附上实验过程及结果。（3）结合扩展练习，总结实验中的主要结论。 实验6 有限冲激响应数字滤波器设计一、实验目的： 1、加深对数字滤波器的常用指标理解。 2、学习数字滤波器的设计方法。 二、实验原理：低通滤波器的常用指标：（1）通带边缘频率；（2）阻带边缘频率；（3）通带起伏；（4）通带峰值起伏，（5）阻带