数字信号处理仿真实验

1、数字信号处理上机实验仿真报告学 院 电子工程学院 班 级 学 号 姓 名 实验一： 信号、 系统及系统响应一、实验目的(1) 熟悉连续信号经理想采样前后的频谱变化关系， 加深对时域采样定理的理解。(2) 熟悉时域离散系统的时域特性。(3) 利用卷积方法观察分析系统的时域特性。(4) 掌握序列傅里叶变换的计算机实现方法， 利用序列的傅里叶变换对连续信号、 离散信号及系统响应进行频域分析。 二、实验原理采样是连续信号数字处理的第一个关键环节。对一个连续信号进行理想采样的过程可用(1.1)式表示。 (1.1)其中为的理想采样，为周期冲激脉冲，即 (1.2)的傅里叶变换为 (1.3)将(1.2)式代入

2、(1.1)式并进行傅里叶变换， (1.4)式中的就是采样后得到的序列， 即的傅里叶变换为 (1.5)比较(1.5)和(1.4)可知 (1.6)为了在数字计算机上观察分析各种序列的频域特性，通常对在上进行m点采样来观察分析。对长度为n的有限长序列，有 (1.7)其中一个时域离散线性时不变系统的输入/输出关系为 (1.8)上述卷积运算也可以转到频域实现 (1.9)三、 实验内容及步骤 (1) 认真复习采样理论、 离散信号与系统、 线性卷积、 序列的傅里叶变换及性质等有关内容， 阅读本实验原理与方法。 (2) 编制实验用主程序及相应子程序。 信号产生子程序， 用于产生实验中要用到的下列信号序列： x

3、a(t)=ae-at sin(0t)u(t)进行采样， 可得到采样序列 xa(n)=xa(nt)=ae-ant sin(0nt)u(n), 0n=0);%产生脉冲信号functionx,n=maichong(n0,n1,n2)n=(n1:n2);x=(n=n0);2.主程序 分析采样序列的特性程序代码如下（见后附程序ex1_1）:a=100;a=200;w0=200;k=-200:200;t=0.001;t=0:t:0.06;n=50;k1=0:1:n;w1max=2*pi*500;w1=w1max*k1/n;w1=w1/pi;xat=a*exp(-a*t).*sin(w0*t).*u(t);

4、xa=xat*exp(-j*t*w1);figure(1);subplot(2,1,1);plot(t,xat);xlabel(t);ylabel(xa(t);title(采样信号a(t);axis(0,0.06,-5,35);subplot(2,1,2);plot(w1,abs(xa);xlabel(w);ylabel(x(jw);title(xa(t)的频谱);%fs=1000hza=100;a=200;w0=200;k=-200:200;fs=1000; w=k/50;xan = xn(a,a,w0,fs); %产生信号xa(n)x=dft(xan,50);figure(2)subplo

5、t(2,1,1)n=0:49;stem(n,xan,.);axis(0,50,-20,50);xlabel(n);ylabel(xa(n);title(采样信号fs=1000hz);subplot(2,1,2);plot(w,abs(x);xlabel(w/pi);ylabel(x(ejw);title(xa(n)的频谱(fs=1000hz);%fs=300hzfs=300; xan = xn(a,a,w0,fs); %产生信号xa(n)x=dft(xan,50);figure(3)subplot(2,1,1)n=0:49;stem(n,xan,.);axis(0,50,-20,50);xla

6、bel(n);ylabel(xa(n);title(采样信号fs=300hz);subplot(2,1,2);plot(w,abs(x);xlabel(w/pi);ylabel(x(ejw);title(xa(n)的频谱(fs=300hz);%fs=200hzfs=200; xan = xn(a,a,w0,fs); %产生信号xa(n)x=dft(xan,50);figure(4)subplot(2,1,1)n=0:49;stem(n,xan,.);axis(0,50,-20,50);xlabel(n);ylabel(xa(n);title(采样信号fs=200hz);subplot(2,1,

7、2);plot(w,abs(x);xlabel(w/pi);ylabel(x(ejw);title(xa(n)的频谱(fs=200hz);仿真结果如下:由图可见，在折叠频率w=,即f=fs/2=500hz处混叠很小。当fs=300hz时，存在较明显的混叠失真；当fs=200时，发生严重的混叠失真。时域离散信号、系统和系统响应分析a. 观察信号xb(n)和系统hb(n)的时域和频域特性；实验代码如下：（见后附程序ex1_2）%信号xbn=64xb=zeros(1,n)xb(1)=1;xn1=0:n-1figure(1);subplot(2,1,1)stem(xn1,xb,.)title(信号xb

8、时域特性)w=linspace(-5*pi,5*pi,1000)xb=xb*exp(-j*xn1*w)figure(1)subplot(2,1,2)plot(w/pi,abs(xb),.)title(信号xb频域特性)%信号hbhb=1,2.5,2.5,1;nhb=0:3figure(2);subplot(2,1,1)stem(nhb,hb,.)title(系统函数hb时域特性)w=linspace(-2.8*pi,2.8*pi,1000)hb=hb*exp(-j*nhb*w)figure(2)subplot(2,1,2)plot(w/pi,abs(hb),.)title(系统函数hb频域特性

9、)%响应yb=conv(xb,hb);nyb=0:length(yb)-1figure(3);subplot(2,1,1)stem(yb,.)title(响应yb时域特性)w=linspace(-2.8*pi,2.8*pi,1000)yb=yb*exp(-j*nyb*w)figure(3)subplot(2,1,2)plot(w/pi,abs(yb),.)title(响应yb频域特性)仿真结果如下:b. 观察系统ha(n)对信号xc(n)的响应特性。实验代码如下(见后附程序ex1_3):%信号xcn=10xc=ones(1,n)xn1=0:n-1figure(1);subplot(2,1,1)

10、stem(xn1,xc,.)title(信号xc时域响应)w=linspace(-5*pi,5*pi,1000)xc=xc*exp(-j*xn1*w)figure(1)subplot(2,1,2)plot(w/pi,abs(xc)%信号han=10ha=ones(1,n)hn1=0:n-1figure(2);subplot(2,1,1)stem(hn1,ha,.)title(系统函数ha时域响应)w=linspace(-2.8*pi,2.8*pi,1000)ha=ha*exp(-j*hn1*w)figure(2)subplot(2,1,2)plot(w/pi,abs(ha),.)title(系

11、统函数ha频域响应)%响应y2=conv(xc,ha);ny2=0:length(y2)-1figure(3);subplot(2,1,1)stem(y2,.)title(响应y2时域响应)w=linspace(-2.8*pi,2.8*pi,1000)y2=y2*exp(-j*ny2*w)figure(3)subplot(2,1,2)plot(w/pi,abs(y2),.)title(响应yb频域响应)仿真结果如下: 卷积定理的验证。 实验代码如下：（见后附程序ex1_4）:%时域卷积定理的验证%信号xcn=10xc=ones(1,n)xn1=0:n-1w=linspace(-5*pi,5*p

12、i,1000)xc=xc*exp(-j*xn1*w)%系统信号han=10ha=ones(1,n)hn1=0:n-1w=linspace(-5*pi,5*pi,1000)ha=ha*exp(-j*hn1*w)%响应y2=conv(xc,ha);ny2=0:length(y2)-1w=linspace(-5*pi,5*pi,1000)y2=y2*exp(-j*ny2*w)figure(1)subplot(2,1,1)plot(w/pi,abs(y2)title(响应yb=dtft(conv(xc,ha)频域特性)y=xc.*hasubplot(2,1,2)plot(w/pi,abs(y)titl

13、e(响应y=xc.*ha频域特性)仿真结果如下：两种方法得到的谱是一样的，即验证了卷积定理。五、思考题(1) 在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同? 它们所对应的模拟频率是否相同? 为什么?答：由可知，若采样频率不同，则其周期t不同，相应的数字频率也不相同；而因为是同一信号，故其模拟频率保持不变。(2) 在卷积定理验证的实验中， 如果选用不同的频域采样点数m值， 例如， 选m=10和m=20， 分别做序列的傅里叶变换， 求得所得结果之间有无差异? 为什么?答：有差异。因为所得图形由其采样点数唯一确定，由频域采样定理可知，若m小

14、于采样序列的长度n，则恢复原序列时会发生时域混叠现象。实验二：用fft作谱分析一、 实验目的：（1） 进一步加深dft算法原理和基本性质的理解（因为fft只是dft的一种快速算法，所以fft的运算结果必然满足dft的基本性质）。（2） 熟悉fft算法原理和fft子程序的应用。（3） 学习用fft对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用fft。二、 实验步骤：（1） 复习dft的定义、性质和用dft作谱分析的有关内容。（2） 复习fft算法原理与编程思想，并对照ditfft运算流图和程序框图，读懂本实验提供的fft子程序。（3） 编制信号产生

15、子程序，产生以下典型信号供谱分析用：=r4(n)=cos(pi/4*n)=sin(pi/8*n)=cos(pi*8*t)+ sin(pi*16*t)+cos(20*pi*t)应当注意，如果给出的是连续信号xa(t),则首先要根据其最高频率确定采样速率fs以及由频率选择采样点数n，然后对其进行软件采样（即计算x(n)=xa(nt),0=n=n-1）,产生对应序列x(n)。对信号x6(t)，频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。对周期序列，最好截取周期的整数倍进行分析，否则有可能产生较大的分析误差。请实验者根据dft的隐含周期性思考这个问题。（4） 编写主程序。下图给出了主程

16、序框图，供参考。本实验提供fft子程序和通用绘图子程序。开始读入长度n调用信号产生子程序产生试验信号调用fft子程序（函数）计算信号的dtf调用绘图子程序（函数）绘制x（k）曲线调用绘图子程序（函数）绘制时间序列波形图结束三实验程序及对应波形(1) 对 2 中所给出的信号逐个进行谱分析。 实验代码如下（见后附程序ex2_1） %信号x1x1=1,1,1,1n1=0:3figure(1);subplot(3,1,1)stem(n1,x1)title(信号x1时域特性)w=linspace(-5*pi,5*pi,1000)x1=x1*exp(-j*n1*w)figure(1)subplot(3,1

17、,2)plot(w/pi,abs(x1)title(dtft分析信号x1频域特性)x11=fftshift(fft(x1,128)figure(1)subplot(3,1,3)stem(abs(x11)title(fft分析信号x1频域特性)%信号x2x2=1,2,3,4,4,3,2,1n2=0:7figure(2);subplot(3,1,1)stem(n2,x2)title(信号x2时域特性)w=linspace(-5*pi,5*pi,1000)x2=x2*exp(-j*n2*w)figure(2)subplot(3,1,2)plot(w/pi,abs(x2)title(dtft分析信号x

18、2频域特性)x22=fftshift(fft(x2,128)figure(2)subplot(3,1,3)stem(abs(x22)title(fft分析信号x2频域特性)%信号x3x3=4 3 2 1 1 2 3 4n3=0:7figure(3);subplot(3,1,1)stem(n3,x3)title(信号x3时域特性)w=linspace(-5*pi,5*pi,1000)x3=x3*exp(-j*n3*w)figure(3)subplot(3,1,2)plot(w/pi,abs(x3)title(dtft信号x3频域特性)x33=fftshift(fft(x3,128)figure(

19、3)subplot(3,1,3)stem(abs(x3)title(fft分析信号x3频域特性)%信号x4=cos(pi/4*n)n4=0:15x4=cos(pi/4*n4)figure(4);subplot(3,1,1)stem(n4,x4)title(信号x4时域特性)w=linspace(-5*pi,5*pi,1000)x4=x4*exp(-j*n4*w)figure(4)subplot(3,1,2)plot(w/pi,abs(x4)title(dtft信号x4频域特性)x44=fftshift(fft(x4,128)figure(4)subplot(3,1,3)stem(abs(x44

20、)title(fft分析信号x4频域特性)%信号x5=sin(pi/8*n)n5=0:15x5=sin(pi/8*n5)figure(5);subplot(3,1,1)stem(n5,x5)title(信号x5时域特性)w=linspace(-5*pi,5*pi,1000)x5=x5*exp(-j*n5*w)figure(5)subplot(3,1,2)plot(w/pi,abs(x5)x55=fftshift(fft(x5,128)figure(5)subplot(3,1,3)stem(abs(x55)title(fft分析信号x4频域特性)%信号x6=cos(8*pi*t)+cos(16*

21、pi*t)+cos(20*pi*t)t=0:0.0001:2.5xa=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t) %原始的连续时间信号xafigure(6);subplot(3,1,1);plot(t,xa);title(连续时间信号曲线xa); ;axis(0 2.5 -5 5) %画连续时间信号曲线t=0.01;n6=0:2.5/t; %建立离散自变量向量x6=cos(8*pi*n6*t)+cos(16*pi*n6*t)+cos(20*pi*n6*t) %采样周期为t的离散时间信号xfigure(6);subplot(3,1,2);stem(n6*t,x6)

22、;title(取样信号);axis(0 2.5 -5 5) %画离散时间信号曲线w=linspace(-5*pi,5*pi,2500); %均匀取点x6=x6*exp(-j*n6*w) %计算dtftfigure(6);subplot(3,1,3);plot(w/pi,abs(x6);title(dtft分析取样信号频谱);仿真结果如下：(2) 令x(n)=x4(n)+x5(n)， 用fft计算 8 点和 16 点离散傅里叶变换，x(k)=dftx(n)(3) 令x(n)=x4(n)+jx5(n)， 重复(2)。代码如下：（见后附程序ex2_2）%信号x4=cos(pi/4*n)n4=0:15

23、x4=cos(pi/4*n4)figure(1);subplot(2,1,1)stem(n4,x4)title(信号x4时域特性)w=linspace(-5*pi,5*pi,1000)x4=x4*exp(-j*n4*w)figure(1)subplot(2,1,2)plot(w/pi,abs(x4)title(dtft分析信号x4频域特性)%信号x5=sin(pi/8*n)n5=0:15x5=sin(pi/8*n5)figure(2);subplot(2,1,1)stem(n5,x5)title(dtft分析信号x5时域特性)w=linspace(-5*pi,5*pi,1000)x5=x5*e

24、xp(-j*n5*w)figure(2)subplot(2,1,2)plot(w/pi,abs(x5)%x(n)=x4(n)+x5(n)xn1=x4+x5nx1=0:length(xn1)-1figure(3);subplot(3,1,1)stem(nx1,xn1)title(xn1=x4+x5时域特性)xn1_8=fftshift(fft(xn1,8)n1_8=0:7figure(3);subplot(3,1,2)stem(n1_8,abs(xn1_8)title(xn1的8点fft变换)xn1_16=fftshift(fft(xn1,16)n1_16=0:15figure(3);subpl

25、ot(3,1,3)stem(n1_16,abs(xn1_16)title(xn1的16点fft变换)%x(n)=x4(n)+jx5(n)xn2=x4+j*x5nx2=0:length(xn2)-1figure(4);subplot(3,1,1)stem(nx2,xn2)title(xn2=x4+j*x5时域特性)xn2_8=fftshift(fft(xn2,8)n2_8=0:7figure(4);subplot(3,1,2)stem(n2_8,abs(xn2_8)title(xn2的8点fft变换)xn2_16=fftshift(fft(xn2,16)n2_16=0:15figure(4);s

26、ubplot(3,1,3)stem(n2_16,abs(xn2_16)title(xn2的16点fft变换)仿真结果如下：四、思考题(1) 在n=8时， x2(n)和x3(n)的幅频特性会相同吗? 为什么? n=16呢?答：在n=8时，x2(n)和x3(n)的幅频特性会相同在n=8时 x2(n)和x3(n)的幅频特性不相同。(2) 如果周期信号的周期预先不知道， 如何用fft进行谱分析?答：如果周期信号的周期预先不知道。实验三： 用窗函数法设计fir数字滤波器一、 实验目的(1) 掌握用窗函数法设计fir数字滤波器的原理和方法。(2) 熟悉线性相位fir数字滤波器特性。(3) 了解各种窗函数对

27、滤波特性的影响。二、实验原理如果所希望的滤波器的理想频率响应函数为 ，则其对应的单位脉冲响应为= (2-1)窗函数设计法的基本原理是用有限长单位脉冲响应序列逼近。由于往往是无限长序列，且是非因果的，所以用窗函数将截断，并进行加权处理，得到： (2-2)就作为实际设计的fir数字滤波器的单位脉冲响应序列，其频率响应函数为： (2-3)式中，n为所选窗函数的长度。 由第七章可知，用窗函数法设计的滤波器性能取决于窗函数的类型及窗口长度n的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度n。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见第七章。 这样选定窗函数

28、类型和长度n后，求出单位脉冲响应，并按式（2-3）求出。是否满足要求，要进行验算。一般在尾部加零使长度满足于2的整数次幂，以便用fft计算。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度n，再次验算，直至满足要求。如果要求线性相位特性，则还必须满足 (2-4)根据上式中的正负号和长度n的奇偶性又将线性相位fir滤波器分成四类。要根据设计的滤波特性正确选择其中一类。例如，要设计线性低通特征，可选择一类，而不能选一类.三、实验内容及步骤(1) 复习用窗函数法设计fir数字滤波器一节内容， 阅读本实验原理， 掌握设计步骤。(2) 编写程序。 编写能产生矩型窗、 升余弦

29、窗、 改进升余弦窗和二阶升余弦窗的窗函数子程序。 编写主程序。 其中幅度特性要求用db表示。三、 实验代码及仿真结果代码如下(见后附程序ex3_1):fcp=25000;fcs=10000;f=fcp-fcs;fs=50000;f=fcp-fcs;b=2*pi*f/fs; n=ceil(11*pi/b); wc=(fcp+fcs)/fs; h=fir1(n,wc,blackman(n+1);n=0:length(h)-1;w=linspace(-pi,pi,10000); h1=h*exp(-j*n*w); h=20*log10(abs(h1); figure(1)stem(n,h,.) xl

30、abel(n);ylabel(h(n)title(h(n)的时域特性曲线) figure(2)subplot(2,1,1)plot(w/pi,h); xlabel(omega/);ylabel(|h(ejomega)|/db)title(h(n)的幅度谱);subplot(2,1,2)plot(w/pi,angle(h1/pi);xlabel(omega/);ylabel(相位) title(h(n)的相位谱);仿真结果如下:以下程序见后附程序(ex3_2)n=input(窗宽度n=); k=input(窗型1矩形窗2.升余弦窗3.改进的升余弦窗4.二阶升余弦窗请选择); subplot(2,

31、2,1); w=pi/5; a=(n-1)/2; n=0:(n-1); m=n-a+eps; h=sin(w*m)./(pi*m); if k=1 b=bartlett(n); else if k=2 b=hanning(n); else if k=3 b=hamming(n); else if k=4 b=blackman(n); end end endendhd=h.*(b);stem(n,hd,.); xlabel(n); ylabel(h(n); title(在矩形窗下的n=64时h(n)函数); subplot(2,2,2); h,m=freqz(hd,1,1024); mag=ab

32、s(h); db=20*log10(mag+eps)/max(mag); plot(m/pi,db); xlabel(w/n); ylabel(20logh(ejw); title(h(n)的幅频特性); grid; pha=angle(h); subplot(2,2,3); plot(m,pha); xlabel(n); ylabel(c); title(h(n)的相频特性); subplot(2,2,4); plot(m,mag); xlabel(w); ylabel(h(ejw); title(h(n)的幅频特性);仿真结果如下:五、思考题(1) 如果给定通带截止频率和阻带截止频率以及阻

33、带最小衰减， 如何用窗函数法设计线性相位低通滤波器? 写出设计步骤。答： 首先确定模拟带通滤波器的技术指标。 确定归一化低通技术要求 设计归一化低通 将低通转化为带通(2) 如果要求用窗函数法设计带通滤波器， 且给定上、 下边带截止频率为1和2，试求理想带通的单位脉冲响应hd(n)。 答： 首先确定所需数字滤波器的类型，确定线性相位数字滤波器的类型。 选择合适的窗函数确定理想数字滤波器的频率响应函数计算理想带通的单位脉冲响应hd(n)实验四：用双线性变换法设计iir数字滤波器一、实验目的(1) 掌握用双线性变换法设计iir数字滤波器的原理和方法。(2) 熟悉iir数字滤波器特性。二、实验内容及

34、步骤(1) 复习用双线性变换法设计iir数字滤波器一节内容， 阅读本实验原理， 掌握设计步骤。(2) 编写程序。 设计产生巴特沃斯滤波器原型，采用双线性变换法设计巴特沃斯数字低通滤波器。 编写主程序。 其中幅度特性要求用db表示。三、上机实验内容 (1) 实验所需设计的iir数字低通滤波器的指标为wp=0.2*pi，ws=0.3*pi，rp=1db，as=15db，滤波器采样频率为100hz。 (2) 采用巴特沃斯滤波器为原型进行设计。(3) 将相同指标要求的iir滤波器和fir滤波器进行比较。四、实验代码及仿真结果(见后附程序ex4_1)fs=100; %采样频率t=1/fs;wp=0.2*

35、pi; ws=0.3*pi; %设置归一化通带和阻带截止频率ap=1;as=15; %设置通带最大和最小衰减fs=fs/fs;n,wc=buttord(wp,ws,ap,as,s); %调用butter函数确定巴特沃斯滤波器阶数b,a=butter(n,wc,s); %调用butter函数设计巴特沃斯滤波器w=linspace(0,pi,400*pi); %指定一段频率值 hf=freqs(b,a,w); %计算模拟滤波器的幅频响应 hf=20*log10(abs(hf);subplot(2,1,1);plot(w/pi,hf); %绘出巴特沃斯模拟滤波器的幅频特性曲线 grid on;title(iir巴特沃斯模拟滤波器);xlabel(频率/hz);ylabel(幅度);d,c= bilinear (b,a,fs); %调用双线性变换法 hz=freqz(d,c,w); %返回频率响应 hz=20*log10(abs(hz