MATLAB

MATLAB在数字信号处理中的应用一、离散时间信号与系统在数字信号处理中，分析的信号是离散时间信号即序列，对应的处理系统也为离散时间系统。1、基本信号的表示MATLAB中，可采用向量表示序列，由于MATLAB中对序列下标默认为从1开始递增，因此要表示离散信号，一般应采用两个向量分别对信号的自变量和因变量进行描述。如,n=-3～4，在MATLAB中表示为n=[-3,-2,-1,0,1,2,3,4];%自变量取值x=[2,1,-1,0,1,4,3,7];%因变量取值说明：向量可用方括号[]表示。当向量取值连续变化时可用冒号运算符“：”简化赋值过程，如可简化为n=[-3：4]或n=-3:4分号“；”表示不回显表达式的值。“%”表示其后内容为注释对象。利用MATLAB，还可对信号的波形进行描述，常采用的绘图语句有stem，plot，subplot,axis，title，xlabel，ylabel，gtext,holdon,holdoff,grid等。其中stem绘制离散图形；plot绘制连续图形；subplot用于绘制子图，应在stem或plot语句前调用；axis指定x和y轴的取值范围，用在stem或plot语句后；title标注图形名称，xlabel,ylabel分别标注x轴和y轴名称；gtext可将标注内容放置在鼠标点击处；holdon和holdoff用于控制对象绘制方式，是在原图上还是在新图上绘制；grid用于绘制网格。以上函数具体使用方法可通过在命令窗口中输入help加函数名进行查阅，如“helpstem”语句。2、序列操作1）信号的加减MATLAB中可用算术运算符“+”和“-”分别实现序列的加法和减法运算，此时两序列的长度必须相等，否则需通过补零法改变信号的长度。例：序列x1(n)、x2(n)，长度分别为n1和n2，其求和的MATLAB代码为：序列相加的长度的具有的长度的具有列的零值）行初始化（产生的长度%;212)(%;2))1))2max((&))2min((((21)(%;1))1))1max((&))1min((((1)(1%;12));(,1(1)(%));2max(),1max(max(:))2min(),1min(min(yyyxnyxnnnnfindyxnyxnnnnfindynlengthyynlengthzerosynynnnnn说明：函数min(n)可获取向量n的最小值，函数max(n)则获取向量n的最大值。函数length(n)可取得向量n的长度。函数find((na)&(nb))可获取向量n在范围anb内的取值。2）信号的乘除MATLAB中两序列对应值的乘、除可分别采用点乘“.*”和点除“./”运算符，运算时两序列的长度需相等，具体处理方法和序列求和差时的一致。例：序列x1(n)、x2(n)，长度分别为n1和n2，其求积的MATLAB代码为：序列相乘的长度的具有的长度的具有列的零值）行初始化（产生的长度%;2*.12)(%;2))1))2max((&))2min((((21)(%;1))1))1max((&))1min((((1)(1%;12));(,1(1)(%));2max(),1max(max(:))2min(),1min(min(yyyxnyxnnnnfindyxnyxnnnnfindynlengthyynlengthzerosynynnnnn3）移位序列x(n)的移位操作不影响向量x，只需对向量n的每个元素加或减去一个移位值。如移位序列等价为，对应的MATLAB语句为：n=n+n0;y=x;4）折叠在这个运算中，序列x(n)的每个样本都对n=0翻转，得到一个折叠后的序列y(n)，MATLAB中可由函数fliplr(x)实现。y(n)=x(-n)对应的代码为：y=fliplr(x);n=-fliplr(n);注意：翻转后向量n的取值为负值。5）样本和与积该运算是将序列中各样本值求和或求积，可分别采用函数sum(x)和prod(x)实现。具体应用时还可指定序列x参加运算的元素范围，即sum(x(n1:n2))prod(x(n1:n2))。6）信号能量MATLAB中，计算序列x(n)能量的可采用下列两种办法：法1：Ex=sum(x.*conj(x));%函数conj(x)是取复数x的共轭复数法2：Ex=sum(abs(x).^2);3、离散系统一个离散时间系统，可将一个序列x(n)变换为另一个序列y(n),即将输入信号处理成输出信号。常用离散系统为线性时不变系统（LTI）。1）计算系统零状态响应-卷积和线性时不变系统（LTI）的输入和输出关系可通过冲激响应h(n)表示，即y(n)为x(n)与h(n)的线性卷积。MATLAB提供卷积函数conv来计算两个有限长度序列的线性卷积，并假定两个序列都从n=0开始。调用方式为：y=conv(x,h)。例：x(n)=[1,2,3],h(n)=[2,2,3],n=[0,1,2],x(n)与h(n)的线性卷积代码为：x=[1,2,3];h=[2,2,3];y=conv(x,h)y=2613129对于起点和终点为任意值的序列进行线性卷积运算，可自行扩展另一个函数conv_m，具体形式为：function[y,ny]=conv_m(x,nx,h,nh)%[y,ny]=线性卷积结果%[x,nx]=第一个信号%[h,nh]=第二个信号nyb=nx(1)+nh(1);%两向量起点之和nye=nx(length(x))+nh(length(h));%两向量长度之和ny=[nyb:nye];y=conv(x,h);2）系统差分方程-时域解系统差分方程的解在MATLAB中可通过调用滤波器函数filter（a,b,x）完成，其中a=[a0,a1,…,aN](a0≠0)，b=[b0,b1,…,bM]分别对应差分方程中y和x的各项系数，x为输入信号序列，即将差分方程看成一个滤波器，x作为滤波器的输入信号，y为滤波器的输出值。nmnxbknyaMmmNkk,)()(003）系统稳定性LTI离散系统稳定性的判别可采用两种方法：方法1：对冲激响应序列h(n)绝对值求和，判别其是否为有限值，若有限，为稳定系统。对应语句为：sum(abs(h))，其中函数abs是取序列h(n)的绝对值。方法2：利用系统函数H(z)的极点分布进行判断，若极点均在单位圆内，则为稳定系统。系统函数H（z）的极点也即为差分方程中y项对应的各系数[a0,a1,…,aN]构成的多项式的极点，而极点可采用求根函数roots获取。二、离散时间傅立叶分析离散信号的频谱可采用离散时间傅立叶变换（DTFT）进行描述。对一离散信号x(n)，其离散时间傅立叶变换式为:()()jjnnXexne频谱为周期函数，周期为2π。分析频谱时，为便于读数，频率轴常以π弧度为单位，同时由于频谱的周期性，只需对一个周期（即或[-，]）进行描述。当x(n)为实值时，为共轭对称函数，此时只需分析半周期[0，]频谱。对于无限长信号x(n)，MATLAB不能直接计算，但也可估计在[0，]的值。1、频谱分析常用命令为描述信号的频谱，在MATLAB中常采用的符号及命令如下：*pi：表示π*函数exp（a+bj）：可实现复指数e(a+bj)的调用*函数angle（x）：获取复数x的相角，可描述信号频谱的相位特性*函数abs（x）：获取复数x的模，可描述信号频谱的幅度特性*函数real（x）：获取复数x的实部*函数image（x）：获取复数x的虚部*函数conj(x):取复数x的共轭复数2、离散时间傅立叶变换（DTFT）例：将[0,π]分为501个等间隔的点，确定x(n)=(0.5)nu(n)的离散傅立叶变换，并画出其模、相角实部、虚部曲线。3、线性时不变（LTI）系统频域表示任意线性时不变系统都可由冲激响应h(n)来表示，相应地在频域中可用频率响应即h(n)的离散时间傅立叶变换H(ej)来表示。由差分方程表示的LTI系统NlmMmlmnxblnyany10)()()(其频率响应为ljlNlMmmjmjeaebeH101)(若在[0，π]等间隔分成K+1点，可计算出对应离散点频率的频谱，具体的MATLAB代码为：4、模拟信号的采样和重构三、Z变换1、Z变换时域信号的卷积运算对应Z域的乘积，而Z域的乘积可转换为多项式的乘积，即而可使用MATLAB的卷积conv函数实现。2、逆Z变换逆Z变换的中心思想：若X(z)是z-1的有理函数，可用部分分式分解的方法将其变成简单因式项（一阶）的和，然后由z变换表写出对应因式项的序列。对有理式3、Z域系统表示离散系统可采用系统函数H(z)来描述，在MATLAB中，对于一确定的H(z)，可分析其零极点及幅度、相位响应。1）零极点MATLAB中提供函数p=roots(c)（poly的逆函数）可求得H(z)分子和分母多项式对应的根，即而可确定H(z)的零极点，其中c为多项式对应的按z-1递增排列的系数序列。函数zplane(z,p)则可画出对应的零极点图，其中z和p分别为零点和极点序列。2）幅度与相位响应MATLAB提供了一个等效的函数freqz，其最简单的调用形式为：[H,w]=freqz(b,a,N)，其中b和a分别表示分子和分母多项式的按z-1递增排列的系数向量，此函数返回该系统的N点频率向量w和N点复数频率向量H。3）系统的单位冲激响应和单位阶跃响应离散系统H(z)的单位冲激响应和单位阶跃响应可采用两种方式。法1：用函数y=filter(b,a,x)，其中b和a分别为H(z)分子和分母多项式对应的按z-1递增排列的系数序列；x可描述为对应的单位冲激序列和单位阶跃序列，y为对应的响应，然后再画其响应曲线。法2：直接调用函数dimpulse(b,a)和dstep(b,a)可绘出系统的单位冲激响应和单位阶跃响应曲线，其中b和a与法1中的一样。4、差分方程求解MATLAB中提供求xic的函数，即xic=filtic(b,a,Y,X)，其中b和a分别为差分方程中x和y对应的系数序列，b=[b0,b1,…,bM]，a=[1,a1,…,aN]；Y和X是初始条件数组,Y=[y(-1),y(-2),…,y(-N)]，X=[x(-1),x(-2),…,x(-M)]。四、离散傅立叶变换（DFT）1、离散傅立叶级数（DFS）2、离散傅立叶变换（DFT）及其性质1）循环折叠性序列x(n)的循环折叠为11),(0),0())((NnnNxnxnxN在MATLAB中，可由x=x(mod(-n,N)+1)得到信号x的折叠循环，其中函数mod是取余操作。2）序列的循环移位)())(()()(~nRmnxnRmnxNNN3）序列的循环卷积序列x1和x2做N点循环卷积时，两序列的长度必须大于或等于N，实现方法为y(n)=sum(x1(m)*x2((n-m)modN))3、快速傅立叶变换（FFT）MATLAB中提供一实现快速傅立叶变换的函数fft，其调用形式之一为X=fft(x,n)，x为输入时域序列，n为FFT点数，X为序列x的频谱。另对于快速傅立叶反变换，MATLAB也提供一内部函数x=ifft(X,n)，x、X和n与函数fft中的一样。五、IIR数字滤波器1、模拟滤波器1）巴特沃斯滤波器MATLAB提供函数[z,p,k]=buttap(N)，来设计N阶归一化巴特沃斯模拟低通滤波器，它返回数组z和p（零点和极点）以及增益k。对于非归一化巴特沃斯滤波器可进行扩展。2）切比雪夫滤波器MATLAB提供两个切比雪夫滤波器设计函数，其一为[z,p,k]=cheb1ap(N,Rp)，可设计阶数为N，通带波动为Rp的归一化切比雪夫Ⅰ型低通滤波器，数组z和p分别返回零点和极点，k为增益；另一函数[z,p,k]=cheb2ap(N,Rs)，可设计阶数为N，阻带衰减为Rs的