数字信号处理实验七-九.doc

实验七 时域采样与频域采样一 实验目的 时域采样理论与频域采样理论是数字信号处理中的重要理论。要求掌握模拟信号采样前后频谱的变化，以及如何选择采样频率才能使采样后的信号不丢失信息；要求掌握频率域采样会引起时域周期化的概念，以及频率域采样定理及其对频域采样点数选择的指导作用。二 实验原理与方法 时域采样定理的要点是：a) 对模拟信号以间隔T进行时域等间隔理想采样，形成的采样信号的频谱是原模拟信号频谱以采样角频率（）为周期进行周期延拓。公式为： b) 采样频率必须大于等于模拟信号最高频率的两倍以上，才能使采样信号的 频谱不产生频谱混叠。 利用计算机计算上式并不方便，下面我们导出另外一个公式，以便用计算机上进行实验。 理想采样信号和模拟信号之间的关系为： 对上式进行傅立叶变换，得到： 在上式的积分号内只有当时，才有非零值，因此： 上式中，在数值上，再将代入，得到： 上式的右边就是序列的傅立叶变换，即 上式说明理想采样信号的傅立叶变换可用相应的采样序列的傅立叶变换得到，只要将自变量用代替即可。 频域采样定理的要点是：a) 对信号x(n)的频谱函数X(ej)在0，2上等间隔采样N点，得到则N点IDFT得到的序列就是原序列x(n)以N为周期进行周期延拓后的主值区序列，公式为： b) 由上式可知，频域采样点数N必须大于等于时域离散信号的长度M(即NM)，才能使时域不产生混叠，则N点IDFT得到的序列就是原序列x(n),即=x(n)。如果NM，比原序列尾部多N-M个零点；如果NM，z则=IDFT发生了时域混叠失真，而且的长度N也比x(n)的长度M短，因此。与x(n)不相同。 在数字信号处理的应用中，只要涉及时域或者频域采样，都必须服从这两个采样理论的要点。 对比上面叙述的时域采样原理和频域采样原理，得到一个有用的结论，这两个采样理论具有对偶性：“时域采样频谱周期延拓，频域采样时域信号周期延拓”。因此放在一起进行实验。三 实验内容及步骤（1）时域采样理论的验证。给定模拟信号， 式中A=444.128，=50，=50rad/s，它的幅频特性曲线如图10.2.1 图10.2.1 的幅频特性曲线 现用DFT(FFT)求该模拟信号的幅频特性，以验证时域采样理论。 安照的幅频特性曲线，选取三种采样频率，即=1kHz，300Hz，200Hz。观测时间选。 为使用DFT，首先用下面公式产生时域离散信号，对三种采样频率，采样序列按顺序用，表示。 因为采样频率不同，得到的，的长度不同， 长度（点数）用公式计算。选FFT的变换点数为M=64，序列长度不够64的尾部加零。X(k)=FFTx(n) ， k=0,1,2,3,-,M-1式中k代表的频率为 。要求： 编写实验程序，计算、和的幅度特性，并绘图显示。观察分析频谱混叠失真。（2）频域采样理论的验证。给定信号如下： 编写程序分别对频谱函数在区间上等间隔采样32和16点，得到： 再分别对进行32点和16点IFFT，得到： 分别画出、的幅度谱，并绘图显示x(n)、的波形，进行对比和分析，验证总结频域采样理论。提示：频域采样用以下方法容易变程序实现。 直接调用MATLAB函数fft计算就得到在的32点频率域采样 抽取的偶数点即可得到在的16点频率域采样，即。 当然也可以按照频域采样理论，先将信号x(n)以16为周期进行周期延拓，取其主值区（16点），再对其进行16点DFT(FFT),得到的就是在的16点频率域采样。四思考题： 如果序列x(n)的长度为M，希望得到其频谱在上的N点等间隔采样，当NM时， 如何用一次最少点数的DFT得到该频谱采样？五实验报告及要求a) 运行程序打印要求显示的图形，。 b) 分析比较实验结果，简述由实验得到的主要结论c) 简要回答思考题d) 附上程序清单和有关曲线。附： 实验程序清单1 时域采样理论的验证程序清单% 时域采样理论验证程序exp2a.mTp=64/1000;%观察时间Tp=64微秒%产生M长采样序列x(n)% Fs=1000;T=1/Fs;Fs=1000;T=1/Fs;M=Tp*Fs;n=0:M-1;A=444.128;alph=pi*50*20.5;omega=pi*50*20.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M); %M点FFTxnt)yn=xa(nT);subplot(3,2,1);tstem(xnt,yn);%调用自编绘图函数tstem绘制序列图box on;title(a) Fs=1000Hz);k=0:M-1;fk=k/Tp;subplot(3,2,2);plot(fk,abs(Xk);title(a) T*FTxa(nT),Fs=1000Hz);xlabel(f(Hz);ylabel(幅度);axis(0,Fs,0,1.2*max(abs(Xk)%=% Fs=300Hz和 Fs=200Hz的程序与上面Fs=1000Hz完全相同。注意：tstem函数是自己定义的一个画序列的函数，而并非matlab自带函数，所以需要将定义的M文件防到当前的目录下面，也就是把tstem.m文件拷贝到当前目录即可。2 频域采样理论的验证程序清单%频域采样理论验证程序exp2b.mM=27;N=32;n=0:M;%产生M长三角波序列x(n)xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=xa,xb;Xk=fft(xn,1024);%1024点FFTx(n), 用于近似序列x(n)的TFX32k=fft(xn,32);%32点FFTx(n)x32n=ifft(X32k);%32点IFFTX32(k)得到x32(n)X16k=X32k(1:2:N);%隔点抽取X32k得到X16(K)x16n=ifft(X16k,N/2);%16点IFFTX16(k)得到x16(n)subplot(3,2,2);stem(n,xn,.);box ontitle(b) 三角波序列x(n);xlabel(n);ylabel(x(n);axis(0,32,0,20)k=0:1023;wk=2*k/1024;%subplot(3,2,1);plot(wk,abs(Xk);title(a)FTx(n);xlabel(omega/pi);ylabel(|X(ejomega)|);axis(0,1,0,200)k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),.);box ontitle(c) 16点频域采样);xlabel(k);ylabel(|X_1_6(k)|);axis(0,8,0,200)n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,.);box ontitle(d) 16点IDFTX_1_6(k);xlabel(n);ylabel(x_1_6(n);axis(0,32,0,20)k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),.);box ontitle(e) 32点频域采样);xlabel(k);ylabel(|X_3_2(k)|);axis(0,16,0,200)n1=0:N-1;subplot(3,2,6);stem(n1,x32n,.);box ontitle(f) 32点IDFTX_3_2(k);xlabel(n);ylabel(x_3_2(n);axis(0,32,0,20)附： 实验程序运行结果1 时域采样理论的验证程序运行结果exp2a.m如图10.3.2所示。由图可见，采样序列的频谱的确是以采样频率为周期对模拟信号频谱的周期延拓。当采样频率为1000Hz时频谱混叠很小；当采样频率为300Hz时，在折叠频率150Hz附近频谱混叠很严重；当采样频率为200Hz时，在折叠频率110Hz附近频谱混叠更很严重。图10.2.22 时域采样理论的验证程序exp2b.m运行结果如图10.3.3所示。图10.3.3该图验证了频域采样理论和频域采样定理。对信号x(n)的频谱函数X(ej)在0，2上等间隔采样N=16时， N点IDFT得到的序列正是原序列x(n)以16为周期进行周期延拓后的主值区序列：由于NM，频域采样定理，所以不存在时域混叠失真，因此。与x(n)相同。简答思考题 先对原序列x(n)以N为周期进行周期延拓后取主值区序列，再计算N点DFT则得到N点频域采样：实验八 应用FFT对信号进行频谱分析1、 实验目的 1、加深对离散信号的DTFT和DFT的及其相互关系的理解。 2、在理论学习的基础上，通过本次实验，加深对快速傅里叶变换的理解，熟悉FFT算法极其程序的编写。 3、熟悉应用FFT对典型信号进行频谱分析的方法。 4、了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。2、 实验原理和方法 一个连续信号的频谱可以用它的傅里叶变换表示 （21） 如果对该信号进行理想采样，可以得到采样序列： （22） 同样可以对该序列进行Z变换，其中T为采样周期 （23） 当得时候，我们就得到了序列的傅里叶变换 （24） 其中称为数字频率，它和模拟频域的关系为 （25） 式中的是采样频率，上式说明数字频率是模拟频率对采样频率的归一化。同模拟域的情况相似，数字频率代表了序列值变化的速率，而序列的傅里叶变换称为序列的频谱。序列的傅里叶变换和对应的采样信号频谱具有下式的对应关系 ： （26） 即序列的频谱是采样信号频谱的周期延拓。从式（26）可以看出，只要分析采样序列的频谱，就可以得到相应的连续信号的频谱。注意：这里的信号必须是带限信号，采样也必须满足Nyquist定理。 在各种信号序列中 ，有限长序列在数字信号处理中占有很重要的地位。无限长的序列也往往可以用有限长序列来逼近。对于有限长的序列我们可以使用离散傅里叶变换（DFT），这一变换可以很好的放映序列的频域特性，并且容易利用快速算法在计算机上实现当序列的长度是N时，我们定义离散傅里叶变化为： （27）其中，它的反变换定义为： （28）根据式（23）和（27）令，则有 （29）可以得到，是Z平面单位圆上幅角为的点，就是见单位圆进行N等分以后第K个点。所以，是变换在单位圆上的等距采样，或者说是序列傅里叶变换的等距采样。时域采样在满足Nyquist定理时，就不会发生频谱混淆；同样地，在频率域进行采样的时候，只要采样间隔足够小，也不会发生时域序列混淆。 DFT时对序列傅里叶变换的等距采样，因此可以用于序列的频谱分析。在运用DFT进行频谱分析的时候可能有三种误差，分析如下： （1）混淆现象 从式（26）中可以看出，序列的频谱是采样信号频谱的周期延拓，周期是，因此当采样速率不满足Nyquist定理，即采样频率小于两倍的信号（这里指的是实信号）频率时，经过采样就会发生频谱混淆。这导致采样后的信号序列频谱不能真实的反映原信号的频谱。所以，在利用DFT分析连续信号频谱的时候，必须注意这一问题。这就告诉我们，在确定信号的采样频率之前，需要对频谱的性质有所了解。在一般的情况下，为了保证高于折叠频率的分量不会出现，在采样之前，先用低通模拟滤波器对信号进行滤波。（2）泄露现象 实际中的信号序列往往很长，甚至是无线长序列。为了方便，我们往往用截短的序列来近似它们，这样可以使用较短的DFT来对信号进行频谱分析。这种截短等价于给原始信号序列乘以一个矩形窗函数，而矩形窗函数的频谱不是有限带宽的，从而它和原信号的频谱进行卷积以后会扩展原信号的频谱。值得一提的是，泄漏时不能和混淆完全分离的，因为泄漏导致频谱的扩展，从而造成混淆。为了减小泄漏的影响，可以选择是党的窗函数使频谱的扩散减小到最小。（3） 栅栏效应 因为DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数。这样就产生了栅栏效应，从某种角度来看，用DFT来观看频谱就好像通过一个栅栏来观看一幅景象，只能在离散点上看到真实的频谱。这样的话就会有一些频谱的峰点或谷点被“栅栏”挡住，不能被我们观察到。减小栅栏效应的一个方法是在原序列的末端补一些零值，从而变动DFT的点数。这种方法的实质是人为地改变了对真是频谱采样的点数和位置，相当于搬动了“栅栏”的位置，从而使得原来被挡住的一些峰点或谷点显露出来。注意，这时候每根谱线对应的频率和原来的已经不相同了。 从上面的分析过程可以看出，DFT可以用于信号的频谱分析，但必须注意可能产生的误差，在应用过程中要尽可能减小和消除这些误差的影响。 快速傅里叶变换FFT并不是与DFT不相同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式（27）进行一次次的分解，使其成为若干小数点DFT的组合，从而减小运算量。常用的FFT是以2为基数，其长度。它的运算效率高，程序比较简单，使用也十分菲娜改变。当需要进行变换的序列的长度不是2的整数次方的时候，为了使用以2为基的FFT，可以用末尾补零的方法，使其长度延长至2的整数次方。IFFT一般可以通过FFT程序来完成，比较式（27）和（28），只要对取共轭，进行FFT运算，然后再取其共轭，并乘以因子，就可以完成IFFT。3、 实验内容及步骤 （一）编制实验用主程序及相应子程序1、在实验之前，认真复习DFT和FFT有关的知识，阅读本实验原理与方法和实验附录部分中和本实验有关的子程序，掌握子程序的原理并学习调用方法。2、编制信号产生子程序及本实验的频谱分析主程序。试验中需要用到的基本信号包括：（1）高斯序列：（2）衰减正弦序列（3）三角波序列，（4）反三角序列，（二）上机实验内容 1、观察高斯序列的时域和频域特性（1）固定信号的参数p=8,改变q的值，使q分别等于2，4，8。观察它们的时域和幅频特性，了解q取不同值的时候，对信号时域特性和幅频特性的影响。（2）固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列时域及幅频特性的影响。注意p等于多少时，会发生明显的泄漏现象，绘制相应的食欲序列和幅频特性曲线。1.高斯序列 n=0:15; p=8;q=8;x=exp(-1*(n-p).2/q); subplot(3,2,2);stem(abs(fft(x); subplot(3,2,1);stem(n,x); p=8;q=4;x=exp(-1*(n-p).2/q); subplot(3,2,3);stem(n,x);subplot(3,2,4);stem(abs(fft(x); p=8;q=2;x=exp(-1*(n-p).2/q); subplot(3,2,5);stem(n,x); subplot(3,2,6);stem(abs(fft(x); n=0:15;p=8;q=8;x=exp(-1*(n-p).2/q); subplot(3,2,2);stem(abs(fft(x); subplot(3,2,1);stem(n,x); p=13;q=8;x=exp(-1*(n-p).2/q); subplot(3,2,3);stem(n,x); subplot(3,2,4);stem(abs(fft(x); p=14;q=8;x=exp(-1*(n-p).2/q); subplot(3,2,5);stem(n,x); subplot(3,2,6);stem(abs(fft(x);2.衰减正弦序列 n=0:15; a=0.1;f=0.0625; x=exp(-a*n).*sin(2*pi*f*n); subplot(3,2,1);stem(n,x); subplot(3,2,2);stem(n,abs(fft(x); n=0:15; a=0.1;f=0.4375; x=exp(-a*n).*sin(2*pi*f*n); subplot(3,2,3);stem(n,x); subplot(3,2,4);stem(n,abs(fft(x); n=0:15; a=0.1;f=0.5625; x=exp(-a*n).*sin(2*pi*f*n); subplot(3,2,5);stem(n,x); subplot(3,2,6);stem(n,abs(fft(x);3三角波序列 for i=1:4 x(i)=i; end for i=5:8 x(i)=9-i; end subplot(2,1,1);stem(x); subplot(2,1,2);stem(abs(fft(x,16) 4.反三角序列 for i=1:4 x(i)=5-i; end for i=5:8 x(i)=i-4; end subplot(2,1,1);stem(x);subplot(2,1,2);stem(abs(fft(x,16)四、思考题 1.实验中的信号序列x（n）和x(n),在单位圆上的Z变换频谱和会相同吗？如果不同，你能说出哪一个低频分量更多一些吗？答：不同在单位圆上的Z变换频谱中，xc(n)的低频分量比xd(n)的多一些。 2.对一个有限长序列进行离散傅里叶变换（DFT）,等价于将该序列周期延拓后进行傅里叶级数展开。因为DFS也只是取其中一个周期来运算，所以FFT在一定条件下也可以用分析周期信号序列。如果正弦信号sin()，.1，用16点的FFT来做DFS运算，得到的频谱是信号本身的真实谱吗？答：当N与进行FFT变换的点数K一样的时候，可以认为DFS与FFT的变换时相等的，这时我们可以用DFS来分析FFT。但是在N与K不相等的时候，DFS与FFT变换不等价。如上面所说得正弦信号sin(2fn)，f=0.1，用16点的FFT来做DFT运算，N=10而K=16 。用16 点FFT得到到的频谱不是真实信号的频谱。五、实验总结： 1通过实验加深对离散信号的DTFT和DFT的及其相互关系的理解。 2在理论学习的基础上，通过本次实验，加深了对快速傅里叶变换的理解，熟悉了FFT算法极其程序的编写。 3了解了应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。实验九 系统响应及系统稳定性一. 实验目的（1）掌握求系统响应的方法。（2）掌握时域离散系统的时域特性。（3）分析、观察及检验系统的稳定性。二实验原理与方法在时域中，描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应，本实验仅在时域求解。在计算机上适合用递推法求差分方程的解，最简单的方法是采用MATLAB语言的工具箱函数filter函数。也可以用MATLAB语言的工具箱函数conv函数计算输入信号和系统的单位脉冲响应的线性卷积，求出系统的响应。系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。重点分析实验系统的稳定性，包括观察系统的暂态响应和稳定响应。系统的稳定性是指对任意有界的输入信号，系统都能得到有界的系统响应。或者系统的单位脉冲响应满足绝对可和的条件。系统的稳定性由其差分方程的系数决定。实际中检查系统是否稳定，不可能检查系统对所有有界的输入信号，输出是否都是有界输出，或者检查系统的单位脉冲响应满足绝对可和的条件。可行的方法是在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数（包括零），就可以断定系统是稳定的19。系统的稳态输出是指当时，系统的输出。如果系统稳定，信号加入系统后，系统输出的开始一段称为暂态效应，随n的加大，幅度趋于稳定，达到稳态输出。注意在以下实验中均假设系统的初始状态为零。三实验内容及步骤（1）编制程序，包括产生输入信号、单位脉冲响应序列的子程序，用filter函数或conv函数求解系统输出响应的主程序。程序中要有绘制信号波形的功能。（2）给定一个低通滤波器的差分方程为 输入信号 a) 分别求出系统对和的响应序列，并画出其波形。 b) 求出系统的单位冲响应，画出其波形。（3）给定系统的单位脉冲响应为 用线性卷积法分别求系统h1(n)和h2(n)对的输出响应，并画出波形。（4）给定一谐振器的差分方程为 令 ，谐振器的谐振频率为0.4rad。 a) 用实验方法检查系统是否稳定。输入信号为时，画出系统输出波形。 b) 给定输入信号为 求出系统的输出响应，并画出其波形。四思考题(1) 如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应? 如何求？ （2）如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号会有何变化，用前面 第一个实验结果进行分析说明。5实验报告要求（1）简述在时域求系统响应的方法。（2）简述通过实验判断系统稳定性的方法。分析上面第三个实验的稳定输出的波形。 （3）对各实验所得结果进行简单分析和解释。（4）简要回答思考题。（5）打印程序清单和要求的各信号波形。五实验程序思考题实验1：系统响应及系统稳定性close all;clear all%=内容1：调用filter解差分方程，由系统对u(n)的响应判断稳定性=A=1,-0.9;B=0.05,0.05; %系统差分方程系数向量B和Ax1n=1 1 1 1 1 1 1 1 zeros(1,50); %产生信号x1(n)=R8(n)x2n=ones(1,128); %产生信号x2(n)=u(n)hn=impz(B,A,58); %求系统单位脉冲响应h(n)subplot(2,2,1);y=h(n);tstem(hn,y); %调用函数tstem绘图title(a) 系统单位脉冲响应h(n);box ony1n=filter(B,A,x1n); %求系统对x1(n)的响应y1(n)subplot(2,2,2);y=y1(n);tstem(y1n,y);title(b) 系统对R8(n)的响应y1(n);box ony2n=filter(B,A,x2n); %求系统对x2(n)的响应y2(n)subplot(2,2,4);y=y2(n);tstem(y2n,y);title(c) 系统对u(n)的响应y2(n);box on%=内容2：调用conv函数计算卷积=x1n=1 1 1 1 1 1 1 1 ; %产生信号x1(n)=R8(n)h1n=ones(1,10) zeros(1,10);h2n=1 2.5 2.5 1 zeros(1,10);y21n=conv(h1n,x1n);y22n=conv(h2n,x1n);figure(2)subplot(2,2,1);y=h1(n);tstem(h1n,y); %调用函数tstem绘图title(d) 系统单位脉冲响应h1(n);box onsubp