数字信号课程设计

1、一、利用DFT分析离散序列频谱应用傅里叶变换DFT，分析各种离散序列x(n)的频谱。离散周期信号的分析 离散周期信号可以展开成离散傅里叶级数，其中傅里叶系数如下所示：式中：N是序列的周期，n为时间离散变量， k为数字频率离散变量，是k次谐波的数字频率。所以离散周期信号的频谱X (k)是一个以N为周期的周期性离散频谱，各谱线之间的间隔为0=2/ N，而且存在着谐波的关系。 离散周期序列在时域与频域都是离散的、以N为周期的序列。正确地在一个周期内选取N，即可以准确地求取周期序列的频谱。离散周期序列频谱的求解步骤： （1）确定离散周期序列的基本周期N； （2）使用fft命令作N点FFT计算。频率分辨

2、率F0=2/N（3）令X(k0 )=1/NX（k）。 1利用DFT计算序列的频谱（1）程序N=100;n=0:N-1;x0=(0.5).n;x1=(n>=0);x=x0.*x1;X=1/N*fft(x,N);omega=2*pi/N*(n-N/2);subplot(2,1,1); stem(omega,abs(fftshift(X); ylabel('Magnitude'); xlabel('Frequency (rad)' );subplot(2,1,2); stem(omega,angle(fftshift(X); ylabel('Phase&

3、#39;); xlabel('Frequency (rad) ' );（2）程序N=256;n=0:N-1; x0=(0.5).n;x1=(n>=0);x=x0.*x1;X=1/N*fft(x,N);omega=2*pi/N*(n-N/2);subplot(2,1,1); stem(omega,abs(fftshift(X); ylabel('Magnitude'); xlabel('Frequency (rad)' );subplot(2,1,2); stem(omega,angle(fftshift(X); ylabel('Ph

4、ase'); xlabel('Frequency (rad) ' );（3）程序N=300;n=0:N-1;x0=(0.5).n;x1=(n>=0);x=x0.*x1;X=1/N*fft(x,N);omega=2*pi/N*(n-N/2);subplot(2,1,1); stem(omega,abs(fftshift(X); ylabel('Magnitude'); xlabel('Frequency (rad)' );subplot(2,1,2); stem(omega,angle(fftshift(X); ylabel('

5、;Phase'); xlabel('Frequency (rad) ' );二、连续系统的分析零输入、零状态法求解系统的响应MATLAB的控制工具箱（control toolbox）里包含了许多可用于分析线性时不变（LTI）系统的函数，使用命令help control可以查看控制工具箱里的这些函数。 对于线性时不变的连续系统，在时域中其数学模型用常系数线性微分方程来描述。如下：在拉普拉斯变换中，可用系统函数H(s)来描述，且根据其零，极点的分布情况，可以决定系统的结构及系统的稳定性。 另外还可以用状态变量分析法进行系统描述。此设计完成用atlab实现几种分析方法的相互转

6、换。在Matlab中，描述系统的传递函数（系统函数）型tf(transfer function),零极点型zp(zero pole)以及状态变量型ss(state space)三种方式可以方便的转换。Matlab相应的语句为：tf2zp系统函数型转换到零极点型tf2ss系统函数型转换到状态变量型zp2tf零极点型转换到系统函数型zp2ss零极点型转换到状态变量型ss2tf状态变量型转换到系统函数型ss2zp状态变量型转换到零极点型1:已知一个因果LTI离散系统的系统函数为： （1）阶跃响应M=30; N=32;omega=0:(N-1)*2*pi/N;H=(0.0675+0.1349*exp(

7、-j.*omega)+0.0675*exp(-2*j.*omega)./(1-1.143*exp(-j.*omega)+0.4128*exp(-2*j.*omega);x=ones(1,M);X=fft(x,N);Y=X.*H ;y=ifft(Y,N);figure(1);subplot(2,1,1); stem(0:M-1,x); xlabel('Time index n'); ylabel('xn');title ('Input signal and output signal');subplot(2,1,2); stem(0:M-1,rea

8、l(y(1:M); xlabel('Time index n'); ylabel('yn');（2）单位脉冲响应h(n)M=30; N=32;omega=0:(N-1)*2*pi/N;H=(0.0675+0.1349*exp(-j.*omega)+0.0675*exp(-2*j.*omega)./(1-1.143*exp(-j.*omega)+0.4128*exp(-2*j.*omega);x=1,zeros(1,M-1);X=fft(x,N);Y=X.*H ;y=ifft(Y,N);figure(1);subplot(2,1,1); stem(0:M-1,x)

9、; xlabel('Time index n'); ylabel('xn');title ('Input signal and output signal');subplot(2,1,2); stem(0:M-1,real(y(1:M); xlabel('Time index n'); ylabel('yn');三、 IIR数字低通滤波器的设计（双线性变换法）从模拟滤波器设计IIR数字滤波器就是要由模拟滤波器的系统函数进一步求得。归根结底是一个由s平面到z平面的变换。这个变换应遵循两个基本的目标。（1）H(z)的频

10、响必须要模仿的频响，也即s平面的虚轴j应该映射到z平面的单位圆上。（2）的因果稳定性，通过映射后仍应在所得到的H(z)中保持，也即s平面的左半平面应该映射到z平面单位圆以内。从模拟滤波器映射(变换)成数字滤波器常用的是脉冲响应不变变换法和双线性变换法两种。其中脉冲响应不变法设计数字滤波器的优点是器频率坐标的转换是线性的，在不存在频率混叠的情况下，能够完全逼近模拟滤波器的频率特性。但是当模拟滤波器的传输函数超出时，映射后必然存在频谱混叠的现象，这是脉冲响应不变法不足的地方，所以一般只用来设计低通、带通滤波器。对于高通、带通滤波器，由于它们在高频部分不衰减，当一定要追求频率响应线性而采用脉冲响应不

11、变法时，必须先对模拟高通、带阻滤波器加一保护滤波器，从而加大了设计成不，一般不采用这种方法。脉冲响应不变法设计滤波器时，从模拟预到数字域的转换是“多对一”的关系，而双线性变换法是“一对一”的关系，从而避免了频率混叠现象，因此本课程设计中采用双线性变换法设计滤波器，下面主要介绍双线性变换法的原理。1. 变换原理双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了克服脉冲响应不变法的多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的平面的一横带里（宽度为，即从到），然后再通过上面讨论过的标准变换关系将此横带变换到整个z平面上去，这样就使s平面与z平面是一一对应

12、的关系，消除了多值变换性，也就消除了频谱混叠现象，基本原理如图2-3所示。图2-3 双线性变换法的映射关系将s平面整个平面压缩到平面的到，可采用以下的变换关系 （2-3）其中C为常数；这样经（2-3）式变为，变为，可将上式写成令，则可得 （2-4）再将平面通过以下标准变化关系映射到z平面 （2-5）这样（2-4）式可表示为 （2-6） （2-7）2变换常数C的选择为了使模拟滤波器与数字滤波器在低频处有较确切的对应关系，即在低频处有，当较小时有由（2-3）式可知因而得 （2-8）则（2-6）和（2-7）式可重新写成： （2-9）即 （2-10）3.逼近情况双线性变换具备模拟域到数字域映射变换的总

13、要求，现分析如下：（1）将代入到(2-9)式则得或 （2-11）由上式可见，当时，；当时，；当时，。这就是说双线性变换把s左半平面映射在单位圆的内部；把s平面的整个轴映射成单位圆，把s右半平面映射在单位圆的外部。(2)令，则由(2-9)式得所以由此得出模拟滤波器的频率和数字滤波器频率的关系式为 (2-12)这一公式的关系如图2-4所示。可以看出，当时，当时，,当时，。这就是说：s平面的原点映射为z平面(1，0)点，而s平面的正虚轴和负虚轴分别映射成z平面单位圆的上半圆和下半圆。图2-4 双线性变换的频率间非线性关系由上所述，可得如下结论：a模拟滤波器中最大和最小值将保留在数字滤波器中，因此模拟

14、滤波器的通带或阻带变换成数字滤波器的通带或阻带。b如果模拟滤波器是稳定的，则通过双线性变换后所得的数字滤波器也一定是稳定的。c由于s平面的整个虚轴映射为z平面上的单位圆，因此双线性变换法确实消除了脉冲响应不变变换法所存在的混叠误差，所以逼近是良好的。但由(2-13)式可见，在频率与间存在严重的非线性。4. 模拟滤波器的数字化由于双线性变换法中，s与z之间有简单的代数关系，故可由模拟系统函数通过代数置换直接得到数字滤波器的系统函数。即 (2-14)可见数字滤波器的极点数等于模拟滤波器的极点数。频率响应也可用直接置换得到 (2-15)这一公式可用于将滤波器的数字域指标，转换为模拟域指标。再者，可在

15、未进行双线性变换前把原模拟系统函数分解成并联或级联子系统函数，然后再对每个子系统函数分别加以双线性变换。就是说，所有的分解，都可以就模拟滤波器系统函数来进行，因为模拟滤波器已有大量图表可供利用，且分解模拟系统函数比较容易。3.2 IIR滤波器的设计由前面IIR滤波器的设计原理确定运用双线性变换法设计IIR滤波器的流程如图3-1所示：数字滤波器指标模拟滤波器指标设计模拟滤波器对应的数字滤波器器双线性变换开始结束图3-1 IIR滤波器设计流程图3.3 IIR低通滤波器的设计1.IIR低通滤波器设计指标：，阻带最小衰减，通带最大衰减。2.IIR低通滤波器的设计在Matlab中设计巴特沃斯型、切比雪夫

16、I型和椭圆型IIR低通滤波器的程序分别如下：%巴特沃斯滤波器设计fp=1000;fs=1200;Rp=1;Rs=10;Fs=22050;wp=2*pi*fp/Fs;%确定数字滤波器指标ws=2*pi*fs/Fs;Wp=2*tan(wp/2);%转换为模拟滤波器指标Ws=2*tan(ws/2);N,Wn=buttord(Wp,Ws,Rp,Rs,'s');%选择巴特沃斯滤波器的最小阶数N和截止频率Wnb,a=butter(N,Wn,'s');%设计巴特沃斯滤波器b1,a1=bilinear(b,a,1);%运用双线性变换为数字滤波器%切比雪夫I型滤波器设计fp=10

17、00;fs=1200;Rp=1;Rs=10;Fs=22050;wp=2*pi*fp/Fs;%确定数字滤波器指标ws=2*pi*fs/Fs;Wp=2*tan(wp/2);%转换为模拟滤波器指标Ws=2*tan(ws/2);N,Wn=cheb1ord(Wp,Ws,Rp,Rs,'s');%选择切比雪夫I型滤波器的最小阶数N和截止频率Wnb,a=cheby1(N,Rp,Wn,'s');%设计切比雪夫I型滤波器b1,a1=bilinear(b,a,1);%运用双线性变换法变为数字滤波器%椭圆型滤波器设计fp=1000;fs=1200;Rp=1;Rs=10;Fs=22050

18、;wp=2*pi*fp/Fs;%确定数字滤波器指标ws=2*pi*fs/Fs;Wp=2*tan(wp/2);%转换为模拟滤波器指标Ws=2*tan(ws/2);N,Wn=ellipord(Wp,Ws,Rp,Rs,'s');%选择椭圆型滤波器的最小阶数N和截止频率Wnb,a = ellip(N,Rp,Rs,Wn,'s');%设计椭圆型滤波器b1,a1=bilinear(b,a,1);%运用双线性变换法变为数字滤波器得到巴特沃斯滤波器、切比雪夫I型滤波器和椭圆型滤波器的幅频特性分别如图3-2、图3-3和图3-4所示。图3-2 IIR低通巴特沃斯型滤波器幅频特性图3-

19、3 IIR低通切比雪夫I型滤波器幅频特性图3-4 IIR低通椭圆型滤波器幅频特性总结经过为期一周的课程设计实践，我完全自主地完成了既定的任务，实现了课程设计的要求，这使我对DSP课程有了更深一步的认识和了解，要想学好它重在实践，要通过不断的实际操作才能更好地学习和运用知识，并且基本掌握了运用MATLAB软件进行滤波器的设计和GUI界面的设计。通过课程设计这一实践环节，加深了我对本学期所学知识的理解，也提高了我运用所学知识分析问题和解决实际问题的能力。在此过程中我也发现自己好多不足之处，首先是自己对所学的基础知识掌握不深，还有对理论知识的运用还有欠缺，通过这次课程设计实践，使我在这几个方面都有所提高。由于时间原因，这次课程设计最大的遗憾是没能实现对语音信号进行实时录入，然后进行采样分析，接下来我会努力尝试实现这项功能。在这个过程中，我也曾经因为实践经验的缺乏失落过，也曾经因设计成功而热情高涨。生活就是这样，汗水预示着结果也见证着收