数字信号处理 area=-1.07578023141031 quad使用 Simpson递归方法,quad8使用 Newton- costes递归方法进行数值积分。为了获得 更精确的结果,它们在所需的区间都计算被积函数。quad8比quad更精确。这两个函数的 使用和 fzero相同。 微分 与积分相反,数值微分十分困难。积分描述了一个函数的整体或宏观的性质,所以积分 对函数的形状在小范围的变化不敏感:而微分则描述了函数在一点处的斜率,是函数的微观 性质,它对函数的微小变化十分敏感,函数的很小的变化,容易产生相邻点斜率的巨大变化。 由于数值微分的固有困难,应当尽量避免使用数值微分,尤其是树言数据获得的数据微分 如果迫切需要,最好先将试验数据进行最小二乘拟合伙这三次样条拟合,然后对拟合函数进 行微分 MATLAB提供了一个有限插分函数dif,可以做数值微分。在应用中需要注意的是,插分后 输出数组比原数组少了一个元素。实际表明,使用有限插分近似将放大噪声,导致极差的结 5、FFT变换 FFT即快速傅立叶变换,是数据分析的基本方法,它是向量X的离散傅立叶变换由基2 的快速变换算法来计算的。如果ⅹ的长度不是精确的2次幂则后面使用0填充,mx)是向 量ⅹ的离散傅立叶变换的逆变换。 在频率轴商会值FFT曲线,要明确FFT结果与实际频率点的关系。设N个数据点,采样频 率为fs,则 Nyquist频率或n=N2+1点与实际频率的关系为 f(num-1)"fs/n 例如:绘制向量x的频谱 =ff n= ength(y)2,%FFT变换是对称的,取前半部分 f=fs*(O: n-1)/n plot(f, abs(y)); plot(f, (180/pi)*unwrap(atan2(imag(y), real(y))); 需要注意的是ft结果为复数矩阵,为了得到幅频特性,可使用abs函数,使用atan2得到相 角,由于有的系统的相角可能大于180°,而相角函数值域在-180180之间,需要使用 unwrap 函数展开折叠的相角,从而得到相频特性。 6、数据分析函数 MA∏LAB提供了很多有用的数据分析函数,这些函数数量极大,需要慢慢发掘。下面 给出一些常用的函数: 函数名 含义
数字信号处理 6 area =-1.07578023141031 quad 使用 Simpson 递归方法，quad8 使用 Newton-costes 递归方法进行数值积分。为了获得 更精确的结果，它们在所需的区间都计算被积函数。quad8 比 quad 更精确。这两个函数的 使用和 fzero 相同。 4、微分 与积分相反，数值微分十分困难。积分描述了一个函数的整体或宏观的性质，所以积分 对函数的形状在小范围的变化不敏感；而微分则描述了函数在一点处的斜率，是函数的微观 性质，它对函数的微小变化十分敏感，函数的很小的变化，容易产生相邻点斜率的巨大变化。 由于数值微分的固有困难，应当尽量避免使用数值微分，尤其是树言数据获得的数据微分。 如果迫切需要，最好先将试验数据进行最小二乘拟合伙这三次样条拟合，然后对拟合函数进 行微分。 MATLAB 提供了一个有限插分函数 diff，可以做数值微分。在应用中需要注意的是，插分后 输出数组比原数组少了一个元素。实际表明，使用有限插分近似将放大噪声，导致极差的结 果。 5、FFT 变换 FFT 即快速傅立叶变换，是数据分析的基本方法，它是向量 X 的离散傅立叶变换由基 2 的快速变换算法来计算的。如果 X 的长度不是精确的 2 次幂则后面使用 0 填充，lfft(x)是向 量 x 的离散傅立叶变换的逆变换。 在频率轴商会值 FFT 曲线，要明确 FFT 结果与实际频率点的关系。设 N 个数据点，采样频 率为 fs，则 Nyquist 频率或 n=N/2+1 点与实际频率的关系为： f=(num-1)*fs/n 例如：绘制向量 x 的频谱： y=fft(x); n=length(y)/2; %FFT 变换是对称的，取前半部分 y=y(1:n); f=fs*(0:n-1)/n; plot(f,abs(y)); figure; plot(f,(180/pi)*unwrap(atan2(imag(y),real(y)))); 需要注意的是 fft 结果为复数矩阵，为了得到幅频特性，可使用 abs 函数，使用 atan2 得到相 角，由于有的系统的相角可能大于 1800，而相角函数值域在-1800~1800 之间，需要使用 unwrap 函数展开折叠的相角，从而得到相频特性。 6、数据分析函数 MATLAB 提供了很多有用的数据分析函数，这些函数数量极大，需要慢慢发掘。下面 给出一些常用的函数： 函数名 含义
数字信号处理 max 最大值 min 最小值 mean 均值 标准方差 sedan 排序分类 元素的总和 元素的乘积 cumrod 元素的累积 cumsum 元素的累加和 Dift 差分函数 Hist 直方图 Tabel 列表 Corr 互相关矩阵 协方差矩阵 Find 查找逻辑 这些函数可以灵活应用,比如在生成Bode图的频率序列时,使用 logspace直接得到了对数 等间距的序列,如果想插入几点自己需要的频率点,可以使用下面的命令 p=logspace(0, 2, 100) [1.5345678] t =sort(x fD) 这样生成的频率点集合xf中不但含有1到100的对数等间隔频率点,还包含了频率f中的 四个频率点,这个频率在绘制有尖锐幅频特性的bode图时十分有用
数字信号处理 7 max 最大值 min 最小值 mean 均值 std 标准方差 sedian 中值 sort 排序分类 sum 元素的总和 prod 元素的乘积 cumrod 元素的累积 cumsum 元素的累加和 Diff 差分函数 Hist 直方图 Tabel 列表 Corr 互相关矩阵 Cov 协方差矩阵 Find 查找逻辑 这些函数可以灵活应用，比如在生成 Bode 图的频率序列时，使用 logspace 直接得到了对数 等间距的序列，如果想插入几点自己需要的频率点，可以使用下面的命令： p=logspace(0,2,100); r=[1.5 3.4 56 78]; t=sort([x f]); 这样生成的频率点集合 xf 中不但含有 1 到 100 的对数等间隔频率点，还包含了频率 f 中的 四个频率点，这个频率在绘制有尖锐幅频特性的 bode 图时十分有用
数字信号处理 四、常微分方程数值解 控制系统的模型很大一部分是常微分方程形式的,求取它们的解析解比较困难,一般使 用数值解。常微分方程数值解一般使用逐步积分的方法实现, Runge-Kutta法是应用最多的 种微分方程数值解的方法。 MATLAB提供了两种 Runge-Kutta法函数 It, x]=ode23(fun? tO, tf, xO, tol, trace) It, x]F=ioe45(fun?, to, tf, xo, tol, trace) 这两种方法格式相同。其中xfun为定义的常微分方程函数名,该函数必须以为输出,以 t、x为输入。输入变量t0、t为积分的启始和中止时间,单位是秒。x0为初始的状态向量。 tol控制结果的精度,可以缺省。一般来说,ode45比ode3运算速度快一些 考虑描述经典的 Var der pol微分方程 a24(1-x x dt 求解中将这高阶方程等价变换为一阶方程组,重新定义变量,变换为 令x1=xx2=d/dt 则dldt=x2 dx2/dt=u(1-x12)x2-x1 把它写成函数: function yp=vdp(t, x) yp(2)=2+(1-x(1)^2)*x(2}x(1)%令 在命令行求解这个方程 [tx]=ode45(dp?0,20,[1;1); plot(t, x(, 1 ), t, x(, 2) %画出x和dx/dt的时域波形+图x和d/dt的时域波形
数字信号处理 8 四、常微分方程数值解 控制系统的模型很大一部分是常微分方程形式的，求取它们的解析解比较困难，一般使 用数值解。常微分方程数值解一般使用逐步积分的方法实现，Runge-Kutta 法是应用最多的 一种微分方程数值解的方法。MATLAB 提供了两种 Runge-Kutta 法函数： [t,x]=ode23(fun?t0,tf,x0,tol,trace) [t,x]=ioe45(fun? ,t0,tf,x0,tol,trace) 这两种方法格式相同。其中 xfun 为定义的常微分方程函数名，该函数必须以 为输出，以 t、x 为输入。输入变量 t0、tf 为积分的启始和中止时间，单位是秒。x0 为初始的状态向量。 tol 控制结果的精度，可以缺省。一般来说，ode45 比 ode23 运算速度快一些。 考虑描述经典的 Var der Pol 微分方程 求解中将这高阶方程等价变换为一阶方程组，重新定义变量，变换为： 令 x1=x x2=dx/dt 则 dx1/dt=x2 dx2/dt=u(1-x12 )x2-x1 把它写成函数： function yp=vdp(t,x) yp(1)=x(2); yp(2)=2*(1-x(1)^2)*x(2)-x(1)； %令 u=2 在命令行求解这个方程： [t,x]=ode45(dp?0,20,[1 ;1]_); plot(t,x(:,1),t,x(:,2)); %画出 x 和 dx/dt 的时域波形+图 x 和 dx/dt 的时域波形
数字信号处理 第二部分: 一数字信号处理与 MATLAB语言基础 11简介 MATLA B是一套用于科学工程计算的可视化高性能语言与软件环境它集数值分析,矩阵 运算,信号处理和图形显示于一体构成一个界面友好的用户环境在这个环境中问题与求解 都能方便地以数学的语言(主要是矩阵形式)或图形方式表示出来 MATLA B的含义是矩阵实验室( Matrix Laboratory),该软件是一个交互式系统,其基本元素 是无需定义维数的矩阵萦研制的初衷主要是方便矩阵的存取但以过几十年的扩充和完善,已己 成为和各类科学研究与工程应用中的标准工具,其典型应用有数值分析,算法分析,自动控制 数字信号处理,图像处理以及模型仿真等等. MATLAB包括了被称为 Toolbox(工具箱)的各类应用问题的求解工具,本书着重介绍其中 的信号处理工具箱 MATLAB信号处理工具箱包含了各类经典的和现代的数字信号处理技术 是一个非常优秀的算法研究与辅助设计工具,它在语音信号处理,生物医学工程,实时控制和 雷达信号处理等多个研究领域得到成功的应用 本章将简要地介绍 MATLAB的基本知识为后面章节中信号处理的理解和学习打好基础, 本书的所有介绍都以 MATLAB6.1版本为基础如使用其他版本,可能有所不同 12 MATLAB快速入门 1.2.1启动 MATLAB有两种方法 (1)在 WINDOW9598/200NT40下,单击任务条上的"开始"按钮,选择"程序"”子菜单,再 选择" MATLAB"程序组,最后单击 MATLAB6.1"菜单项,就可以启动 MATLAB,由于在安装了 MATLAB软件后会自动在桌面上添加一个图标" MATLAB6.",故也可以直接双击该图标以 MATLAB ()在 WINDOWS95/982000NT40下,单击任务条上”开始"按钮,选择"运行"菜单项,然后 在"运行"对话框内输入" MATLAB"命令,然后点击"确定"按钮,就可以 MATLAB 以上两种操作的结果都会出现一个 MATLAB的命令窗口,其命令提示符为"? 122 MATLAB命令窗口与基本的矩阵操作 MA∏LAB命令窗口中,在 MATLAB提示符下要输入一个4*4矩阵可以按如下方式输入 命令 X=[1234;543234567654] 或 X=[1234 5432 456 7654] 以上两种输入方式的效果是一样的,命令未尾的分号用于禁止显示该命令的执行结 果如果去掉分号 MATLAB就会显示如下的命令的执行结果 X=[123454323456;7654 234
数字信号处理 9 第二部分： 一 数字信号处理与 MATLAB 语言基础 1.1 简介 MATLAB 是一套用于科学工程计算的可视化高性能语言与软件环境.它集数值分析,矩阵 运算,信号处理和图形显示于一体,构成一个界面友好的用户环境,在这个环境中,问题与求解 都能方便地以数学的语言(主要是矩阵形式)或图形方式表示出来. MATLAB 的含义是矩阵实验室(Matrix Laboratory),该软件是一个交互式系统,其基本元素 是无需定义维数的矩阵萦研制的初衷主要是方便矩阵的存取,但以过几十年的扩充和完善,已 成为和各类科学研究与工程应用中的标准工具,其典型应用有:数值分析,算法分析,自动控制, 数字信号处理,图像处理以及模型仿真等等. MATLAB 包括了被称为 Toolbox(工具箱)的各类应用问题的求解工具,本书着重介绍其中 的信号处理工具箱.MATLAB 信号处理工具箱包含了各类经典的和现代的数字信号处理技术, 是一个非常优秀的算法研究与辅助设计工具,它在语音信号处理,生物医学工程,实时控制和 雷达信号处理等多个研究领域得到成功的应用. 本章将简要地介绍 MATLAB 的基本知识,为后面章节中信号处理的理解和学习打好基础, 本书的所有介绍都以 MATLAB6.1 版本为基础,如使用其他版本,可能有所不同. 1.2.1 MATLAB 快速入门 1.2.1 启动 MATLAB 有两种方法 (1)在 WINDOW95/98/2000/NT4.0 下,单击任务条上的"开始"按钮,选择"程序"子菜单,再 选择"MATLAB"程序组,最后单击"MATLAB6.1"菜单项,就可以启动 MATLAB,由于在安装了 MATLAB 软件后,会自动在桌面上添加一个图标"MATLAB6.1",故也可以直接双击该图标以 MATLAB. (2)在 WINDOWS95/98/2000/NT4.0 下,单击任务条上"开始"按钮,选择"运行"菜单项,然后 在"运行"对话框内输入"MATLAB"命令,然后点击"确定"按钮,就可以 MATLAB. 以上两种操作的结果都会出现一个 MATLAB 的命令窗口,其命令提示符为"?". 1.2.2 MATLAB 命令窗口与基本的矩阵操作 MATLAB 命令窗口中, 在 MATLAB 提示符下要输入一个 4*4 矩阵可以按如下 方式输入 命令: X=[1 2 3 4; 5 4 3 2 ;3 4 5 6;7 6 5 4]; 或 X=[1 2 3 4 5 4 3 2 3 4 5 6 7 6 5 4 ]; 以上两种输入方式的效果是一样的,命令未尾的分号用于禁止显示该命令的执行结 果,如果去掉分号,MATLAB 就会显示如下的命令的执行结果: X=[1 2 3 4;5 4 3 2;3 4 5 6;7 6 5 4] X= 1 2 3 4
数字信号处理 5432 3456 要察看X变量的内容,也可以直接输入法X,按回车 MATLAB的显示与上面的结果似如果要 引用X的某几行(以第二,三行为例,则 X([2,3],) ans- 5432 3456 如果要引用X的某几列(以第二,三,四列为例),则 234 456 654 要表示一个递增右递减的序列可以用冒号,其格式为 此式表示产生从n到m且步长为s的一系列值当s省略时,默认步长为1,在前面的取x 的列的例子中我们己经用到了该方法现在如果要倒过来取列为43,2,且行为3,2,1的值,那么 命令如下 x(3:-1:1,4:-1:2) ans=6 5 4 234 432 求一个矩阵的转置,命令如下 ans- 1537 2446 3355 4264 要求矩阵X的逆命令如下(本例中的X的矩阵已为奇异阵求逆的结果已无意义) MATLAB除了基本的标量运算外,还有各种矩阵运算与数组运算 矩阵与矩阵,矩阵与标量的加法减法,和乘法在形式上与标量运算类似只是当矩阵与矩阵 作相应的运算时,要求维数相互匹配有时需要对矩阵进行必要的转置操作,这属于线性代数 的基础知 MATLAAB为矩阵的除法提供了两种运算,左除()和右除(,如果a为一非奇异矩阵则ab 和bla分别等价于 a/b=inv(a)b b/a=b*inv(a)
数字信号处理 10 5 4 3 2 3 4 5 6 7 6 5 4 要察看 X 变量的内容,也可以直接输入法 X,按回车,MATLAB 的显示与上面的结果似,如果要 引用 X 的某几行(以第二,三行为例),则 X([2,3],:) ans= 5 4 3 2 3 4 5 6 如果要引用 X 的某几列(以第二,三,四列为例),则 x(：，2：4) ans= 2 3 4 4 3 2 4 5 6 6 5 4 要表示一个递增右递减的序列,可以用冒号,其格式为 n：s：m 此式表示产生从 n 到 m 且步长为 s 的一系列值,当 s 省略时,默认步长为 1,在前面的取 x 的列的例子中我们已经用到了该方法,现在,如果要倒过来取列为 4,3,2,且行为 3,2,1 的值,那么 命令如下: x(3：-1：1,4：-1：2) ans= 6 5 4 2 3 4 4 3 2 求一个矩阵的转置,命令如下: x' ans= 1 5 3 7 2 4 4 6 3 3 5 5 4 2 6 4 要求矩阵 X 的逆,命令如下(本例中的 X 的矩阵已为奇异阵,求逆的结果已无意义): y=inv(x) MATLAB 除了基本的标量运算外,还有各种矩阵运算与数组运算. 矩阵与矩阵,矩阵与标量的加法,减法,和乘法在形式上与标量运算类似,只是当矩阵与矩阵 作相应的运算时,要求维数相互匹配,有时需要对矩阵进行必要的转置操作,这属于线性代数 的基础知识. MATLAAB 为矩阵的除法提供了两种运算,左除(\)和右除(/),如果 a 为一非奇异矩阵,则 a\b 和 b\a 分别等价于: a/b=inv(a)*b b/a=b*inv(a)