数字信号处理实验报告

山东建筑大学实验报告学院: 信息与电气工程学院 班级: 姓名: 学号: 课程: 数字信号处理 实验日期: 2011 年 10 月 30日 成绩: 实验一 信号、系统及系统响应1. 实验目的(1) 熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。(2) 熟悉时域离散系统的时域特性。(3) 利用卷积方法观察分析系统的时域特性。(4) 掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。 2. 实验原理与方法采样是连续信号数字处理的第一个关键环节。对一个连续信号进行理想采样的过程可用(1.1)式表示。 (1.1)其中为的理想采样，为周期冲激脉冲，即 (1.2)的傅里叶变换为 (1.3)将(1.2)式代入(1.1)式并进行傅里叶变换， (1.4)式中的就是采样后得到的序列， 即的傅里叶变换为 (1.5)比较(1.5)和(1.4)可知 (1.6)为了在数字计算机上观察分析各种序列的频域特性，通常对在上进行M点采样来观察分析。对长度为N的有限长序列，有 (1.7)其中一个时域离散线性非移变系统的输入/输出关系为 (1.8)上述卷积运算也可以在频域实现 (1.9)3. 实验内容及步骤(1) 认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容，阅读本实验原理与方法。(2) 编制实验用主程序及相应子程序。 信号产生子程序， 用于产生实验中要用到的下列信号序列：a. 采样信号序列：对下面连续信号：进行采样， 可得到采样序列： ， 其中A为幅度因子，a为衰减因子，0是模拟角频率，T为采样间隔。这些参数都要在实验过程中由键盘输入，产生不同的和。b. 单位脉冲序列：c. 矩形序列： 系统单位脉冲响应序列产生子程序。本实验要用到两种FIR系统。a. b. 有限长序列线性卷积子程序，用于完成两个给定长度的序列的卷积。 可以直接调用MATLAB语言中的卷积函数conv。 conv用于两个有限长度序列的卷积，它假定两个序列都从n=0 开始。调用格式如下： y=conv (x, h)其中参数x和h是两个已赋值的行向量序列。 图1.1 的幅频特性曲线 (3) 运行实验程序，完成下述实验内容： 分析采样序列的特性。产生采样信号序列，使A=444.128， ，。图1.1给出了连续信号的幅频特性曲线。a. 取采样频率kHz，即T=1ms。观察所得采样的幅频特性和图1.1中的，在折叠频率附近有无明显差别。b. 改变采样频率，Hz，观察的变化，并做记录(打印曲线)； 进一步降低采样频率，Hz，观察频谱混叠是否明显存在，说明原因，并记录(打印)这时的曲线。 时域离散信号、系统和系统响应分析。a. 观察信号和系统的时域和频域特性；利用线性卷积求信号通过系统的响应，比较所求响应和的时域及频域特性，注意它们之间有无差别，绘图说明，并用所学理论解释所得结果。b. 观察系统对信号的响应特性。利用线性卷积求系统响应，并判断图形及其非零值序列长度是否与理论结果一致，对，说出一种定性判断图形正确与否的方法。调用序列傅里叶变换数值计算子程序，求得，观察特性曲线，定性判断结果的正确性。改变的长度，取，重复该实验。注意参数变化的影响，说明变化前后的差异，并解释所得结果。 卷积定理的验证。将实验中的信号换成，使，重复实验a，打印曲线；对主程序作简单修改，按式（1.9）计算，并绘出曲线，与前面直接对进行傅里叶变换所得幅频特性曲线进行比较，验证时域卷积定理。实验源代码：产生采样信号序列，使A=444.128， ，A=444.128;a=50*2(1/2)*pi;m=a;N=50;fs=1000;n=0:N-1;t=n/fs;x=A*exp(-a*t).*sin(m*t);y=fft(x,N);f=n*fs/N;stem(f,abs(y),.);xlabel(f/Hz);ylabel(幅度);grid on;由图像可看出，与图1.1中的相比，在折叠频率附近无明显差别。Fs=300Hz时，结果为：Fs =200Hz时，结果为：分析：由以上三个图可明显看出当fs为1000Hz、300Hz及200Hz时没有在Fs/2处发生频谱混叠，因为三个抽样频率均大于原信号频率的两倍，满足抽样定理，因此 时域离散信号、系统和系统响应分析。ai观察信号和系统的时域和频域特性N=32;n=0:50;xb=1 zeros(1,50);y=fft(xb,N);subplot(1,2,1);stem(n,xb,.);grid on;title(xb时域特性);subplot(1,2,2);stem(0:length(y)-1,y,.);title(xb频域特性);grid on;结果：N=50;n=0:N-1;xb=1 2.5 2.5 1 zeros(1,46);y=fft(xb,N);subplot(2,1,1);stem(n,xb,.);grid on;title(hb时域特性);subplot(2,1,2);stem(0:length(y)-1,y,.);title(hb频域特性);grid on;结果：ii利用线性卷积求信号通过系统的响应N=100;n=0:50;xb=1 zeros(1,50);hb=1 2.5 2.5 1 zeros(1,47);yn=conv(xb,hb)subplot(2,1,1);stem(0:length(yn)-1,yn,.);grid on;title(yn时域特性);subplot(2,1,2);y=fft(yn,N);stem(0:length(y)-1,y,.);title(yn频域特性);grid on;iii比较所求响应和的时域及频域特性，注意它们之间有无差别，绘图说明，并用所学理论解释所得结果由以上图可看出yn和hb的时域和频域特性均相同，这是由于yb是冲激函数，任何信号经过yb后保持不变，因此yn和hb的时域和频域特性均相同，因为这是两个相同的信号。bi观察系统对信号的响应特性。xcn=ones(1,10) zeros(1,40);han=xcn;y=conv(xcn,han);stem(0:length(y)-1,y,.);title(系统响应y(n)时域特性);grid on;结果：经计算，图形及其非零值序列长度与理论结果一致。ii说出一种定性判断图形正确与否的方法：的序列长度为N1+N2-1,N1、N2分别为和的序列长度，且不为零的部分关于n=n1成偶对称分布。iii调用序列傅里叶变换数值计算子程序，求得，观察特性曲线，定性判断结果的正确性。N=100;xcn=ones(1,10) zeros(1,40);han=xcn;y=conv(xcn,han);stem(0:length(y)-1,y,.);title(系统响应y(n)时域特性);grid on;y2=fft(y,N);stem(0:length(y2)-1,abs(y),.);title(系统响应y(n)的频域幅频特性);grid on;结果：iv当N=5时，y(n)的时域特性曲线为：xcn=ones(1,5) zeros(1,45);han=ones(1,10) zeros(1,40);y=conv(xcn,han);stem(0:length(y)-1,y,.);title(系统响应y(n)时域特性);grid on; 频域特性曲线为：xcn=ones(1,5) zeros(1,45);han=ones(1,10) zeros(1,40);y=conv(xcn,han);stem(0:length(y)-1,y,.);title(系统响应y(n)时域特性);grid on;y2=fft(y,N);stem(0:length(y2)-1,abs.);title(系统响应y(n)的频域幅频特性);grid on;分析：由图可见时域特性和频域特性均有所变化，但整体特性没变。这是由于当N变化时，所对应的序列卷积时对应项相乘时发生变化，因此y(n)的非零值所对应的n及非零值的幅度都会相应改变，但扔保持轴对称特性，这是有序列本身决定的。而对于频域特性，N的变化导致的频谱Sa函数发生变化，从而=发生改变，但总体的趋势不会不会改变，这是由两个Sa函数的乘积所决定的。 卷积定理的验证。程序：A=1;a=0.4;m=2.0734;N=50;fs=1;n=0:N-1;t=n/fs;x=A*exp(-a*t).*sin(m*t);y=fft(x,N);f=n*fs/N;subplot(2,1,1);stem(0:length(x)-1,x,.);title(xa(n)的时域特性);grid on;subplot(2,1,2);stem(f,abs(y),.);grid on;xlabel(f/Hz);ylabel(幅度);title(xa(n)的频域幅频特性);subplot(3,1,3);stem(f,angle(y),.);xlabel(f/Hz);ylabel(相位);title(xa(n)的频域相频特性);grid on;y(n)的时域特性：程序：A=1;a=0.4;m=2.0734;N=50;fs=1;n=0:N-1;t=n/fs;x=A*exp(-a*t).*sin(m*t);hb=1 2.5 2.5 1 zeros(1,47);m=conv(hb,x);subplot(3,1,1);stem(0:length(m)-1,m,.);title(h(n)的时域特性);grid on;y(n)的频域特性y=fft(m,2*N);subplot(3,1,2);stem(0:length(y)-1,abs(y),.);grid on;ylabel(幅度);title(h(n)的频域幅频特性);subplot(3,1,3);stem(0:length(y)-1,angle(y),.);ylabel(相位);title(h(n)的频域相频特性);grid on;验证：程序：A=1;a=0.4;m=2.0734;N=50;fs=1;n=0:N-1;t=n/fs;x=A*exp(-a*t).*sin(m*t);y1=fft(x,N);f=n*fs/N;subplot(4,1,1);stem(f,y1,.);grid on;xlabel(f/Hz);title(xa(n)的频域特性);hb=1 2.5 2.5 1 zeros(1,47);y2=fft(hb,50);subplot(4,1,2);stem(0:length(y2)-1,y2,.);title(hb(n)的频域特性);grid on;y3=y1.*y2;subplot(4,1,3);stem(0:length(y3)-1,abs(y3),.);title(h(n)的幅频特性)ylabel(幅度);grid on;subplot(4,1,4);stem(0:length(y3)-1,angle(y3),.);title(h(n)的相频特性);ylabel(相位)grid on;结果：分析：与前面h(n)的的幅频和相频特性曲线相比较，发现频域相乘后的曲线与h(n)的曲线完全相同，由此验证了时域卷积定理。4. 思考题(1) 在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？它们所对应的模拟频率是否相同？为什么？答：采样率不同，采样序列的傅里叶变换频谱的数字频率度量不同。但它们所对应的模拟频率相同。(2) 在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如， 选M=10和M=20，分别做序列的傅里叶变换，求得所得结果之间有无差异？为什么？答：所得的傅氏变换序列长度会不一样，因为k的取值范围不一样。实验二 用FFT作谱分析 1. 实验目的及原理(1) 进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质)。(2) 熟悉FFT算法原理和FFT子程序的应用。(3) 学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。2. 实验步骤(1) 复习DFT的定义、性质和用DFT作谱分析的有关内容。(2) 复习FFT算法原理与编程思想，并对照DIT-FFT运算流图和程序框图， 读懂本实验提供的FFT子程序。(3) 编制信号产生子程序， 产生以下典型信号供谱分析用： (4) 编写主程序。(5) 按实验内容要求，上机实验，并写出实验报告。 3. 上机实验内容(1) 对 2 中所给出的信号逐个进行谱分析。下面给出针对各信号的FFT变换区间N以及对连续信号x6(t)的采样频率Fs，供实验时参考。 x1(n)，x2(n)，x3(n)，x4(n)，x5(n)：N=8，16x6(t)：Fs=64Hz，N=16，32，64(2) 令x(n)=x4(n)+x5(n)，用FFT计算8点和16点离散傅里叶变换。 X(k)=DFT x(n) (3) 令x(n)=x4(n)+jx5(n)， 重复(2)。实验源代码：I.（1）x1(n )= R4(n)的波形，N=8，N=16点的DFT幅频特性图和相频特性图。subplot(3,2,1);xn=ones(1,4); n=0:3;stem(n,xn,.);title(a)x(n)的波形);xlabel(n);ylabel(x(n);axis(0,4 0,1.5)Xk8=fft(xn,8);Xk16=fft(xn,16);subplot(3,2,3);k=0:7;stem(k,abs(Xk8),.);title(b)8点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,8,0,6)subplot(3,2,5);stem(k,angle(Xk8),.);title(d)8点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,8,-3.5,3.5)subplot(3,2,4);k=0:15;stem(k,abs(Xk16),.);title(c)16点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,16,0,6)subplot(3,2,6);stem(k,angle(Xk16),.);title(e)16点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,16,-3.5,3.5)结果为图2-1 图2-1（2）x2 (n )= 1 2 3 4 4 3 2 1 的波形，N=8，N=16点的DFT幅频特性图和相频特性图。subplot(3,2,1);xn=1 2 3 4 4 3 2 1 ;n=0:7;stem(n,xn,.);title(a)x(n)的波形);xlabel(n);ylabel(x(n);axis(0 8 0 5);Xk8=fft(xn,8);Xk16=fft(xn,16);subplot(3,2,3)；k=0:7;stem(k,abs(Xk8),.);title(b)8点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,8,0,21)subplot(3,2,5);stem(k,angle(Xk8),.); title(d)8点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,8,-3.5,3.5)subplot(3,2,4);k=0:15;stem(k,abs(Xk16),.);title(c)16点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,16,0,21)subplot(3,2,6);stem(k,angle(Xk16),.); title(e)16点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,16,-3.5,3.5)结果为图2-2图2-2（3）x3 (n )= 4 3 2 1 1 2 3 4 的波形，N=8，N=16点的DFT幅频特性图和相频特性图。subplot(3,2,1);xn=4 3 2 1 1 2 3 4 ;n=0:7;stem(n,xn,.);title(a)x(n)的波形);xlabel(n);ylabel(x(n);axis(0 8 0 5);Xk8=fft(xn,8);Xk16=fft(xn,16);subplot(3,2,3);k=0:7;stem(k,abs(Xk8),.);title(b)8点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,8,0,21)subplot(3,2,5);stem(k,angle(Xk8),.); title(d)8点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,8,-3.5,3.5)subplot(3,2,4);k=0:15;stem(k,abs(Xk16),.);title(c)16点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,16,0,21)subplot(3,2,6);stem(k,angle(Xk16),.); title(e)16点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,16,-3.5,3.5)结果为图2-3图2-3（4）x4 (n )=cos(n.pi/4)的波形，N=8，N=16点的DFT幅频特性图和相频特性图。subplot(3,2,1);n=0:10;xn=cos(pi*n/4);stem(n,xn,.);title(a)x(n)的波形);xlabel(n);ylabel(x(n);axis(0,11,-1.2,1.2);Xk8=fft(xn,8);Xk16=fft(xn,16);subplot(3,2,3);k=0:7;stem(k,abs(Xk8),.);title(b)8点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,8,0,8)subplot(3,2,5);stem(k,angle(Xk8),.); title(d)8点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,8,-3.5,3.5)subplot(3,2,4);k=0:15;stem(k,abs(Xk16),.);title(c)16点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,16,0,8)subplot(3,2,6);stem(k,angle(Xk16),.); title(e)16点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,16,-3.5,3.5)结果为图2-4图2-4(5)x5 (n )=sin(n.pi/8)的波形，N=8，N=16点的DFT幅频特性图和相频特性图。subplot(3,2,1);n=0:20;xn=sin(pi*n/8);stem(n,xn,.);title(a)x(n)的波形);xlabel(n);ylabel(x(n);axis(0,21,-1.2,1.2);Xk8=fft(xn,8);Xk16=fft(xn,16);subplot(3,2,3);k=0:7;stem(k,abs(Xk8),.);title(b)8点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,8,0,9)subplot(3,2,5);stem(k,angle(Xk8),.); title(d)8点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,8,-3.5,3.5)subplot(3,2,4);k=0:15;stem(k,abs(Xk16),.);title(c)16点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,16,0,9)subplot(3,2,6);stem(k,angle(Xk16),.); title(e)16点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,16,-3.5,3.5)结果为图2-5图2-5(6)x6 (n )=cos(8*pi*t)+cos(16*pi)+cos(20*pi*t)的波形，N=16，N=32,N=64点的DFT幅频特性图和相频特性图。(其中Fs=64Hz，即t=n/64)subplot(4,2,1);t=-0.8:0.01:0.8;xt=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);plot(t,xt);title(a)x(t)的波形);xlabel(t);ylabel(x(t); axis(-0.82 0.82 -3.2 3.2)subplot(4,2,2);n=0:20;xn=cos(n*pi/8)+cos(n*pi/4)+cos(5*n*pi/16);stem(n,xn,.);title(b)x(n)的波形);xlabel(n);ylabel(x(n);axis(0,21,-3.2,3.2);Xk16=fft(xn,16);Xk32=fft(xn,32);Xk64=fft(xn,64);subplot(4,2,3);k=0:15;stem(k,abs(Xk16),.);title(c)16点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,16,0,12)subplot(4,2,4);stem(k,angle(Xk16),.); title(d)16点DFT的相频特性图);xlabel(k);ylabel(相位);axis(0,16,-1.5,1.5)subplot(4,2,5);k=0:31;stem(k,abs(Xk32),.);title(e)32点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,32,0,16)subplot(4,2,6);stem(k,angle(Xk32),.); title(f)32点DFT的相频特性图);xlabel(k);ylabel(相位);axis(0,32,-2.5,2.5)subplot(4,2,7);k=0:63;stem(k,abs(Xk64),.);title(g)64点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,64,0,16)subplot(4,2,8);stem(k,angle(Xk64),.); title(h)64点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,64,-2.5,2.5)结果为图2-6图2-6II.令x(n)=x4(n)+x5(n)，用FFT计算8点和16点离散傅里叶变换。并画图。subplot(3,2,1);n=0:20;xn=cos(pi*n/4)+sin(pi*n/8);stem(n,xn,.);title(a)x(n)的波形);xlabel(n);ylabel(x(n);axis(0,21,-2.2,2.2);Xk8=fft(xn,8);Xk16=fft(xn,16);subplot(3,2,3);k=0:7;stem(k,abs(Xk8),.);title(b)8点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,8,0,6)subplot(3,2,5);stem(k,angle(Xk8),.); title(d)8点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,8,-3.5,3.5)subplot(3,2,4);k=0:15;stem(k,abs(Xk16),.);title(c)16点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,16,0,10)subplot(3,2,6);stem(k,angle(Xk16),.); title(e)16点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,16,-3.5,3.5)结果为图2-7图2-7III.令x(n)=x4(n)+jx5(n)，用FFT计算8点和16点离散傅里叶变换。并画图。subplot(3,2,1);n=0:20;xn=cos(pi*n/4)+j*sin(pi*n/8);stem(n,xn,.);title(a)x(n)的波形);xlabel(n);ylabel(x(n);axis(0,21,-2.2,2.2);Xk8=fft(xn,8);Xk16=fft(xn,16);subplot(3,2,3);k=0:7;stem(k,abs(Xk8),.);title(b)8点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,8,0,6)subplot(3,2,5);stem(k,angle(Xk8),.); title(d)8点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,8,-3.5,3.5)subplot(3,2,4);k=0:15;stem(k,abs(Xk16),.);title(c)16点DFT的幅频特性图);xlabel(k);ylabel(幅度);axis(0,16,0,10)subplot(3,2,6);stem(k,angle(Xk16),.); title(e)16点DFT的相频特性图)xlabel(k);ylabel(相位);axis(0,16,-3.5,3.5)以上程序执行结果为下图2-8图2-84. 思考题(1) 在N=8时，x2(n)和x3(n)的幅频特性会相同吗？为什么？N=16呢？N8时一样，N16时不一样。因为DFT变换可以看成是将该序列进行周期延拓后的傅里叶级数变换的主值序列。当N=8时，两序列进行周期延拓后序列相同，所以其傅里叶级数变换的主值序列也相同，进而DFT变换也相同。而当N=16时，两序列进行周期延拓后序列不相同，所以其傅里叶级数变换的主值序列也相同，进而DFT变换也不相同(2) 如果周期信号的周期预先不知道，如何用FFT进行谱分析？通过的不同取值，然后分析频谱，从而不断缩小的猜测范围，最后得出值。