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

1、数字信号处理实验指导数字信号处理实验指导书信息与控制工程学院电子信息工程系前 言一、实验目的和基本要求数字信号处理实验是和数字信号处理课程同步开设的非独立设课实验，是理论教学的深化和补充。通过实验，使学生巩固和加深对自动控制原理理论知识的理解，进一步培养学生独立分析问题和解决问题的能力，同时注意培养学生综合设计能力、创新能力和实事求是、严谨认真的科学作风以及良好的实验习惯，为今后工作学习打下良好的基础。通过实验学生应达到以下基本要求：（1）通过实验验证数字信号处理的基本理论，并进一步巩固和加深对基本知识的理解。（2）能根据实验指导以及相关资料，综合运用所学知识，深入钻研有关问题，学会自己独立设

2、计实验，分析问题、解决问题，培养一定的实验研究能力和创新能力。（3）能正确使用实验设备，掌握实验原理，熟练运用计算机处理问题。（4）能独立撰写实验报告，准确分析实验结果，及时发现及解决实验中的问题。二、实验项目汇总序号实验项目名称学时实验类型实验要求每组学生1离散信号的MATLAB产生和显示2验证必修12离散系统的差分方程及冲激响应2验证必修13离散系统的频率响应分析2验证必修14用FFT作信号的频谱分析2验证选修15IIR数字滤波器的设计2综合选修16FIR数字滤波器的设计2综合选修1三、实验报告与考核方式要求学生每人独立完成实验，实验结束后按照学院标准格式，自行完成实验报告并上交。按照学校

3、教务处对学生实验考核有关文件精神以及实验过程考勤、操作技能、实验结果和实验报告综合考核。实验成绩占该课程平时成绩的1/3计入总成绩。实验一 常见离散信号的MATLAB产生和图形显示实验目的：1.熟悉MATLAB使用环境及使用方法。2.加深对常用离散时间信号的理解。实验原理：1 单位抽样序列在MATLAB中可以利用zeros()函数实现。注意：MATLAB中数组下标从1开始如果在时间轴上延迟了k个单位，得到即： 2单位阶跃序列在MATLAB中可以利用ones()函数实现。 3正弦序列在MATLAB中 4复正弦序列在MATLAB中5指数序列在MATLAB中实验内容：编制程序产生上述离散时间信号，绘

4、出其图形，并回答问题。参考MATLAB程序清单：% Program P1.1% Generation of a Unit Sample Sequence clf;n = -10:20;% Generate a vector from -10 to 20u = zeros(1,10) 1 zeros(1,20);% Generate the unit sample sequencestem(n,u);% Plot the unit sample sequencexlabel(Time index n);ylabel(Amplitude);title(Unit Sample Sequence);a

5、xis(-10 20 0 1.2);% Program P1.2 (参见教材P_9页的实指数序列MATLAB程序)% Generation of a real exponential sequenceclf;n = 0:35; a = 1.2; K = 0.2;x = K*a.n;stem(n,x);xlabel(Time index n);ylabel(Amplitude);% Program P1.3% Generation of a sinusoidal sequencen = 0:40; f = 0.1; phase = 0; A = 1.5; arg = 2*pi*f*n - pha

6、se; x = A*cos(arg);clf; % Clear old graphstem(n,x); % Plot the generated sequenceaxis(0 40 -2 2);grid; title(Sinusoidal Sequence);xlabel(Time index n);ylabel(Amplitude);axis;实验思考题：1.运行程序P1.1,以产生单位抽样序列。a) 其中命令clf, axis, title, xlabel和ylabel的作用是什么？b) 修改程序，产生带有延时11个样本的延迟单位抽样序列。2.运行程序P1.2，以产生实指数序列。a) 哪个

7、参数控制该序列的增长或衰减率？哪个参数控制该序列的振幅？b) 若参数a小于1，会发生什么情况？将参数a改为0.9，将参数k改为20，再次运行程序P1.2.3.运行程序P1.3，以产生正弦序列。a) 该序列的频率是多少？怎样改变它？哪个参数控制该序列的相位？哪个参数控制该序列的振幅？该序列的周期是多少？b) 用plot命令代替stem命令，重新运行程序，显示的图形与原图形有什么区别？实验二 离散系统的差分方程、冲激响应和卷积分析实验目的：1.学习用差分方程判断系统是否为线性的方法。2.进一步加深对离散系统冲激响应的理解。3.熟悉离散系统的卷积分析方法。实验原理：离散系统其输入、输出关系可用以下差

8、分方程描述：输入信号分解为冲激信号，。系统单位冲激响应，则系统响应为如下的卷积计算式： 当时，hn是有限长度的（n：0，M），称系统为FIR系统；反之，称系统为IIR系统。在MATLAB中，可以用函数y=Filter(p,d,x) 求解差分方程，也可以用函数 y=Conv(x,h)计算卷积。实验内容：一、判断如下系统是否是线性系统假定系统由下列差分方程描述：y(n)-0.4y(n-1)+0.75y(n-2)=2.2403x(n)+2.4908x(n-1)+2.2403x(n-2)试用MATLAB程序仿真该系统，输入三个不同的输入序列x1(n),x2(n)和x(n)=ax1(n)+bx2(n),

9、 计算相应的输出响应y1(n),y2(n)和y(n).从而判断该系统是否是线性系统。% Program P2.1% Generate the input sequencesclf;n = 0:40;a = 2;b = -3;x1 = cos(2*pi*0.1*n);x2 = cos(2*pi*0.4*n);x = a*x1 + b*x2;num = 2.2403 2.4908 2.2403;den = 1 -0.4 0.75;ic = 0 0; % Set zero initial conditionsy1 = filter(num,den,x1,ic); % Compute the outp

10、ut y1ny2 = filter(num,den,x2,ic); % Compute the output y2ny = filter(num,den,x,ic); % Compute the output ynyt = a*y1 + b*y2; d = y - yt; % Compute the difference output dn% Plot the outputs and the difference signalsubplot(3,1,1)stem(n,y);ylabel(Amplitude);title(Output Due to Weighted Input: a cdot

11、x_1n + b cdot x_2n);subplot(3,1,2)stem(n,yt);ylabel(Amplitude);title(Weighted Output: a cdot y_1n + b cdot y_2n);subplot(3,1,3)stem(n,d);xlabel(Time index n);ylabel(Amplitude);title(Difference Signal);二、线性时不变离散时间系统的冲激响应的计算系统：y(n)-0.4y(n-1)+0.75y(n-2)=2.2403x(n)+2.4908x(n-1)+2.2403x(n-2)用MATLAB函数y=im

12、pz(num,den,N)来计算系统的冲激响应。参考程序清单：% Program P2.2% Compute the impulse response yclf;N = 40;num = 2.2403 2.4908 2.2403;den = 1 -0.4 0.75;y = impz(num,den,N);% Plot the impulse responsestem(y);xlabel(Time index n); ylabel(Amplitude);title(Impulse Response); grid;运行程序P2.2，计算上例系统的冲激响应。运行结果：三、用conv命令计算离散卷积%

13、 Program P2.3clf;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 by Convolution); grid;x1 = x zeros(1,8);y1 = filter(h,1,x1);subplot(2,1,2);stem(n,y1);xlabel(

14、Time index n); ylabel(Amplitude);title(Output Generated by Filtering); grid;运行结果：参考:计算任意两个有限长序列的卷积和：% Program 2.4% Illustration of Convolution%a = input(Type in the first sequence = );b = input(Type in the second sequence = );c = conv(a, b);M = length(c)-1;n = 0:1:M;disp(output sequence =);disp(c)st

15、em(n,c)xlabel(Time index n); ylabel(Amplitude);运行结果：Type in the first sequence = -2 0 2 -1 3Type in the second sequence = 1 2 0 -1output sequence = -2 -4 2 5 1 4 1 -3实验思考题：1运行程序P2.1，对由加权输入得到的y(n)与在相同系数下输出y1(n) 和y2(n)相加得到的yt(n)进行比较,这两个序列是否相等?该系统是线性系统吗?用几组不同的权系数a和b的值以及几组不同的输入频率,重新运行程序进一步验证.2.运行程序2.2,若

16、想画出系统冲激响应的前45个样本,怎样修改程序?实验三 离散系统的频率响应分析和零、极点分布实验目的：1.加深对离散系统的频率响应分析理解。2.加深对离散系统的零、极点分布的概念理解。实验原理：离散系统的时域方程为其变换域分析方法如下：频域 系统的频率响应为 Z域: 系统的转移函数为 分解因式，其中和称为零、极点。在MATLAB中，可以用函数z，p，K=tf2zp（num，den）求得有理分式形式的系统函数的零、极点；用函数zplane（z,p）绘出零、极点分布图；用函数zplane（num,den）直接绘出有理分式形式的系统函数的零、极点分布图。用函数 r，p，k=residuez（num，

17、den）完成部分分式展开计算；用函数sos=zp2sos（z，p，K）完成将高阶系统分解为2阶系统（即“二阶节” ）的级联。实验内容：求系统 的零、极点，并用因式分解形式表示H(z),并确定收敛域。实验要求：编程实现系统参数输入，绘出零、极点分布图和幅度频率响应曲线。程序清单： % Program 3.1num = input(Type in the numerator coefficients = );den = input(Type in the denominator coefficients = );z,p,k = tf2zp(num,den);m = abs(p);disp(Zero

18、s are at);disp(z);disp(Poles are at);disp(p);disp(Gain constant);disp(k);disp(Radius of poles);disp(m);sos = zp2sos(z,p,k);disp(Second-order sections);disp(real(sos);subplot(2,1,1);zplane(num,den);w = 0:pi/255:pi;h = freqz(num,den,w);subplot(2,1,2);plot(w/pi,abs(h);gridxlabel(w/pi);ylabel(Amplitude)

19、;运行结果：Type in the numerator coefficients = Type in the denominator coefficients = Zeros are atPoles are atGain constantRadius of polesSecond-order sections其零极点图为：实验思考题：运行程序P3.1, 解释程序的运行结果，并利用运行结果，将该系统的系统函数写成几个“二阶节”级联的形式，并画出其级联形式的结构流图。实验四 用FFT作信号的频谱分析实验目的：1.加深对离散信号的DTFT和DFT及其相互关系的理解。2.对快速傅里叶变换FFT算法的理

20、解。3.理解用FFT对典型信号的频谱分析方法。实验原理：序列xn 的DTFT定义： N点序列xn 的DFT定义：在MATLAB中，对形式为的DTDFT可以用函数H=Freqz（num，den，w）计算；可以用函数U=fft（u，N）和u=ifft（U，N）计算N点序列的DFT正、反变换。实验内容：1、计算N点序列u(n)的M点DFT。程序要求输入数据N，M。M必须大于N。这里我们取N=8，M=16。然后再取M=32再运行一遍程序。比较两次的结果。程序清单：% Program 3.1% Illustration of DFT Computation% Read in the length N o

21、f sequence and the desired% length M of the DFTN = input(Type in the length of the sequence = );M = input(Type in the length of the DFT = );% Generate the length-N time-domain sequenceu = ones(1,N);% Compute its M-point DFTU = fft(u,M);% Plot the time-domain sequence and its DFTt = 0:1:N-1;stem(t,u)

22、title(Original time-domain sequence)xlabel(Time index n); ylabel(Amplitude)pausesubplot(2,1,1)k = 0:1:M-1;stem(k,abs(U)title(Magnitude of the DFT samples)xlabel(Frequency index k); ylabel(Magnitude)subplot(2,1,2)stem(k,angle(U)title(Phase of the DFT samples)xlabel(Frequency index k); ylabel(Phase);运

23、行结果：原时间序列(8点)：对其进行16点DFT：对其进行32点DFT：2、求正弦序列的DFT：%Program 3.2close all;t=linspace(1e-3,100e-3,10);xn=sin(100*2*pi*t);N=length(xn);WNnk=dftmtx(N);Xk=xn*WNnk;figure;subplot(1,2,1),stem(1:N,xn),title(时域离散序列x(n);subplot(1,2,2),stem(1:N,abs(Xk),title(x(n)的DFT变换结果);运行结果：3求实指数序列x(n)=an的DFT:%Program 3.3clear

24、N=8;a=0.7;n=0:7;xn=a.n;Xk=fft(xn,N);subplot(3,1,1)stem(n,xn,.k); xlabel(n); ylabel(x(n);title(实指数序列);axis(0,8,0,1.5)subplot(3,1,2)stem(n,abs(Xk),.k);xlabel(k); ylabel(absX(k);title(实指数序列8点DFT幅度);axis(0,8,0,5)subplot(3,1,3)stem(n,angle(Xk),.k);xlabel(k); ylabel(angleX(k);title(实指数序列8点DFT相位);axis(0,8,

25、-1.5,1.5)实验思考题：1. 简要回答离散信号的DTFT与DFT之间的关系.2. 运行程序P3.3, 通过实验结果来分析实数序列的DFT的圆周共轭对称性质.实验五 IIR数字滤波器设计实验目的：1. 加深对IIR型数字滤波器的理解。2. 掌握IIR数字滤波器的结构和设计方法。实验原理：MATLAB中滤波器的分析和实现可用以下函数完成：1、freqs函数：模拟滤波器的频率响应；2、freqz函数：数字滤波器的频率响应；3、buttord函数：ButterWorth滤波器阶数的选择；4、butter函数：ButterWorth滤波器设计；5、cheb1ord函数：Chebyshev型滤波器阶

26、数计算；6、cheby1函数：Chebyshev型滤波器设计；7、impinvar函数：脉冲响应不变法设计数字滤波器；8、bilinear函数：双线性变换法设计数字滤波器。实验内容：1、模拟滤波器的频率响应(freqs函数)系统传递函数为的模拟滤波器，在MATLAB中可以用以下程序来实现：%Program 4.1b=0.2 0.3 1; a=1 0.4 1; w=logspace(-1,1); %产生从10-1到101之间等间距点，即50个频率点freqs(b,a,w) %根据输入的参数绘制幅度谱和相位谱2、数字滤波器的频率响应(freqz函数):系统传递函数为的数字滤波器在MATLAB中可以

27、用以下程序来实现：%Program 4.2a=1 0.4 1; b=0.2 0.3 1; freqz(b,a,128) %根据输入的参数绘制幅度谱和相位谱，得到0到之间128个点处的频率响应3、ButterWorth模拟和数字滤波器（1）buttord函数：ButterWorth滤波器阶数的选择。调用格式：n,Wc=buttord(Wp,Ws,Rp,Rs)，在给定滤波器性能的情况下（通带临界频率Wp、阻带临界频率Ws、通带内最大衰减Rp和阻带内最小衰减Rs），计算ButterWorth滤波器的阶数n和截止频率Wn。相同参数条件下的模拟滤波器则调用格式为：n,Wc=buttord(Wp,Ws,R

28、p,Rs,s) （2）butter函数：ButterWorth滤波器设计。调用格式：b,a=butter(n,Wc)，根据阶数n和截止频率Wc计算ButterWorth滤波器分子分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）。相同参数条件下的模拟滤波器则调用格式为：b,a=butter(n,Wc,s) 采样频率为1Hz，通带临界频率fp =0.2Hz，通带内衰减小于1dB（Rp=1）；阻带临界频率fs=0.3Hz，阻带内衰减大于25dB（As=25）。设计一个数字滤波器满足以上参数。%Program 4.3n,Wc=buttord(0.2,0.3,1,25); b,a=butter

29、(n,Wc); freqz(b,a,512,1); 4、Chebyshev模拟和数字滤波器 （1）cheb1ord函数：Chebyshev型滤波器阶数计算。调用格式：n,Wc=cheb1ord(Wp,Ws,Rp,Rs)，在给定滤波器性能的情况下（通带临界频率Wp、阻带临界频率Ws、通带内波纹Rp和阻带内衰减Rs），选择Chebyshev型滤波器的最小阶n和截止频率Wc。（2）cheby1函数：Chebyshev型滤波器设计。调用格式：b,a=cheby1(n,Rp,Wc)，根据阶数n、通带内波纹Rp和截止频率Wc计算Chebyshev型滤波器分子分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）。用Chebyshev型实现上例中的滤波器%Program 4.4n,Wn=cheb1ord(0.2,0.3,1,25); b,a=cheby1(n,1,Wc