数字信号处理上机实验.docx

1、数字信号处理上机实验 数字信号处理上机实验报告院系: 班级: 学号：姓名： 实验一：信号、系统及系统响应1实验二：用FFT作谱分析7实验三：用双线性变换法设计IIR数字滤波器15四：总结18 实验一：信号、系统及系统响应1.实验目的（1）熟悉联系信号经理想采样后的频谱变化关系，加深对时域采样定理的理解（2）熟悉时域离散系统的时域特性。（3）利用卷积方法观察分析系统的时域特性。（4）掌握序列傅里叶变换的计算机实现方法， 利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。2.实验原理与方法(1) 时域采样。(2) LTI系统的输入输出关系。3.实验内容及步骤(1) 认真复习采样理论、

2、 离散信号与系统、 线性卷积、 序列的傅里叶变换及性质等有关内容， 阅读本实验原理与方法。(2) 编制实验用主程序及相应子程序。 信号产生子程序， 用于产生实验中要用到的下列信号序列： a. xa(t)=A*e-at *sin(0t)u(t)b. 单位脉冲序列：xb(n)=(n)c. 矩形序列： xc(n)=RN(n), N=10 系统单位脉冲响应序列产生子程序。 本实验要用到两种FIR系统。a. ha(n)=R10(n);b. hb(n)=(n)+2.5(n-1)+2.5(n-2)+(n-3) 有限长序列线性卷积子程序用于完成两个给定长度的序列的卷积。 可以直接调用MATLAB语言中的卷积函

3、数conv。 conv用于两个有限长度序列的卷积， 它假定两个序列都从n=0 开始。 调用格式如下：y=conv (x, h)(3)调通并运行实验程序，完成下述实验内容 分析采样序列的特性。采样信号xa(n)的参数为A=444.128,a=50,=50。a. 取采样频率fs=1 kHz,，即T=1 ms。 b. 改变采样频率，fs=300 Hz，观察|X(ej)|的变化，并做记录(打印曲线)；进一步降低采样频率，fs=200 Hz，观察频谱混叠是否明显存在，说明原因， 并记录(打印)这时的|X(ej)|曲线。程序代码如下：function y=x(t,A,a,omiga0)y=A*exp(-a

4、*t).*sin(omiga0*t).*u(t);end%分析采样序列特性A=444.128;a=50*20.5*pi;omiga0=50*20.5*pi;fs1=1000;fs2=300;fs3=200;w=linspace(-2*pi,2*pi,10000);n=0:49;x1=x(n/fs1,A,a,omiga0);x2=x(n/fs2,A,a,omiga0);x3=x(n/fs3,A,a,omiga0);X1=x1*exp(-j*n*w);X2=x2*exp(-j*n*w);X3=x3*exp(-j*n*w);subplot(3,2,1);stem(n,x1,.);ylabel(y);

5、xlabel(n);title(时间函数);subplot(3,2,2);plot(w/pi,abs(X1),r);ylabel(|X(jf)|);xlabel(omega/pi);title(频谱图);text(1.5,1200,f=1000Hz);subplot(3,2,3);stem(n,x2,.);ylabel(y);xlabel(n);title(时间函数);subplot(3,2,4);plot(w/pi,abs(X2),r);ylabel(|X(jf)|);xlabel(omega/pi);title(频谱图);text(1.5,500,f=300Hz);subplot(3,2,

6、5);stem(n,x3,.);ylabel(y);xlabel(n);title(时间函数);subplot(3,2,6);plot(w/pi,abs(X3),r);ylabel(|X(jf)|);xlabel(omega/pi);title(频谱图);text(1.5,250,f=200Hz);图1 对xa不同取样频率得到的取样图和频谱图由图可知，采样频率不同，所得到的采样后信号和其傅里叶变换都不同。fs=1000Hz时，频谱的混叠效应很小，fs=300Hz时，混叠效应加大，fs=200Hz时，混叠效应进一步加大。这是因为采样频率越来越接近临界采样频率，所以会造成频谱的混叠。 时域离散信号

7、、 系统和系统响应分析。 a. 观察信号xb(n)和系统hb(n)的时域和频域特性； 利用线性卷积求信号xb(n)通过系统hb(n)的响应y(n)， 比较所求响应y(n)和hb(n)的时域及频域特性， 注意它们之间有无差别， 绘图说明， 并用所学理论解释所得结果。b. 观察系统ha(n)对信号xc(n)的响应特性。 程序代码如下：a. clear all;xb=double(impact(0,0,1);xc=single(R(1,1,10);ha=single(R(1,1,10);hb=impact(0,0,3)+2.5*impact(1,0,3)+2.5*impact(2,0,3)+impa

8、ct(3,0,3);Xb,w=DFT(xb,2);Hb,w=DFT(hb,4);y=conv(xb,hb);Y,w=DFT(y,5);subplot(3,2,1);stem(0:1,xb,.);ylabel(xb);xlabel(n);title(时间函数);subplot(3,2,2);plot(w/pi,Xb,r);ylabel(|Xb(jf)|);xlabel(omega/pi);title(频谱图);text(1.5,700,f=1000Hz);subplot(3,2,3);stem(hb,.);ylabel(hb);xlabel(n);title(时间函数);subplot(3,2,

9、4);plot(w/pi,Hb,r);ylabel(|Hb(jf)|);xlabel(omega/pi);title(频谱图);text(1.5,700,f=300Hz);subplot(3,2,5);stem(y,.);ylabel(y);xlabel(n);title(时间函数);subplot(3,2,6);plot(w/pi,Y,r);ylabel(|Y(jf)|);xlabel(omega/pi);title(频谱图);text(1.5,700,f=300Hz);图2 xb、hb的时域频域特性(n)的傅里叶变换恒为1。由图可知，xb(n)*hb(n)=hb(n)，即y(n)=hb(n

10、)，所以y(n)和hb(n)的频谱也是完全一样的。b. clear all;xc=single(R(1,1,10);%xc=R10(n)ha=single(R(1,1,10);yn=conv(ha,xc);Y,w=DFT(yn,19);subplot(2,2,1);stem(yn,.);ylabel(yn);xlabel(n);title(时间函数);subplot(2,2,2);plot(w/pi,abs(Y),r);ylabel(|Y(jf)|);xlabel(omega/pi);title(频谱图); xc=single(R(1,1,5);%xc=R5(n)ha=single(R(1,1

11、,10);yn=conv(ha,xc);Y,w=DFT(yn,14);subplot(2,2,3);stem(yn,.);ylabel(yn);xlabel(n);title(时间函数);subplot(2,2,4);plot(w/pi,abs(Y),r);ylabel(|Y(jf)|);xlabel(omega/pi);title(频谱图);图3 xc长度为10和5时与ha的卷积时序图和频域图y(n)的长度与理论计算的相同。两序列卷积后的长度为n1+n2-1，用此方法可以快速验证两个序列的卷积结果是否正确。两次卷积的长度分别为10+10-1=19和10+5-1=14由图可以判断，卷积的结果正

12、确。 卷积定理的验证。A=1,a=0.4,=2.0734，T=1。程序代码如下：n=0:49;xa=x(n,1,0.4,2.0374);hb=impact(0,0,3)+2.5*impact(1,0,3)+2.5*impact(2,0,3)+impact(3,0,3);Xa,w=DFT(xa,50);yn=conv(xa,hb);Y=DFT(yn,53);figure(1);subplot(2,2,1);stem(yn,.);ylabel(yn);xlabel(n);title(时间函数);subplot(2,2,3);plot(w/pi,abs(Y),r);ylabel(|Y(jf)|);x

13、label(omega/pi);title(频谱图);subplot(2,2,4);plot(w/pi,angle(Y),r);ylabel(相位);xlabel(omega/pi);title(Y(jf)的相位);Hb=DFT(hb,4);Yn=Hb.*Xa;figure(2);subplot(1,2,1);plot(w/pi,abs(Yn),r);ylabel(|Yn(jf)|);xlabel(omega/pi);title(频谱图);subplot(1,2,2);plot(w/pi,angle(Yn),r);ylabel(相位);xlabel(omega/pi);title(Yn(jf)

14、的相位);先对xa,hb求卷积再转换为频谱图如下：图4 xa时域图、xa与hb卷积频域图和对应相位图Hb.*Xa的频谱图如下图5 Hb.*Xa的频谱图对比图4图5即可验证卷积定理的正确性。4、思考题(1)在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？它们所对应的模拟频率是否相同？为什么？答：数字频率度量不相同，但他们所对应的模拟频率相同。由公式w=*Ts可知，采样间隔Ts的变化会引起数字频率w的变化，但是不会引起模拟频率的变化。(2)在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶

15、变换，求得的结果有无差异？答：有差异，采样点数不一样，得到的傅里叶变换的点数也会不一样。实验二：用FFT作谱分析1、实验目的(1) 进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法， 所以FFT的运算结果必然满足DFT的基本性质)。(2) 熟悉FFT算法原理和FFT子程序的应用。(3) 学习用FFT对连续信号和时域离散信号进行谱分析的方法， 了解可能出现的分析误差及其原因， 以便在实际中正确应用FFT。2、实验步骤(1) 复习DFT的定义、 性质和用DFT作谱分析的有关内容。(2) 复习FFT算法原理与编程思想， 并对照DIT-FFT运算流图和程序框图， 读懂本实验

16、提供的FFT子程序。(3) 编制信号产生子程序， 产生以下典型信号供谱分析用： (4) 编写主程序。(5) 按实验内容要求，上机实验， 并写出实验报告。3、实验内容(1) 对 2 中所给出的信号逐个进行谱分析。程序代码如下：x1=1 1 1 1 0 0 0 0 ;x2=1 2 3 4 4 3 2 1;x3=4 3 2 1 1 2 3 4;fs=64;N1=8;N2=16;N3=32;N4=64;% x4=cos(0.25*pi*n);% x5=sin(pi/8)*n);% x7=cos(0.25*pi*n)+sin(pi/8)*n);% x8=cos(0.25*pi*n)+j*sin(pi/8

17、)*n);%x1作图i=1;f1=fft(x1,N1);figure(i);i=i+1;subplot(1,3,1);stem(0:N1-1,x1,.);title(x1n的波形);ylabel(x1n);xlabel(n);subplot(1,3,2);stem(0:N1-1,abs(f1),.);title(x1n的8点FFT);ylabel(|X(k)|);xlabel(k);f1=fft(x1,N2);subplot(1,3,3);stem(0:N2-1,abs(f1),.);title(x1n的16点FFT);ylabel(|X(k)|);xlabel(k);%x2作图f2=fft(

18、x2,N1);figure(i);i=i+1;subplot(1,3,1);stem(0:N1-1,x2,.);title(x2n的波形);ylabel(x2n);xlabel(n);subplot(1,3,2);stem(0:N1-1,abs(f2),.);title(x2n的8点FFT);ylabel(|X(k)|);xlabel(k);f2=fft(x2,N2);subplot(1,3,3);stem(0:N2-1,abs(f2),.);title(x2n的16点FFT);ylabel(|X(k)|);xlabel(k);%x3作图f3=fft(x3,N1);figure(i);i=i+

19、1;subplot(1,3,1);stem(0:N1-1,x3,.);title(x3n的波形);ylabel(x3n);xlabel(n);subplot(1,3,2);stem(0:N1-1,abs(f3),.);title(x3n的8点FFT);ylabel(|X(k)|);xlabel(k);f3=fft(x3,N2);subplot(1,3,3);stem(0:N2-1,abs(f3),.);title(x3n的16点FFT);ylabel(|X(k)|);xlabel(k);%x4作图n=0:N1-1;x4=sin(pi/8)*n);f4=fft(x4,N1);figure(i);

20、i=i+1;subplot(2,2,1);stem(0:N1-1,x4,.);title(x4n的8点波形);ylabel(x4n);xlabel(n);subplot(2,2,2);stem(0:N1-1,abs(f4),.);title(x4n的8点FFT);ylabel(|X(k)|);xlabel(k);n=0:N2-1;x4=cos(0.25*pi*n);f4=fft(x4,N2);subplot(2,2,3);stem(0:N2-1,x4,.);title(x4n的16点波形);ylabel(x4n);xlabel(n);subplot(2,2,4);stem(0:N2-1,abs

21、(f4),.);title(x4n的16点FFT);ylabel(|X(k)|);xlabel(k);%x5作图figure(i);i=i+1;n=0:N1-1;x5=cos(0.25*pi*n);f5=fft(x5,N1);subplot(2,2,1);stem(0:N1-1,x5,.);title(x5n的8点波形);ylabel(x5n);xlabel(n);subplot(2,2,2);stem(0:N1-1,abs(f5),.);title(x5n的8点FFT);ylabel(|X(k)|);xlabel(k);n=0:N2-1;x5=cos(0.25*pi*n);f5=fft(x5

22、,N2);subplot(2,2,3);stem(0:N2-1,x5,.);title(x5n的16点波形);ylabel(x5n);xlabel(n);subplot(2,2,4);stem(0:N2-1,abs(f5),.);title(x5n的16点FFT);ylabel(|X(k)|);xlabel(k);%x6作图figure(i);i=i+1;n=0:N2-1;x=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);f6=fft(x,N2);subplot(3,2,1);stem(0:N2-1,x,.);title(x6n的16点波形);y

23、label(x6n);xlabel(n);subplot(3,2,2);stem(0:N2-1,abs(f6),.);title(x6n的16点FFT);ylabel(|X(k)|);xlabel(k);n=0:N3-1;x=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);f6=fft(x,N3);subplot(3,2,3);stem(0:N3-1,x,.);title(x6n的32点波形);ylabel(x6n);xlabel(n);subplot(3,2,4);stem(0:N3-1,abs(f6),.);title(x6n的32点FFT);

24、ylabel(|X(k)|);xlabel(k);n=0:N4-1;x=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);f6=fft(x,N4);subplot(3,2,5);stem(0:N4-1,x,.);title(x6n的64点波形);ylabel(x6n);xlabel(n);subplot(3,2,6);stem(0:N4-1,abs(f6),.);title(x6n的64点FFT);ylabel(|X(k)|);xlabel(k);运行结果如下:图6 x1的8、16点谱分析图7 x2的8、16点谱分析图8 x3的8、16点谱分析图9

25、x4的8、16点谱分析图10 x5的8、16点谱分析图11 x6的16、32、64点谱分析(2)令x7(n)=x4(n)+x5(n)，用FFT计算8点和16点离散傅里叶变换。程序代码如下：%x7作图figure(i);i=i+1;n=0:N1-1;x7=sin(pi/8)*n)+cos(0.25*pi*n);f7=fft(x7,N1);t1=max(x7);t2=max(f7);subplot(3,2,1);stem(0:N1-1,x7,.);title(x7n的8点波形);ylabel(x7n);xlabel(n);subplot(3,2,2);stem(0:N1-1,abs(f7),.);

26、title(x7n的8点FFT);ylabel(|X(k)|);xlabel(k);n=0:N2-1;x7=sin(pi/8)*n)+cos(0.25*pi*n);f7=fft(x7,N2);t1=max(x7);t2=max(f7);subplot(3,2,3);stem(0:N2-1,x7,.);title(x7n的16点波形);ylabel(x7n);xlabel(n);subplot(3,2,4);stem(0:N2-1,abs(f7),.);title(x7n的16点FFT);ylabel(|X(k)|);xlabel(k);k=conj(f7);x4=(k+f7)/2;subplo

27、t(3,2,5);stem(0:N2-1,abs(x4),.);title(恢复的X4(K);ylabel(|Re(X7(k)|);xlabel(k);x5=(k-f7)/2;subplot(3,2,6);stem(0:N2-1,abs(x5),.);title(恢复的X5(K);ylabel(|jIm(X7(k)|);xlabel(k);运行结果图12 左1：x7的8点时域图 右1：x7的8点频谱图 左2：x7的16点时域图 右2：x7的16点频谱图 左3：对ifft(Re(X7) 的到的X4(k)频谱 右3：对ifft(Im(X7) 的到的X5(k)频谱 将上图 左3和右3 与图4图5中的

28、16点频谱图对比可知两者对应完全相等，验证了DFT的如下性质：xnep Re(Xn) xnop jIm(Xn)。(3)令x8(n)=x4(n)+j*x5(n)，重复(2)。程序代码如下：%x8作图figure(i);i=i+1;n=0:N1-1;x8=cos(0.25*pi*n)+j*sin(pi/8)*n);f8=fft(x8,N1);subplot(3,2,1);stem(0:N1-1,x8,.);title(x8n的8点波形);ylabel(x8n);xlabel(n);subplot(3,2,2);stem(0:N1-1,abs(f8),.);title(x8n的8点FFT);ylabel(|X(k)|);xlabel(k);n=0:N2-1;x8=cos(0.25*pi*n)+j*sin(pi/8)*n);f8=fft(x8,N2);subplot(3,2,3);stem(0:N2-1,x8,.);title(x8n的16点波形);ylabel(x8n);xlabel(n);subplot(3,2,4);stem(0:N2-1,abs(f8),.);title(x8n的16点FFT);ylabel(|X(k)|);xlabel(k);k(1)=conj(f8(1);for m=2:N2k(