中南大学数字信号处理实验报告.doc

中南大学数字信号处理实验报告-课程名称数字信号处理指导教师李宏学 院信息科学与工程学院专业班级姓 名学 号目录实验一 常见离散时间信号的产生和频谱分析3一、实验目的3二、实验原理3三、 实验内容6实验二 数字滤波器的设计12一、实验目的12二、实验原理12三、实验内容16实验一 常见离散时间信号的产生和频谱分析一、实验目的（1）熟悉MATLAB应用环境，常用窗口的功能和使用方法；（2）加深对常用离散时间信号的理解；（3）掌握简单的绘图命令；（4）掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对离散信号进行频域分析。二、实验原理（1）常用离散时间信号a）单位抽样序列 如果在时间轴上延迟了k个单位，得到即： b）单位阶跃序列 c）矩形序列 d）正弦序列e）实指数序列f）复指数序列（2）离散傅里叶变换：设连续正弦信号为这一信号的频率为，角频率为，信号的周期为。如果对此连续周期信号进行抽样，其抽样时间间隔为T，抽样后信号以表示，则有，如果令为数字频率，满足，其中是抽样重复频率，简称抽样频率。为了在数字计算机上观察分析各种序列的频域特性，通常对在上进行M点采样来观察分析。 对长度为N的有限长序列x(n), 有其中 通常应取得大一些，以便观察谱的细节变化。取模可绘出幅频特性曲线。（3）用DFT进行普分析的三种误差三种误差：混叠现象、泄露现象、栅栏效应 a) 混叠现象 当采样频率小于两倍信号（这里指是信号）最大频率时，经过采样就会发生频谱混叠，这使得采样后的信号序列频谱不能真实地反映原信号的频谱。所以在利用DFT分析连续信号的频谱时，必须注意这一问题。避免混叠现象的唯一方法是保证采样速率足够高，使频谱交叠现象不致出现。也就是说，在确定采样频率之前，必须对信号的性质有所了解，一般在采样前，信号通过一个防混叠低通滤波器。b) 泄漏现象 实际中的信号序列往往很长,为了方便我们往往用截短的序列来近似它们,这样可以使用较短的DFT来对信号进行频谱分析,这种截短等价于给原信号序列乘以一个矩形窗函数。泄漏是不能与混叠完全分离开的,因为泄漏导致频谱的扩散，从而造成混叠。为了减小泄漏的影响，可以选择适当的窗函数，使频谱的扩散减到最小。c) 栅栏效应 因为DFT是对单位圆上Z变换的均匀采样，所以他不可能将频谱视为一个连续函数。这样就产生了栅栏效应，就一定意义上看，DFT来观看频谱就好像通过一个尖桩的栅栏来观看一个图景一样，只能在离散点上看到真实频谱，这样就可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所挡住，不能被我们观察到。减小栅栏效应的一个方法就是借助在原序列的末端添补一些零值，从而变动DFT的点数。这一方法实际上是人为地改变了对真实谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得频谱的峰点或者谷点暴露出来。当然，这是每根谱线所对应的频率和原来的不同了。 综上所述，DFT可以用于信号的频谱分析，但必须注意可能产生的误差，在应用过程中要尽可能减少和消除这些误差的影响。 3、 实验内容（1）复习常用离散时间信号的有关内容在时间上依次出现的数值序列,例如，0.5，1，2，-1，0，5，。相邻两个数之间的时间间隔可以是相等的，也可以是不等的。在前一情况下，设时间间隔为T秒,则离散信号可用符号x(nT)来表示。在间隔T归一化为1的条件下,T可以省略,即将x(nT)表示为x(n)。x(n)既可表示整个序列, 也可表示离散信号在nT瞬间的值。（2）用MATLAB编程产生上述任意3种序列（长度可输入确定，对(d) (e) (f)中的参数可自行选择），并绘出其图形；subplot(2,2,1)n=-20:20;xn=heaviside(n);xn(n=0)=1;plot(n,xn);stem(n,xn);axis(-20 20 0 1.2);title(单位阶跃序列);xlable(n);ylable(u(n);subplot(2,2,2);n=-20:20;N=5;xn=heaviside(n)-heaviside(n-N);xn(n=0)=1;xn(n=N)=0;plot(n,xn);stem(n,xn);axis(-20 20 0 1.2);title(矩阵序列);xlabel(n);ylabel(R_N(n);box onsubplot(2,2,3);n=-40:40;A=2;w=pi/8;f=pi/4;xn=A*sin(w.*n+f);plot(n,xn);stem(n,xn);axis(-40 40 -4.2 4.2)title(正弦序列);xlabel(n);ylabel(x(n);box on（3） 混叠现象 对连续信号其中，进行采样，分别取采样频率，观察的变化，并做记录(打印曲线),观察随着采样频率降低频谱混叠是否明显存在，说明原因。subplot(2,2,1);f01=500;fs=2500;n=1/fs;N=length(t);t=0:n:0.1;x=sin(2*pi*f01*t);plot(t,x);X=fft(x,N);plot(fs/N*(0:(N/2-1),abs(X(1:(N/2);title(取样频率为2500HZ);subplot(2,2,2);f01=500;fs=2000;n=1/fs;N=length(t);t=0:n:0.1;x=sin(2*pi*f01*t);plot(t,x);X=fft(x,N);plot(fs/N*(0:(N/2-1),abs(X(1:(N/2);title(取样频率为2000HZ);subplot(2,2,3);f01=500;fs=1200;n=1/fs;N=length(t);t=0:n:0.1;x=sin(2*pi*f01*t);plot(t,x);X=fft(x,N);plot(fs/N*(0:(N/2-1),abs(X(1:(N/2);title(取样频率为1200HZ);subplot(2,2,4);f01=500;fs=800;n=1/fs;N=length(t);t=0:n:0.1;x=sin(2*pi*f01*t);plot(t,x);X=fft(x,N);plot(fs/N*(0:(N/2-1),abs(X(1:(N/2);title(取样频率为800HZ); （4）截断效应给定，截取一定长度的信号，为窗函数，长度为N，。分别取N=6，8，12，计算的N点DFT变换，画出其幅频特性曲线；做2N点DFT变换，分析当N逐渐增大时，分析是否有频谱泄露现象、主瓣的宽度变化？如何减小泄露？subplot(2,2,1);n=50;Rn=ones(1,n);Wn=Rn;xn=cos(pi./4)*(0:n-1);yn=xn.*Wn;N=20;Y=fft(yn,2*N);plot(2*pi/N*(0:N/2-1),abs(Y(1:N/2);title(截断效应N=20);subplot(2,2,2);n=50;Rn=ones(1,n);Wn=Rn;xn=cos(pi/4)*(0:n-1);yn=xn.*Wn;N=12;Y=fft(yn,2*N);plot(2*pi/N*(0:N/2-1),abs(Y(1:N/2);title(截断效应N=12);subplot(2,2,3);n=50;Rn=ones(1,n);Wn=Rn;xn=cos(pi./4)*(0:n-1);yn=xn.*Wn;N=8;Y=fft(yn,2*N);plot(2*pi/N*(0:N/2-1),abs(Y(1:N/2);title(截断效应N=8);subplot(2,2,4);n=50;Rn=ones(1,n);Wn=Rn;xn=cos(pi./4)*(0:n-1);yn=xn.*Wn;N=6;Y=fft(yn,2*N);plot(2*pi/N*(0:N/2-1),abs(Y(1:N/2);title(截断效应N=6);（5）栅栏效应给定，分别计算在频率区间上的16点、32点、64点等间隔采样，绘制采样的幅频特性图，分析栅栏效应，如何减小栅栏效应？n=0:1:10;xn=ones(1,4),zeros(1,7);xk16=fft(xn,16);xk32=fft(xn,32);xk64=fft(xn,64);subplot(2,2,1);stem(n,xn,.);title(x(n);xlable(n);ylable(x(n);k=0:15;wk=2*k/16;subplot(2,2,2);stem(wk,abs(xk16),.);title(16点DFT的幅频特性图);xlable(omega/pi);ylable(幅度);k=0:31;wk=2*k/32;subplot(2,2,3);stem(wk,abs(xk32),.);title(32点DFT的幅频特性图);xlable(omega/pi);ylable(幅度);k=0:63;wk=2*k/64;subplot(2,2,4);stem(wk,abs(xk64),.);title(64点DFT的幅频特性图);xlable(omega/pi);ylable(幅度);实验二 数字滤波器的设计一、实验目的（1）熟悉用双线性变换法和脉冲响应不变法设计IIR数字滤波器的原理与方 法；（2）学会调用MATLAB信号处理工具中滤波器设计函数，设计各种IIR滤波 器，学会根据滤波需求确定滤波器指标参数；（3）掌握用窗函数法设计FIR数字滤波器的原理和方法；二、实验原理（一）IIR滤波器 模拟滤波器设计巴特沃兹滤波器的振幅平方函数为 （1）其传输函数为 （2） （3）首先确定技术指标：i. 通带中允许的最大衰减和通带截止频率；ii. 阻带允许的最小衰减和阻带起始频率。由式（8-11）可得： （4） （5）得到 （6） （7）再利用上面两式得到令， 则 （8）已知，可由式(8)求出滤波器的阶数N。求出的N可能有小数部分一般取大于等于N的最小整数。关于3dB截止频率，有时在技术指标中给出，如果没有给出可以按照式（6）或式（7）求出。根据以上所述，巴特沃兹滤波器的设计步骤为：i. 根据要求，由式(8)求出阶数N；ii. 由式（6）或式(7)求出3dB截止频率；iii. 由式（3）求出N个极点；iv. 由式（2）写出传递函数。双线性变换法和脉冲响应不变法实现将数字频率和模拟频率相互映射，具体见教材相关介绍。（二）FIR滤波器 设所希望得到的滤波器的理想频率响应为。那么FIR滤波器的设计就在于寻找一个传递函数去逼近。在这种逼近中最直接的一种方法是从单位取样响应序列着手，使逼近理想的单位取样响应。我们知道可以从理想频率响应通过傅里叶反变换来得到，即： （9）但是一般来说，这样得到的单位取样响应往往都是无限长序列；而且是非因果的。以一个截止频率为的线性相应位理想低通为例来说明。设低通滤波器的时延为，即： （10）则 这是一个以为中心的偶对称的无限长非因果序列。这样一个无限长的序列怎样用一个有限长序列去近似呢？最简单的办法就是直接截取它的一段来代替它。例如把到的一段截取来作为，但是为要保证所得到的是线性相位滤波器。必须满足的对称性，所以时延应该取长度的一半，即这种直接截取的办法可以形象地想象为，好比是通过一个“窗口”所看到的一段。中表达为和一个“窗口函数”的乘积。在这里，窗口函数就是矩形脉冲函数，即 但是一般来说，窗口函数并不一定是矩形函数，可以在矩形以内还对作一定的加权处理，因此，一般可以表示为 这里就是窗口函数。这种对理想单位取样响应加窗的处理对频率响应会产生以下三点影响：a) 使理想特性不连续的边沿加宽，形成一过渡带，过渡带的宽度取决于窗口频谱的主瓣宽度。b) 在过渡带两旁产生肩峰和余振，它们取决于窗口频谱的旁瓣；旁瓣越多，余振也越多；旁瓣相对值越大，肩峰则越强。c) 增加截取长度N，只能缩小窗口频谱的主瓣宽度而不能改变旁瓣的相对值；旁瓣与主瓣的相对关系只决定于窗口函灵敏的形状。因此增加N，只能相对应减小过渡带宽。而不能改变肩峰值。肩峰值的大小直接决定通带内的平稳和阻带的衰减，对滤波器性能有很大关系。例如矩形窗的情况下，肩峰达8.95%，致使阻带最小衰减只有21分贝，这在工程上往往是不够的。怎样才能改善阻带的衰减特性呢？只能从改善窗口函数的形状上找出路，所以希望的窗口频谱中应该减少旁瓣，使能量集中在主瓣，这样可以减少肩峰和余振，提高阻带的衰减。而且要求主瓣宽度尽量窄，以获得较陡的过渡带，然而这两个要求总不能兼得，往往需要用增加主瓣宽度带换取决瓣的抑制，于是提出了海明窗、凯塞-贝塞尔窗、切比雪夫窗等窗口函数。三、实验内容(1) 设计IIR数字滤波器，绘制幅频响应特性曲线。基于脉冲响应和双线性变换设计低通数字滤波器：要求其通带截至频率100Hz，阻带截至频率200Hz，通带衰减小于2dB，阻带衰减大于15dB，采样频率Fs=500HZ。rp=2;rs=15;Fs=500;Ts=1/Fs; wp=100*2*pi*Ts; %利用公式，求数字域通带截止频率wpws=200*2*pi*Ts; %利用公式，求数字域通带截止频率ws%T=1s时的模拟滤波器指标N,wc=buttord(wp,ws,rp,rs,s);%计算响应的模拟滤波器阶数N和截止频率wcB,A=butter(N,wc,s);%计算响应的模拟滤波器系统函数Bz,Az=impinvar(B,A);%用脉冲响应不变法将模拟滤波器转化为数字滤波器H,W=freqz(Bz,Az);plot(W/pi,abs(H),k);grid on;xlabel(频率w/pi);ylabel(频率相应幅度/dB);title(低通滤波器的幅频相应曲线(脉冲响应不变法);legend(Rp=2dB,Rs=15dB,Fs=500HZ);Rp=2;Rs=15;Fs=500;Ts=1/Fs; %双线性变换时频率的预畸变wp=100*2*pi*Ts; %利用公式，求数字域通带截止频率wpws=200*2*pi*Ts; %利用公式，求数字域通带截止频率wswp1=2*tan(wp/2)/Ts; %利用公式，进行预畸变ws1=2*tan( ws/2)/Ts;N,Wn=buttord(wp1,ws1,Rp,Rs,s); %选择滤波器的最小阶数NZ,P,K=buttap(N); %计算Butterworth低通原型模拟滤波器Bap,Aap=zp2tf(Z,P,K); %零极点增益到传递函数模型的转变b,a=lp2lp(Bap,Aap,Wn); %把模拟滤波器原型转换成截止频率为Wn的低通滤波器bz,az=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换H,W=freqz(bz,az);plot(W/pi,abs(H),k);grid on;xlabel(频率w/pi);ylabel(频率相应幅度/dB);title(低通滤波器的幅频相应曲线（双线性法）);legend(Rp=2dB,Rs=15dB,Fs=500HZ);（2）设计FIR数字滤波器，绘制幅频响应特性曲线。(a) 用矩形窗函数法设计一个线性