第六章-matlab在信号与系统中的应用.ppt

第六章matlab在信号与系统中的应用 6 1连续信号和系统的MATLAB描述 提示 严格来讲 只有符号运算才能描述连续信号或系统 1 单位冲激函数 2 单位阶跃函数 3 复指数函数 6 1 1连续信号 Clear t0 0 tf 5 dt 0 05 t1 1 t t0 dt tf 1 单位脉冲信号 在t1 t0 t1 tf 处有一持续时间为dt 面积为1的脉冲信号 其余时间均为零 t t0 dt tf st length t n1 floor t1 t0 dt 求t1对应的样本序号x1 zeros 1 st 把全部信号先初始化为零x1 n1 1 dt 给出t1处的脉冲信号subplot 2 2 1 stairs t x1 grid 绘图 注意为何用stairs命令axis 0 5 0 1 1 dt m文件内容 q601 m 2 单位阶跃信号 信号从t0到tf 在t1 t0 t1 tf 前为0 到t1处有一跃变 以后为1 程序前几句即求t st n1的语句与上同 只把x1处改为x2x2 zeros 1 n1 1 ones 1 st n1 1 产生阶跃信号subplot 2 2 3 stairs t x2 grid 绘图axis 0 5 0 1 1 3 复数指数信号u 0 5 w 10 x3 exp u j w t subplot 2 2 2 plot t real x3 grid 绘图 subplot 2 2 4 plot t imag x3 grid 绘图 6 1 2LTI连续系统分析 n阶线性时不变连续系统的微分方程 已知y及其各阶导数的初始值为y 0 y 1 0 y n 1 0 求系统的零输入响应 解 方程的解为 p1 p2 pn是方程a1 n a2 n 1 an an 1 0的根 C1 Cn由y及其各阶导数的初始值来确定 时域 由不同输入时的响应来描述 变换域 传递函数 数字系统 差分方程 C1 C2 Cn y0y0 y 0 p1C1 p2C2 pnCn Dy0 即V C Y0其解为C V Y0式中V为范德蒙矩阵 在MATLAB的特殊矩阵库中有vander 调用方法 V vander p 求任意线性系统的零输入分量 a input 输入分母系数向量a a1 a2 书上取 3 5 7 1 n length a 1 disp 输入初始条件向量 长度应比a小一 Y0 y0 Dy0 D2y0 Y0 input Y0 书上分别取 1 0 0 0 1 0 0 0 1 p roots a 求系统的极点V rot90 vander p c V Y0 用范德蒙特矩阵求对应的系数dt input dt 书上取0 2 tf input tf 书上取8 t 0 dt tf y zeros 1 length t fork 1 ny y c k exp p k t end 将各分量叠加plot t y gridon LTI连续系统分析 零输入响应的m文件 q602 m 编程要点 1 设置已知条件 a a1 a2 Y0 y0 Dy0 D2y0 零输入时与b b1 b2 无关 可不设置 2 根据a a1 a2 求p p1 p2 pn p roots a 3 由p p1 p2 pn 计算范德蒙矩阵V vander p 4 由V及Y0计算C C1 Cn C V Y0 Y0 1 0 0 Y0 0 1 0 Y0 0 0 1 n阶LTI系统的冲激响应 n阶微分方程 写成系统函数为 冲击响应就是H s 的拉普拉斯反变换 可以把H s 展开为极点留数式 其反变换为 有输入 求任意线性系统的冲击函数a input 多项式分母系数向量a 书上取poly 0 1 2j 1 2j 2 5 b input 多项式分子系数向量b 书上取 8 3 1 r p residue b a 求留数disp 解析式h t r i exp p i t disp 给出时间数组t 0 dt tf dt input dt 给定时间数组tf input tf t 0 dt tf h zeros 1 length t h的初始化fori 1 length a 1 根数为a的长度减1h h r i exp p i t 叠加各根分量endplot t h grid LTI连续系统分析 冲激响应的m文件 q603 m 编程要点 1 设置系数值a a1 a2 an 1 b b0 b1 b2 bm 1 2 求留数 展开为极点留数式 r p residue b a 3 求冲激响应 h r i exp p i t 最后各分量叠加即为所求 课后练习 掌握LTI系统零状态响应的matlab求解方法 P 133例6 5 6 2傅里叶分析 图示的周期性方波 其傅里叶级数为分别计算直到9次谐波 并做图 用matlab演示谐波合成情况 信号特征分析 6 2 1信号的的分解 以方波为例 P 137例6 7 方波表为多阶正弦波之和t 0 01 2 pi 设定一个时间数组 有101个点y sin t figure 1 plot t y pause 频率为w 1 f 1 2 的正弦基波y sin t sin 3 t 3 plot t y pause 叠加三次谐波 用1 3 5 7 9次谐波叠加y sin t sin 3 t 3 sin 5 t 5 sin 7 t 7 sin 9 t 9 plot t y 傅里叶分析 方波信号分解 q607 m 为了绘制三维曲面 要把各次波形数据存为一个三维数组 因此必须重新定义y 重编程 y zeros 10 max size t x zeros size t fork 1 2 19x x sin k t k y k 1 2 x end 将各波形迭合绘出pause figure 1 plot t y 1 9 gridline 0 pi 0 5 pi 4 pi 4 加上方波幅度线及标注text pi 0 5 pi 4 pi 4 将各半波形绘成三维网格图 看出增加谐波阶次对方波逼近程度的影响halft ceil length t 2 pause figure 2 mesh t 1 halft 1 10 y 1 halft pause 只用正半周波形clc 6 2 2信号的频谱分析 以方波为例 将积分上下限定为0 10s 并将t分成N等份 用求和代替积分 这样 傅立叶变换式可写为求和可以用f t 行向量乘以e j tn列向量来实现 式中 t是t的增量 在程序中 用dt表示 P 143例6 11 求不同 处的F值 都用同一公式 这就可以利用MATLAB中的元素群运算 将 设为一个行数组 代入上式 则可写为 程序中 用w表示 其中 F是与w等长的行向量 t 是列向量 w是行向量 t w是一矩阵 其行数与t相同 列数与w相同 这样 此式就完成了傅里叶变换 类似地也可得到傅里叶逆变换表示式为 方波的频谱分析m文件 q611 m 演示 1 频域样本点数可取得与时域样本点数不同 2 若要求的频谱太宽 而时域样本点数又取得太少 会发生频率泄漏clear tf 10 N input 取时间分隔的点数N dt 10 N t 1 N dt 给出时间分割 时域采样f ones 1 N 2 zeros 1 N 2 给出信号 此处是方波 wf input 需求的频谱宽度wf Nf input 需求的频谱点数Nf w1 linspace 0 wf Nf dw wf Nf 1 F1 f exp j t w1 dt 求付利叶变换w fliplr w1 w1 2 Nf 补上负频率F fliplr F1 F1 2 Nf 补上负频率区的频谱subplot 1 2 1 plot t f linewidth 1 5 gridaxis 0 10 0 1 1 subplot 1 2 2 plot w abs F linewidth 1 5 grid 算得的时域信号波形及其频谱图如右 下图为采样周期较低时的情况 有明显的频率泄漏 N 256 wf 40 Nf 64 N 64 wf 40 Nf 256 6 2 3信号滤波 计算幅度为1 宽度为5s的矩形脉冲 同例6 11 通过下列滤波器的响应 1 理想低通滤波器 2 低通滤波器解 滤波器输出的频谱Y j F j H j 其时间响应y t 是Y j 的傅里叶反变换 P 145例6 12 1 理想低通滤波器的截止角频率 c 10 故只取F j 中 0 10的部分 用MATLAB语言表述 输出频率分量对应的 的下标数组为n2 find w wc w wc 其对应的频率数组为w2 w n2 频段内的频谱数组为F2 F n2 它就是滤波后的频谱数组Y2 其逆变换即y2 F2 exp j n2 t pi dw 2 三阶低通滤波器的频率响应滤波器的输出为 方波经过理想低通滤波器滤波q611g 调用上题的程序 但数据已确定F2f input 理想低通滤波器的带宽 1 秒 n2 find w F2f 对中段频谱求付利叶逆变换subplot 1 2 2 plot t f t y2 linewidth 1 5 grid 画出原波形及滤波后的波形 信号滤波 理想低通m文件 q612a m 方波经过三阶巴特沃斯低通滤波器后波形q611g 调用上题的程序 但数据已确定 截止频率为wc 10的三阶巴特沃斯低通滤波器转移函数为 H s 500 s 3 20s 2 200s 1000H freqs 500 1 20 200 1000 w 如用信号处理工具箱 语句为 b a butter wn 阶数 H freqs b a w Y H F figure 3 subplot 1 2 1 plot w abs Y linewidth 1 5 grid 画出滤波后的频谱y Y exp j w t pi dw 对中段频谱求付利叶逆变换subplot 1 2 2 plot t f t y linewidth 1 5 grid 画出原波形及滤波后的波形 理想低通 3阶巴特奥斯低通 6 4线性时不变系统的模型 1 连续系统 状态空间型设x为状态变量 u为输入 y为输出 系统的状态方程为 如果系统是n阶的 输入有nu个 输出有ny个 则A为n n阶 B为n nu阶 C为ny n阶 而D为ny nu阶矩阵 对单输入单输出 SISO 系统 ny nu 1 6 4 1模型的典型表达式 6 3离散信号和系统 放在数字信号处理部分讲 传递函数型 单输入单输出 SISO n阶系统的传递函数为知道分子系数矢量f f 1 f 2 f m 1 分母系数矢量g g 1 g 2 g n 1 就惟一地确定了系统的模型 注意系统的阶次n 而对物理可实现的系统 必有n m 零极增益型 对传递函数分子和分母进行因式分解 可得令z z 1 z 2 z m 为系统的零点矢量 p p 1 p 2 p n 为系统的极点矢量 k为系统增益 它是一个标量 可以看出 H s 有m个零点 n个极点 物理可实现系统的n m 系统的模型将由矢量z p及增益k惟一确定 故称为零极增益模型 零极增益模型通常用于描述SISO系统 并可以推广到MISO系统 极点留数型 如果零极增益模型中的极点都是单极点 将它分解为部分分式 可得其中p p 1 p 2 p n 仍为极点矢量 而r r 1 r 2 r n 为对应于各极点的留数矢量 p r两个矢量及常数h惟一地决定了系统的模型 典型模型表达式的比较 比较一下这四种情况下模型系数的总个数 假定都是SISO系统 阶数为n 则状态空间型有n2 2n 1个系数 传递函数型为m n 1个 不含g 1 注意 由于m n 系数的数目小于等于2n 1 零极增益型的系数个数为n m 1 而极点留数型为2n 1 因此 传递函数法的待定系数最少 而状态空间法的待定系数最多 这说明了状态空间法中有许多冗余的系数 事实上 同一个系统可以有无数个状态空间矩阵A B C D的组合来描述 其他描述方法则都是惟一的 2 离散系统 以上四种表示模型的方法可以全部推广至离散系统 只是将系数矩阵后面加小写字母d 便有 状态空间型x n 1 Adx n Bdu n y n Cdx n Ddu n 传递函数型零极增益型 极点留数型数字信号处理模型二阶环节型P 153页表6 1列出了连续和离散线性系统的模型式 6 4 2模型转换 MATLAB中各种模型转换的函数 模型转换函数调用格式 传递函数型 极点留数型 r p h residue f g 传递函数型 零极点增益型 z p k tf2zp f g f g zp2tf z p k f g residue r p h 传递函数型 状态空间型 A B C D tf2ss f g f g ss2tf A B C D 状态空间型 零极点增益型 z p k ss2zp A B C D A B C D zp2ss z p k 例子 已知描述系统的微分方程为 1 2 求出它的传递函数模型 零极增益模型 极点留数模型和状态空间模型 解 1 f 2 5 3 g 2 3 5 9 2 f 1 3 2 g 1 5 7 3 可用工具箱函数tf2zp tf2ss及residue函数求得 由传递函数模型转为其他模型f input 传递函数分子系数数组f f 1 f nf 书上取 2 5 3 g input 传递函数分母系数数组g g 1 g 2 g ng 书上取 2 3 5 9 printsys f g s disp 转为零极增益模型 z p k tf2zp f g disp 转为零极留数模型 r p h residue f g disp 转为状态空间模型 A B C D tf2ss f g printsys A B C D 模型转换m文件 q617 m 系统的串联 并联和反馈 系统的串联 由图所示 YB WBUB WBWAUA WU故W s WA s WB s 多项式相乘由卷积函数conv实现 其表示式为 f conv fA fB g conv gA gB 系统的并联Y WAU WBU WA WB U WU故W s WA s WB s f polyadd conv fA gB conv fB gA g conv gA gB p 159例6 19 故MATLAB表达式为 f conv fA gB g polyadd conv fA