数字信号处理实验指导书

1、实验一 离散系统的时域分析一. 实验目的（1）掌握求系统响应的方法。（2）掌握时域离散系统的时域特性。二.实验原理与方法在时域中，描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应，本实验仅在时域求解。在计算机上适合用递推法求差分方程的解，最简单的方法是采用MATLAB语言的工具箱函数filter函数。也可以用MATLAB语言的工具箱函数conv函数计算输入信号和系统的单位脉冲响应的线性卷积，求出系统的响应。三实验内容及结果（1）编制程序，包括产生输入信号、单位脉冲响应序列的子程序，用fi

2、lter函数或conv函数求解系统输出响应的主程序。程序中要有绘制信号波形的功能。给定一个低通滤波器的差分方程为： 输入信号： a) 分别求出系统对两个输入信号的响应序列，并画出其波形。 b) 求出系统的单位冲响应，画出其波形。% A=1,-0.9;% B=0.05,0.05;% x1n=1 1 1 1 1 1 1 1 zeros(1,25);% x2n=ones(1,32);% y1n=impz(B,A,33);% n=0:length(y1n)-1;% figure(1);% stem(n,y1n);% title('1102303005 y1(n)');% box on;

3、% xlabel('n');% ylabel('y1(n)');% x2n=ones(1,33);% y2n=filter(B,A,x2n);% n=0:length(y2n)-1;% figure(2);% stem(n,y2n);% title('1102303005 y2(n)');% box on;% xlabel('n');% ylabel('y2(n)');% hn=filter(B,A,x1n);% n=0:length(hn)-1;% figure(3);% stem(n,hn);% title(&

4、#39;1102303005 h(n)');% box on% xlabel('n');%ylabel('h(n)');（2）给定系统的单位脉冲响应为：用线性卷积法分别求两个系统对输入信号的输出响应，并画出波形。% x1n=1 1 1 1 1 1 1 1 zeros(1,24);% h1n=1 1 1 1 1 1 1 1 1 1 zeros(1,22);% h2n=1 2.5 2.5 1 zeros(1,28);% y21n=conv(h1n,x1n);% y22n=conv(h2n,x1n);% n=0:length(y21n)-1;% figure(

5、4);% stem(n,y22n);% title('1102303005 x1(n)');% box on;% xlabel('n');% ylabel('x1(n)');% n=0:length(y21n)-1;% figure(5);% stem(n,y21n);% title('1102303005 y21(n)');% box on;% xlabel('n');% ylabel('y21(n)');% n=0:length(y22n)-1;% figure(6);% stem(n,y22n

6、);% title('1102303005 y22(n)');% box on;% xlabel('n');% ylabel('y22(n)');（3）测量人耳辨别回声的最小时间设计一个混合声音数字系统，利用MATLAB语言实现声音的混合。A=1,-0.9;B=0.05,0.05;x1n=1 1 1 1 1 1 1 1 zeros(1,25);x2n=ones(1,32);xn=1 zeros(1,32);y1n=filter(B,A,33);h=filter(B,A,xn);load mtlb.matyn=conv(h,mtlb)n=0:len

7、gth(yn)-1;figure(1);stem(n,yn);title('1102303005 yn)');box on;xlabel('n');ylabel('yn');四思考题（1）简述在时域求系统响应的方法。有两种，一种是通过解差分方程求得系统的输出，注意要合理的选择初始条件，一种是已知系统的单位脉冲响应，通过球输出信号和系统的单位脉冲的线性卷积求得性输出。（2）简述判断系统稳定性的方法。如果信号经过低通滤波器，把信号的高频部分滤去，时域信号的的剧烈变化将被平滑。实验二 离散系统的频域分析一、实验目的 （1）熟悉对离散系统的频率响应分析方

8、法。 （2）加深对零、极点分布的概念理解。二、实验原理 离散系统的时域方程为： 其变换域分析方法如下： 频域: 系统的频率响应为: Z域: 系统函数为: 分解因式: ， 其中和称为零、极点。三、实验内容及结果 （1）求差分方程 所对应的系统的频率响应。程序：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,2,1);plot(w/pi,real(h);grid;title(1102303005 实部)；xlabel('omega/p

9、i');ylabel('幅度');subplot(2,2,2);plot(w/pi,imag(h);grid;title('1102303005 虚部)；xlabel('omega/pi');ylabel('Amplitude');subplot(2,2,3);plot(w/pi,abs(h);grid;title(1102303005 幅度谱)；xlabel('omega/pi');ylabel('幅值');subplot(2,2,4);plot(w/pi,angle(h);grid;title(

10、1102303005 相位谱)；xlabel('omega/pi');ylabel('弧度'); 结果：（2）求系统的零、极点和幅度频率响应和相位响应。程序：clc;clear ;num=0.0528 0.0797 0.1295 0.1295 0.0797 0.0528;den=1 -1.8107 2.4947 -1.8801 0.9537 -0.2336;figure(1)zplane(num,den);grid;z,p,k=tf2zp(num,den);disp('零点=');disp(z);disp('极点=');disp(

11、p);disp('增益系数=');disp(k);k=256; %采样频率为K=256Hzw=0:pi/k:pi;h=freqz(num,den,w);figure(2)subplot(2,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('Ampl

12、itude');subplot(2,2,3);plot(w/pi,abs(h);grid;title('幅度谱');xlabel('omega/pi');ylabel('幅值');subplot(2,2,4);plot(w/pi,angle(h);grid;title('相位谱');xlabel('omega/pi');ylabel('弧度'); 结果：的零、极点和幅度频率响应和相位响应。实验三 FFT算法的应用一、实验目的 （1）加深对离散信号的DFT的理解。 （2）熟悉FFT算法的应用。

13、二、实验原理 N点序列的DFT和IDFT变换定义式如下： , 利用旋转因子具有周期性，可以得到快速算法（FFT）。 在MATLAB中，可以用函数和 计算N点序列的DFT正、反变换。三、实验内容及结果 （1）观察一句话的频谱特点：将自己录制的一段话均匀分段，并用DFT计算每段信号的频谱。程序：Fs=8000;T=4;x=wavrecord(T*Fs,Fs);L=length(x);t=(1:L)/Fs;subplot(411);plot(t,x);xlabel('t/s');ylabel('x(t)');title('1102303005)；axis(0

14、4 -0.7 0.7);subplot(413);specgram(x,256,Fs,256,250);xlabel('t/s');ylabel('f/Hz');title('1102303005)；结果： （2）2N点实数序列 N=64。一次算出，并绘出 的图形。程序：N=64;n=0:2*N-1;x1=cos(2*pi*7*n/N);x2=0.5*cos(2*pi*19*n/N);x=x1+x2;subplot(2,1,1);stem(n,x);xlabel('t/T');ylabel('x(n)');title(&#

15、39;1102303005)；xk=fft(x,N);disp(xk);xk=abs(xk);k=0:N-1;subplot(2,1,2);stem(k,xk);xlabel('k');ylabel('X(k)');title('1102303005)；结果： （3）已知某序列在单位圆上的N=64等分样点的Z变换为:。用N点IFFT程序计算出。程序：N=64;k=0:N-1;xk=1./(1-0.8*exp(2*pi*k/N)*-1*j);subplot(2,1,1)stem(k,xk)xlabel('k');ylabel('X(

16、k)');title('1102303005)；axis(0 63 0 6)x=ifft(xk,N);subplot(2,1,2)stem(k,x)xlabel('t/T');ylabel('x(n)');axis(0 63 0 1.5)title('1102303005)；结果：四、思考题（1）如果周期信号的周期预先不知道，如何用FFT进行谱分析？ 答：周期信号的周期预先不知道时,可先截取M点进行DFT,再将截取长度扩大1倍截取,比较结果,如果二者的差别满足分析误差要求,则可以近似表示该信号的频谱,如果不满足误差要求就继续将截取长度加倍

17、,重复比较,直到结果满足要求。（2）如何选择FFT的变换区间？（包括非周期信号和周期信号） 答：一、对于非周期信号：有频谱分辨率F，而频谱分辨率直接和FFT的变换区间有关，因为FFT能够实现的频率分辨率是2/N.因此有最小的N>2/F。就可以根据此式选择FFT的变换区间。 二、对于周期信号，周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。实验四 IIR数字滤波器的设计1、实验目的（1）掌握双线性变换法及脉冲相应不变法设计IIR数字滤波器的具体设计方法;（2）熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。2、

18、实验原理在MATLAB中，可以用下列函数辅助设计IIR数字滤波器：（1）利用buttord和cheb1ord确定低通原型巴特沃斯和切比雪夫滤波器的阶数和截止频率；如：求阶数N, 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）lp2h

19、p，lp2bp，lp2bs可以完成低通滤波器到高通、带通、带阻滤波器的转换；（4）使用bilinear可以对模拟滤波器进行双线性变换，求得数字滤波器的传输函数系数；（5）利用impinvar可以完成脉冲响应不变法的模拟滤波器到数字滤波器的转换。3、实验内容利用MATLAB编程，用脉冲响应不变法和双线性变换法设计一个数字低通滤波器，指标要求如下：通带边缘频率：，通带峰值起伏：； 阻带边缘频率：，最小阻带衰减：。画出两种设计方法下的损耗函数和幅频特性。%脉冲响应不变法Fs=400;T=1/Fs;wp=0.3*pi;ws=0.8*pi;Wp=wp*Fs;Ws=ws*Fs;Ap=1;As=40;N,W

20、c=buttord(Wp,Ws,Ap,As,'s');B,A=butter(N,Wc,'s');D,C=impinvar(B,A,Fs);Hz=freqz(D,C,W);f=(W*Fs)/(2*pi);figure(1)plot(f,abs(Hz);grid on;title('1102303005脉冲响应不变法、幅频特性');xlabel('频率/Hz');ylabel('归一化幅值');figure(2)a=20*log10(abs(Hz);plot(f,a);grid on;title('110230

21、3005脉冲响应不变法、衰减函数')xlabel('频率/Hz');ylabel('幅度/dB');%双线性变化Fs=400;T=1/Fs;wp=0.3*pi;ws=0.8*pi;Wp=(2/T)*tan(wp/2);Ws=(2/T)*tan(ws/2);Ap=1;As=40;N,Wc=buttord(Wp,Ws,Ap,As,'s');B,A=butter(N,Wc,'s');D,C=bilinear(B,A,Fs);Hz,W=freqz(D,C);f=(W*Fs)/(2*pi);figure(1)plot(f,abs(H

22、z);grid on;title('1102303005双线性变换法、幅频特性');xlabel('频率/Hz');ylabel('归一化幅值');figure(2)a=20*log10(abs(Hz);plot(W/pi,a);grid on;title('1102303005双线性变换法、幅频特性');xlabel('频率/Hz');ylabel('幅度/dB');实验五 FIR数字滤波器的设计1、实验目的： （1）加深对数字滤波器的常用指标理解。 （2）学习数字滤波器的设计方法。 2、实验原理

23、：（1）低通滤波器的常用指标：（a）通带边缘频率；（b）阻带边缘频率；（c）通带起伏；（d）通带峰值起伏，（e）阻带起伏，最小阻带衰减。（2）在MATLAB中，熟悉函数fir1、kaiserord 、remezord、remez的使用； B = fir1(n,Wn,'high','noscale')设计滤波器；n,Wn,beta,ftype = kaiserord(f,a,dev)估计滤波器阶数；（默认情况下，过滤器的第一通带中心为精确的一个窗口。使用一个尾随'noscale'参数来防止这种缩放remezord、remez一起使用设计等波纹线性相

24、位FIR滤波器）3、实验内容：利用MATLAB编程,用窗函数法设计FIR数字滤波器，指标要求如下：通带边缘频率：，通带峰值起伏：。阻带边缘频率：，最小阻带衰减：。clear allclcwp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;width1 =wp1-ws1;width2 =ws2-wp2;width=min(width1,width2)N=ceil(6.62*pi/width)n =0: (N-1);a =(N-1)/2+0.00001;wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;hd=(sin(wc2*(n-a)-si

25、n(wc1*(n-a)./(pi*(n-a);win =(hanning(N)'h=hd.*win;freqz(h);subplot(2,1,1)title('1102303005 幅频特性')xlabel('频率')ylabel('db')subplot(2,1,2)title('1102303005 相频特性')xlabel('频率')ylabel('db')实验六 数字信号处理在双音多频拨号系统中的应用1、 DTMF的组成在电话中，数字09的中每一个都用两个不同的单音频传输，所用的8个

26、频率分成高频带和低频带两组，低频带有四个频率：679Hz,770Hz,852Hz和941Hz；高频带也有四个频率：1209Hz,1336Hz,1477Hz和1633Hz.。每一个数字均由高、低频带中各一个频率构成，例如1用697Hz和1209Hz两个频率，信号用表示，其中。这样8个频率形成16种不同的双频信号。2、 电话中的双音多频（DTMF）信号的产生与检测（1）双音多频信号的产生假设时间连续的DTMF信号用表示，式中f1和f2是按照上表选择的两个频率，f1代表低频带中的一个频率，f2代表高频带中的一个频率。采用数字方法产生DTMF信号，规定用8KHz对DTMF信号进行采样，得到时域离散信号

27、：因为采样频率是8000Hz，因此要求每125ms输出一个样本，得到的序列再送到D/A变换器和平滑滤波器，输出便是连续时间的DTMF信号。DTMF信号通过电话线路送到交换机。（2）双音多频信号的检测在接收端，要对收到的双音多频信号进行检测，检测两个正弦波的频率是多少，以判断所对应的十进制数字或者符号。显然这里仍然要用数字方法进行检测，因此要将收到的时间连续 DTMF信号经过A/D变换，变成数字信号进行检测。检测的方法有两种，一种是用一组滤波器提取所关心的频率，根据有输出信号的2个滤波器判断相应的数字或符号。另一种是用DFT（FFT）对双音多频信号进行频谱分析，由信号的幅度谱，判断信号

28、的两个频率，最后确定相应的数字或符号。3、 根据实验原理利用数字信号处理方法进行双音多频信号的产生与解调。clear all;clc; tm=1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68;N=205;K=18,20,22,24,31,34,38,42; f1=697,770,852,941;f2=1209,1336,1477,1633;TN=input('键入10学号=');TNr=0;for l=1:10;d=fix(TN/10(10-l);TN=TN-d*10(10-l); for p=1:4; for q=1:4; if tm(p,q)=ab

29、s(d); break,end end if tm(p,q)=abs(d); break,end end n=0:1023; x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000); sound(x,8000); pause(0.1) X=goertzel(x(1:205),K+1); val = abs(X); subplot(5,2,l); stem(K,val,'.');grid;title('1102303005');xlabel('k');ylabel('|H(k)|');axis(10 50 0 120) limit = 80;