基于 MATLAB 的数字信号处理的模拟与仿真.doc

《基于 MATLAB 的数字信号处理的模拟与仿真.doc》由会员分享，可在线阅读，更多相关《基于 MATLAB 的数字信号处理的模拟与仿真.doc（16页珍藏版）》请在三一办公上搜索。

1、基于 MATLAB 的数字信号处理的模拟与仿真王靖斌山东农业大学信息科学与工程学院，山东泰安（371018）摘要：本文介绍了数字信号处理中各种算法基本原理，包括 Z 变换、DFT 变换及其快速 算法 FFT 等信号变换算法和包括 IIR 与 FIR 数字滤波器的设计与实现，详细论述了利用- 16 -MATLAB软件对信号变换和设计数字滤波器的原理 步骤和实现方法，并给出了基于MATLAB 的模拟与仿真。文章最后给出了实验结果，并就其结果做了进一步的解释和说明。 关键词：Z 变换，离散傅立叶变换 DFT，数字滤波器，系统结构仿真1引言数字信号处理（Digital Signal Processin

2、g,DSP）是一门交叉性的学科，它的理论基础涉 及信息、通信、雷达、航空航天及生物医药等众多学科领域，其成果又为这些学科的发展起 着重要的促进作用。MATLAB 是数字信号处理技术实现的重要手段。MATLAB 提供了信号处理工具箱， MATLAB 信号处理工具箱随着信号处理理论与方法的发展而发展，同时又为信号处理理论 与方法的实现与工程应用提供了有效的帮助。本文基于 MATLAB 实现了数字信号处理的信 号变换模拟与滤波器设计仿真。2信号变换的各种算法2.1 Z 变换的算法原理及 MATLAB 实现若序列为 x(n)，则 z 变换为X ( z) = （x n）z nn=这种变换将离散系统的差分

3、方程转化为简单的代数方程。以一个简单的差分方程为例，用 MATLAB 实现 z 变换： 系统的差分方程为：y(n)-2y(n-1)+3y(n-2)=4u(n)-5u(n-1)+6u(n-2)-7u(n-3)其初始条件为 x(-1)=1,x(-2)=-1,y(-1)=-1,y(-2)=1，求系统的输出 y(n).MATLAB 程序如下所示：clear all; close all; clc;b=4,-5,6,-7; a=1 -2 3; x0=1 -1; y0=-1 1;xic=filtic(b,a,y0,x0)bxplus=1; axplus=1 -1; ayplus=conv(a,axplus

4、)byplus=conv(b,bxplus)+conv(xic,axplus)R,P,K=residuez(byplus,ayplus) Mp=abs(P)Ap=angle(P)*180/piN=100; n=0:N-1; xn=ones(1,N); yn=filter(b,a,xn,xic); plot(n,yn)运行结果如下：xic =-1616-7ayplus =1-35-3byplus =-1227-170R =-5.5000 - 1.0607i-5.5000 + 1.0607i-1.0000P =1.0000 + 1.4142i1.0000 - 1.4142i1.0000K =0Mp

5、 =1.73211.73211.0000Ap =54.7356-54.73560输出相应曲线如下图所示图 1.差分方程的输出结果2.2 离散傅立叶变换（DFT）2.2.1 DFT 的算法原理及 MATLAB 实现 对 N 点有限长序列 x(n),其正变换为N 1 j 2 nkN 1NDFT x(n) = X (k ) =其逆变换为n=0x(n)eN=n= 0x(n)W nk , k = 0,1N 11 N 1j 2 nk1 N 1NIDFT X (k ) = x(n) =N n =0X (k )e N=N n=0X (k )W nk , n = 0,1N 1下面用 MATLAB 实现傅立叶变换

6、及其逆变换：离散傅立叶变换的 MATLAB 实现functionXk=dft(xn,N)Xk=在 0=k=N-1 间的 DFT 系数数组xn=N 点有限长度序列N=DFT 的长度n=0:1:N-1;%n 的行向量 k=0:1:N-1;%k 的行向量 WN=exp(-j*2*pi/N);%Wn 因子nk=n*k;%产生一个含 nk 值的 N 乘 N 维矩阵WNnk=WN.nk;%DFT 矩阵Xk=xn*WNnk;%DFT 系数的行向量 离散傅立叶逆变换的 MATLAB 实现：Functionxn=idft(Xk,N)xn=在 0=k=N-1 间的 N 点有限长度序列 Xk=在 0=k=N-1 间

7、的 DFT 系数数组 N=DFT 的长度n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k; WNnk=WN.（-nk）;xn=(Xk*WNnk)/N;以一个有限序列为例，计算它的傅立叶变换：若 x(n)=sin(n/4)是一个 N=32 的有限序列，实现它的傅立叶变换的 MATLAB 程序如下：clear all; close all; clc; N=32; n=0:N-1;xn=cos(pi*n/6);k=0:N-1;WN=exp(-j*2*pi/N); nk=n*k; WNnk=WN.nk; Xk=xn*WNnk; figure(1);stem(n,

8、xn) figure(2); stem(k,abs(Xk);得到的结果若下图图 2.有限长序列图 3.DFT 运行结果2.2.2 直接计算 DFT 存在的问题及改进方法 有限长序列的重要特点是其频域也可离散化成有限长序列，即可进行离散傅立叶变换。DFT 的计算在数字信号处理中非常有用，但是由于 DFT 的计算量太大，即使使用计算机也很难对问题进行实时处理，所以在很长的时间里并没有得到广泛应用。快速傅立叶变换（FFT）的提出使 DFT 的计算大大简化，从而使 DFT 的运算在实际中真正得到了广泛的应 用。现有的 FFT 可以分为两大类：时间抽取法（DIT-FFT）和频率抽取(DIP-FFT)，其

9、算法 思想基本一致，只是划分方式略有差异。现以按时间抽取(DIT)的基-2FFT 算法为例介绍 FFT 的基本原理并给出相应的 MATLAB实现方法。先设序列点数为 N=2L，L 为整数。如果不满足这个条件，可以认为的加上若干零值， 使之达到这一要求。将 N=2L 的序列 x(n)(n=0,1,，N-1)先按 n 的奇偶分成两组：x(2r)=x1(r)x(2r+1)=x2(r)r=0,1,N/2-1则其 DFT 可化为两个部分：nkNN1 12 rkN 1 22X (k ) =n=0x(n)WN =r = 0x(2r)WN+( 2 r +1) kr = 0x(2r + 1)WNN 1N 122

10、= x1 (r )WN+ WN x2 (r )WNrkkrkr = 02r = 02k= X1 (k ) + WN X 2 (k ), k = 0,1N 1这里 X1(k)和 X2(k)都是 N/2 点的 DFT，而 X(k)却有 N 个点，所以求 X(k)可分为前 N/2 和后N/2 点：k前一部分 X ( k ) =X1 ( k ) + WN X2 ( k )后一部 N N( N +k ) N分X 2 +k =X1 2 +k + WN 2X2 2 +k k= X1 (k ) WN X 2 (k )只要求出 k=0,1,2,N/2-1 区间内各点的 X1(k)和 X2(k)，即可求出 k=0

11、,1,2,N-1 区间全部 X(k)值，这样就大大节省了计算量。Matlab 提供一个称为ff t 函数来计算x(n) 的DFT ,ff t 函数是用机器语言,而不是以 Matlab 指令写成的,因此,它的计算速度非常快.如果N 是2 的整数幂,就采用基22FF T 算法, 如果N 不是2 的某个幂,就采用较慢的分裂基算法,如果N 是某个素数,那么ff t 函数就蜕化 为原始的DF T 算法.对于刚才的例子，用FFT实现如下：clear all; close all; clc; N=32; n=0:N-1;xn=cos(pi*n/6); y=fft(xn,N); figure(1); stem

12、(n,xn) figure(2); stem(n,y);得到的结果如下图：图 4.有限长序列图 5.FFT 运行结果由图可以看到结果与 DFT 的结果完全相同3数字滤波器的 MATLAB 仿真为了处理信号，必须设计和实现滤波器系统。滤波器的设计结果受滤波器的类型（即IIR 或 FIR）和其实现形式（结构）的影响。 一个数字滤波器可以用系统函数表示为：M k bk z NH ( z) = k =0 = Y ( z)k1 k =1a z kX ( z)由此式可得出表示输入输出关系的常系数线性差分方程为：NMy(n) = ak y(n k ) + bk x(n k )k =0k =0可见数字滤波器的

13、功能就是把输入序列 x(n)通过一定的运算变换成输出序列 y(n)。不同的运 算处理方法决定了滤波器实现结构的不同。3.1 IIR 滤波器的结构及 MATLAB 实现IIR 滤波器的特点（1） 单位冲激响应 h(n)是无限长的；（2） 系统函数 H(z)在有限 Z 平面（1|z|0 处收敛，极点全部在 z=0 处（稳定系统）；（3） 结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率 抽样结构）也包含有反馈的递归部分。FIR 滤波器实现的基本结构包括横截型结构、级联型结构、频率抽样型结构和快速卷积 型结构。3.2.1 横截型滤波器的结构设 FIR 滤波器的单位冲激相应 h(n

14、)为一个长度为 N 的序列，则滤波器系统函数为：N 1H ( z) = h(n) z nn= 0表示这一系统输入输出关系的差分方程为：N 1y(n) = h(m) x(n m)n= 0直接由差分方程得出的实现结构如下图【1】所示：X(n)z-1z-1z-1h(0)h(1)h(2)h(N-1)图 13 差分方程实现结构【1】3.2.2 级联型滤波器的结构2 k将 H(z)分解成实系数二阶因子的乘积形式：Ny(n)N 1 2 H ( z) = h(n) z n = 0 k+ 1kz 1 + z 2n= 0k =1这时 FIR 滤波器可用二阶节的级联结构来实现，每个二阶节用横截型结构实现如下图【1】

15、所示：B01B02X(n)y(n)z-1z-1B11B12z-1B21z-1B22图 14 FIR 滤波器级联结构【1】3.2.3 频率抽样型滤波器的结构若 FIR 滤波器的冲激响应为有限长（N 点）序列 h(n),则可化成1N 1H (z)=(1-z N )H (k ) K 1N k =0 1 WN z这种结构由两个部分级联而成，即：H (z)=H (k )( Nk K 1H c z)=1-z1 WN z其频率结构如下：x(n)z-1W0N-1W1N-ej2pi/NH(0)z-1H(1)z-1y(n)H(N-1)N-eWN-1j2pi(N-1)/Nz-1图 15 FIR 滤波器的频率抽样结构

16、现以一个 32 点的系统用 matlab 实现,其频率样本定义如下：|H(k)|=1(k=0,1,2)0.5(k=3)0(k=4,5,15)实现程序如下：clear all close all clc M=32;alpha=(M-1)/2; magHk=1,1,1,0.5,zeros(1,25),0.5,1,1; k1=0:15;k2=16:M-1;angHk=-alpha*2*pi/M*k1,alpha*2*pi/M*(M-k2); H=magHk.*exp(j*angHk);h=real(ifft(H,M); M=length(h); H=fft(h,M);magH=abs(H),phaH

17、=angle(H);if (M=2*floor(M/2) L=M/2-1;A1=1,-1,0;1,1,0;C1=real(H(1),real(H(L+2);elseendL=M/2-1; C1=real(H(1);k=1:L; B=zeros(L,2);A=ones(L,3);A(1:L,2)=-2*cos(2*pi*k/M); A=A;A1; B(1:L,1)=cos(phaH(2:L+1);B(1:L,2)=-cos(phaH(2:L+1)-(2*pi*k/M); C=2*magH(2:L+1),C1;结果如下图所示图 16 结果 A图 17 结果 B图 18 结果 C可以看出求出的 FI

18、R 系统频率收养响应结构只有 4 个增益系数非零。因此频率抽样形式为：H（z)=1 z 32320.9952 + 0.9952z 121 1.9616z 1 + z 20.9808 + 0.9808z 1+ 21 108478z 1 + z 20.9569 + 0.9569 z 1+ 21 1.6629 z 1 + z 2+1 1 z 1实现 32 点 FIR 系统频率抽样型机构时，注意到 H(0)=1，故一阶子系统不需要乘法运算， 而 3 个二阶子系统中的每个子系统需要 3 次乘法运算，所以当每个样本经过有的二阶子系统 时，共需要 9 次乘法运算，连同一阶子系统和二阶子系统总共需要 13 次

19、加法运算。而实现 相同点数的线性运算结构时，每个输出样本需要 16 次乘法和 31 次加法运算，因此 FIR 系统结构的频率抽样结构比线性相位结构复杂度更低，效率更高。3.2.4 快速卷积型滤波器的结构 快速卷积结构如下【1】X(n)L 点 DFTL 点 IDFTy(n)H(n)L 点 DFT图 19 FIR 快速卷积结构【1】这里的 DFT 和 IDFT 都采用 FFT 算法，当序列足够长时，比直接计算线性卷积要快的 多。4结论本文给出了数字信号处理中的各种算法的原理以及 MATLAB 模拟，并且结合滤波器的 基本结构实现了 MATLAB 的滤波器仿真。参考文献1 程佩青数字信号处理教程M北

20、京：清华大学出版社，20022 郭仕剑，王宝顺，贺志国，杨可心M。MATLAB7.x数字信号处理。北京：人民邮电出版社，2006。3 赵红怡,张常年. 数字信号处理及其MATLAB 实现M。北京:化学工业出版社,2002. 4 薛年喜. MATLAB 在数字信号处理中的应用M.北京:清华大学出版社,2003.Digital Signal Processing Simulation and Emulation Based on MATLABWang JingbinCollege of Information Science and Engineering, Shandong Agricultur

21、al University, Taian, PRC (271018)AbstractThis paper introduce basic principles of various algorithms in digital signal processing, including Ztransform, DFT Transform and its fast algorithm FFT, and digital filter design and realization including IIR and FIR, detailedly discusses principle, steps a

22、nd implementation method of signal transformation and digital filter design. Simulation and emulation using Matlab is given. At last experimental results is given, and explains the reasons of the different effect.Keywords: Z transform; Discrete Fourier Transform; digital filter; System Structure Simulation