第4章MATLAB在信号处理中的应用_Matlab及在电子信息类课程中的应用.ppt

第4章MATLAB在信号处理中的应用 4 1信号及其表示4 2信号的基本运算4 3信号的能量和功率4 4线性时不变系统4 5线性时不变系统的响应4 6线性时不变系统的频率响应4 7傅里叶 Fourier 变换4 8IIR数字滤波器的设计方法4 9FIR数字滤波器设计 4 1信号及其表示 4 1 1连续时间信号的表示连续时间信号 时间变化连续 如y x t 离散时间信号 序列 时间离散 如x nT x t t nT 例4 1 用MATLAB命令绘出连续时间信号x t e 0 707tsin 2t 3 关于t的曲线 t的范围为 0 30s 并以0 1s递增 4 1 2工具箱中的信号产生函数 4 1信号及其表示 1 Sawtooth 产生锯齿波或三角波信号 格式一 x sawtooth t 功能 产生周期为2 幅值从 1 1的锯齿波 在2 的整数倍处值为 1 1 这一段波形斜率为1 格式二 x sawtooth t width 功能 产生三角波 width在0 1之间 例4 2 产生周期为0 02的三角波 2 Square 产生方波信号 格式一 x square t 功能 产生周期为2 幅值从 1 1的方波 格式二 x square t duty 功能 产生制定周期的方波 duty为正半周期比例 例4 3 产生周期为0 02的方波 4 1信号及其表示 3 Sinc 产生sinc函数波形函数sinc 格式一 y sinc x 功能 sinc x 用于计算sinc函数 例4 4 产生sinc函数波形 4 Rectpuls 产生非周期方波信号函数 格式一 y rectpuls t 功能 产生非周期方波信号 方波的宽度为时间轴的一半 格式二 y rectpuls t w 功能 产生指定宽度为w的非周期方波 例4 5 产生非周期方波 4 1信号及其表示 5 Tripuls 产生非周期三角波信号函数Tripuls 格式一 y tripuls t 功能 产生非周期三角波信号 三角波的宽度为时间轴的一半 格式二 y tripuls t w s 功能 产生指定宽度为w的非周期方波 斜率为s 1 s 1 例4 6 6 Chirp 产生线性调频信号函数chirp 格式 y chirp t f0 t1 f1 功能 产生一个线性扫频 频率随时间线性变化 信号 其时间轴的设置由数组t定义 时刻0的瞬时频率为f0 时刻t1的瞬时频率为f1 默认情况下 f0 0Hz t1 1 f1 100Hz 例4 7 绘制一个线性调频信号 4 1 3离散时间信号的表示 在MATLAB中 离散时间信号x n 的表示 需用一个向量x表示序列幅值 用另一个等长的定位时间变量n 才能完整地表示一个序列 例4 10 绘制离散时间信号的棒状图 其中x 1 1 x 0 1 x 1 2 x 2 1 x 3 0 x 4 1 MATLAB源程序为 n 3 5 定位时间变量x 0 0 1 1 2 1 1 0 0 stem n x grid 绘制棒状图line 3 5 0 0 画x轴线xlabel n ylabel x n 运行结果如图4 11所示 图4 11离散时间信号图形 4 1 4几种常用离散时间信号的表示 1 单位脉冲序列 直接实现 x zeros 1 N x 1 n0 1 函数实现 利用单位脉冲序列 n n0 的生成函数impseq 即function x n impseq n0 ns nf n ns nf x n n0 0 2 单位阶跃序列 直接实现 n ns nf x n n0 0 函数实现 利用单位阶跃序列u n n0 的生成函数stepseq 即function x n stepseq n0 ns nf n ns nf x n n0 0 3 实指数序列 直接实现 n ns nf x a n 函数实现 利用实指数序列x n an的生成函数rexpseq 即function x n rexpseq a ns nf n ns nf x a n 4 复指数序列 直接实现 n ns nf x exp sigema jw n 函数实现 利用复指数序列x n e j n的生成函数cexpseq 即function x n cexpseq sigema w ns nf n ns nf x exp sigema j w n 5 正 余 弦序列 直接实现 n ns nf x cos w n sita 函数实现 利用正 余 弦序列x n cos n 的生成函数cosseq 即function x n cosseq w ns nf sita n ns nf x cos w n sita 4 2信号的基本运算 4 2 1信号的相加与相乘y n x1 n x2 n y n x1 n x2 n MATLAB实现 y x1 x2 y x1 x2例4 11 信号的相加与相乘 4 2 2序列移位与周期延拓运算 序列移位 y n x n m MATLAB实现 y x ny nx m序列周期延拓 y n x n M MATLAB实现 ny nxs nxf y x mod ny M 1 例4 12 4 2 4两序列的卷积运算 两序列卷积运算 MATLAB实现 y conv x1 x2 序列x1 n 和x2 n 必须长度有限 4 2 5两序列的相关运算 两序列相关运算 MATLAB实现 y xcorr x1 x2 4 2 3序列翻褶与序列累加运算 序列翻褶 y n x n MATLAB可实现 y fliplr x 序列累加的数学描述为 MATLAB实现 y cumsum x 例4 13求序列x n 3e 0 2n的翻褶序列y n x n 4 3信号的能量和功率 1 信号能量 数字定义 MATLAB实现 E sum x conj x 或E sum abs x 2 数字定义 2 信号功率 MATLAB实现 P sum x conj x N 或E sum abs x 2 N 例4 16非周期三角波信号能量的MATLAB计算和余弦信号的平均功率计算 4 4线性时不变系统 4 4 1系统的描述 1 常系数线性微分 差分方程 2 系统传递函数 3 零 极点增益模型 连续系统 连续系统 离散系统 离散系统 4 极点留数模型 离散系统 连续系统 5 二次分式模型 连续系统 离散系统 6 状态空间模型 连续系统 离散系统 4 4 2系统模型的转换函数 在MATLAB中 用sos ss tf zp分别表示二次分式模型 状态空间模型 传递函数模型和零 极点增益模型 其中sos表示二次分式 g为比例系数 sos为L 6的矩阵 即 4 15 1 ss2tf函数格式 num den ss2tf A B C D iu 功能 将指定输入量iu的线性系统 A B C D 转换为传递函数模型 num den 2 zp2tf函数格式 num den zp2tf z p k 功能 将给定系统的零 极点增益模型转换为传递函数模型 z p k分别为零点列向量 极点列向量和增益系数 线性系统模型的变换函数 例4 19 求离散时间系统 的零 极点向量和增益系数 在命令窗口输入 num 2 3 den 1 0 4 1 num den eqtflength num den 使长度相等 z p k tf2zp num den 屏幕显示为z 0 1 5000p 0 2000 0 9798i 0 2000 0 9798ik 2 4 4 3系统互联与系统结构 MATLAB实现函数series 格式 A B C D series A1 B1 C1 D1 A2 B2 C2 D2 或 num den series num1 den1 num2 den2 将系统1 系统2级联 可得到级联连接的传递函数形式为 1 系统的级联 MATLAB实现函数parallel 格式 A B C D parallel A1 B1 C1 D1 A2 B2 C2 D2 或 num den parallel num1 den1 num2 den2 2 系统的并联 将系统1 系统2并联 可得到并联连接的传递函数形式为 3 两个系统的反馈连接函数feedback格式 A B C D feedback A1 B1 C1 D1 A2 B2 C2 D2 sign 或 num den feedback num1 den1 num2 den2 sign 将系统1和系统2进行反馈连接 sign表示反馈方式 默认值为 1 当sig 1时表示正反馈 当sig 1时表示负反馈 例4 20 求两个单输入单输出子系统 的级联 并联和反馈后系统的传递函数 MATLAB源程序为 num1 1 den1 1 1 系统1num2 2 den2 1 2 系统2 nums dens series num1 den1 num2 den2 实现两个系统级联 nump denp parallel num1 den1 num2 den2 实现两个系统并联 numf denf feedback num1 den1 num2 den2 实现两个系统反馈程序运行结果为 nums 002 dens 132nump 034 denp 132numf 012 denf 134因此 各系统的传递函数分别为 4 5线性时不变系统的响应 4 5 1线性时不变系统的时域响应 1 连续LTI系统的响应 2 离散LTI系统的响应 用MATLAB中的卷积函数conv 来实现 例4 24P147 用MATLAB中的卷积函数conv 来实现 例4 25已知LTI离散系统的单位冲击响应为 h n 0 5n n 0 1 2 14 求输入信号序列x n 1 5 n 4 的系统响应 格式 y x lsim a b c d u t 功能 返回连续LTI系统 2 对任意输入的离散LTI系统响应函数dlsim 格式 y x dlsim a b c d u 功能 返回离散LTI系统 对任意输入时系统的输出响应y和状态记录x 其中u给出每个输入的时序列 一般情况下u为一个矩阵 t用于指定仿真的时间轴 它应为等间隔 例4 26P150 对输入序列u的响应y和状态记录x 例4 27P151 3 时域响应函数 1 对任意输入的连续LTI系统响应函数lsim 4 5 2LTI系统的单位冲激响应 1 求连续LTI系统的单位冲激响应函数impulse 格式 Y T impulse sys 或impulse sys 功能 返回系统的响应Y和时间向量T 自动选择仿真的时间范围 其中sys可为系统传递函数 零极增益模型或状态空间模型 例4 28P152 2 求离散系统的单位冲激响应函数dimpulse 格式 y x dimpulse num den 功能 返回项式传递函数 的单位冲激响应y向量和时间状态历史记录x向量 4 5 3时域响应的其它函数1 求连续LTI系统的零输入响应函数initial 格式 y t x initial a b c d x0 功能 计算出连续时间LTI系统由于初始状态x0所引起的零输入响应y 其中x为状态记录 t为仿真所用的采样时间向量 例4 30P153 2 求离散系统的零输入响应函数dinitial 格式 y x n dinitial a b c d x0 功能 计算离散时间LTI系统由初始状态x0所引起的零输入响应y和状态响应响应x 取样点数由函数自动选取 n为仿真所用的点数 例4 31P154 3 求连续系统的单位阶跃响应函数step 格式 Y T step sys 功能 返回系统的单位阶跃响应Y和仿真所用的时间向量T 自动选择仿真的时间范围 其中sys可为系统传递函数 TF 零极增益模型 ZPK 或状态空间模型 SS 4 求离散系统的单位阶跃响应函数dstep 格式 y x dstep num den 功能 返回多项式传递函数G z num z den z 表示的系统单位阶跃响应 4 6线性时不变系统的频率响应 1 求模拟滤波器Ha s 的频率响应函数freqs 格式 H freqs B A W 功能 计算由向量W rad s 指定的频率点上模拟滤器系统函数Ha s 的频率响应Ha j 结果存于H向量中 向量B和A分别为模拟滤波器系统函数Ha s 的分子和分母多项式系数 例4 32 已知某模拟滤波器的系统函数 求该模拟滤波器的频率响应 MATLAB源程序如下 B 1 A 12 61313 41422 61311 W 0 0 1 2 pi 5 freqs B A W 图4 30模拟滤波器的频率响应 例4 33 已知某滤波器的系统函数为 求该滤波器的频率响应 MATLAB源程序为 B 10000000 1 A 1 freqz B A 该程序运行所绘出的幅频与相频性曲线如图4 31所示 图4 31滤波器幅度和相位曲线 2 求数字滤波器H z 的频率响应函数freqz 格式 H freqz B A W 功能 计算由向量W rad 指定的数字频率点上 通常指 0 范围的频率 数字滤波器H z 的频率响应H ejw 结果存于H向量中 向量B和A分别为数字滤波器系统H z 的分子和分母多项式系数 3 滤波函数filter格式 y filter B A x 功能 对向量x中的数据进行滤波处理 即差分方程求解 产生输出序列向量y B和A分别为数字滤波器系统函数H z 的分子和分母多项式系数向量 例4 34 设系统差分方程为 MATLAB源程序为 B 1 A 1 0 8 n 0 31 x 0 8 n y filter B A x subplot 2 1 1 stem x subplot 2 1 2 stem y 该程序运行所得结果如图4 32所示 求该系统对信号 的响应 图4 32系统对信号的响应 4 7傅里叶 Fourier 变换 4 7 1连续时间 连续频率 傅里叶变换 4 7 2连续时间 离散频率 傅里叶级数 正变换 逆变换 正变换 逆变换 4 7 3时间离散 连续频率 序列傅里叶变换 4 7 4离散时间 离散频率 离散傅里叶级数 4 7 5离散时间 离散频率 离散傅里叶变换 DFT 正变换 逆变换 正变换 逆变换 正变换 逆变换 1 一维快速正傅里叶变换函数fft格式 X fft x N 功能 采用FFT算法计算序列向量x的N点DFT变换 当N缺省时 fft函数自动按x的长度计算DFT 当N为2整数次幂时 fft按基 2算法计算 否则用混合算法 2 一维快速逆傅里叶变换函数ifft格式 x ifft X N 功能 采用FFT算法计算序列向量X的N点IDFT变换 例4 37 用快速傅里叶变换FFT计算下面两个序列的卷积 并测试直接卷积和快速卷积的时间 图4 35快速卷积框图 MATLAB程序 部分 线性卷积xn sin 0 4 1 15 对序列x n 赋值 M 15hn 0 9 1 20 对序列h n 赋值 N 20yn conv xn hn 直接调用函数conv计算卷积 园周卷积L pow2 nextpow2 M N 1 Xk fft xn L Hk fft hn L Yk Xk Hk yn ifft Yk L 图4 36x n h n 及其线性卷积波形 4 8IIR数字滤波器的设计方法 1 数字滤波器的频率响应函数 幅度响应 相位响应 图4 37理想低通 高通 带通 带阻数字滤波器幅度特性 2 滤波器的技术指标幅度响应指标 相位响应指标 图4 38数字低通滤波器的幅度特性 通带要求 阻带要求 通带最大衰减 阻带最小衰减 4 8 1冲激响应不变法 2 MATLAB信号处理工箱中的专用函数impinvar 格式 BZ AZ impinvar B A Fs 功能 把具有 B A 模拟滤波器传递函数模型转换成采样频率为Fs Hz 的数字滤波器的传递函数模型 BZ AZ 采样频率Fs的默认值为Fs 1 1 冲激响应不变法设计IIR数字滤波器的基本原理 例4 38 MATLAB源程序如下 num 1 模拟滤波器系统函数的分子den 1 sqrt 5 2 sqrt 2 1 模拟滤波器系统函数的分母 num1 den1 impinvar num den 求数字低通滤波器的系统函数程序的执行结果如下 num1 0 00000 09420 21580 0311den1 1 0000 2 00321 9982 0 76120 1069 MATLAB信号处理工具箱中的专用双线性变换函数bilinear 格式 numd dend bilinear num den Fs 功能 把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型 4 8 2双线性变换法 双线性变换利用频率变换关系 例4 39 MATLAB源程序如下 num 1 模拟滤波器系统函数的分子den 1 sqrt 3 sqrt 2 1 模拟滤波器系统函数的分母 num1 den1 bilinear num den 1 求数字滤波器的传递函数运算的结果如下 num1 0 05330 15990 15990 0533den1 1 0000 1 33820 9193 0 1546 4 8 3IIR数字滤波器的频率变换设计法 1 IIR数字滤波器的频率变换设计法的基本原理根据滤波器设计要求 设计模拟原型低通滤波器 然后进行频率变换 将其转换为相应的模拟滤波器 高通 带通等 最后利用冲激响应不变法或双线性变换法 将模拟滤波器数字化成相应的数字滤波器 图4 39IIR数字滤波器MATLAB设计步骤流程图 1 MATLAB的典型设计 利用在MATLAB设计IIR数字滤波器可分以下几步来实现 1 按一定规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标 2 根据转换后的技术指标使用滤波器阶数函数 确定滤波器的最小阶数N和截止频率Wc 3 利用最小阶数N产生模拟低通滤波原型 4 利用截止频率Wc把模拟低通滤波器原型转换成模拟低通 高通 带通或带阻滤波器 5 利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器 例4 40 设计一个数字信号处理系统 它的采样率为Fs 100Hz 希望在该系统中设计一个Butterworth型高通数字滤波器 使其通带中允许的最小衰减为0 5dB 阻带内的最小衰减为40dB 通带上限临界频率为30Hz 阻带下限临界频率为40Hz MATLAB源程序设计如下 把数字滤波器的频率特征转换成模拟滤波器的频率特征wp 30 2 pi ws 40 2 pi rp 0 5 rs 40 Fs 100 N Wc buttord wp ws rp rs s 选择滤波器的最小阶数 Z P K buttap N 创建Butterworth低通滤波器原型 A B C D zp2ss Z P K 零 极点增益模型转换为状态空间模型 AT BT CT DT lp2hp A B C D Wc 实现低通向高通的转变 num1 den1 ss2tf AT BT CT DT 状态空间模型转换为传递函数模型 运用双线性变换法把模拟滤波器转换成数字滤波器 num2 den2 bilinear num1 den1 100 H W freqz num2 den2 求频率响应plot W Fs 2 pi abs H grid 绘出频率响应曲线xlabel 频率 Hz ylabel 幅值 程序运行结果如图4 40所示 2 MATLAB的直接设计 图4 39IIR数字滤波器MATLAB设计步骤流程图 例4 42 试设计一个带阻IIR数字滤波器 其具体的要求是 通带的截止频率 wp1 650Hz wp2 850Hz 阻带的截止频率 ws1 700Hz ws2 800Hz 通带内的最大衰减为rp 0 1dB 阻带内的最小衰减为rs 50dB 采样频率为Fs 2000Hz MATLAB源程序设计如下 wp1 650 wp2 850 ws1 700 ws2 800 rp 0 1 rs 50 Fs 2000 wp wp1 wp2 Fs 2 ws ws1 ws2 Fs 2 利用Nyquist频率进行归一化 N wc ellipord wp ws rp rs z 求滤波器阶数 num den ellip N rp rs wc stop 求滤波器传递函数 H W freqz num den 绘出频率响应曲线plot W Fs 2 pi abs H grid xlabel 频率 Hz ylabel 幅值 该程序运行后的幅频响应曲线如图4 42所示 4 9FIR数字滤波器设计 格式 w boxcar M 功能 返回M点矩形窗序列 MATLAB信号处理工具箱中的窗函数法设计FIR数字滤波器的专用命令fir1 格式 B fir1 N wc 功能 设计一个具有线性相位的N阶 N点 的低通FIR数字滤波器 返回的向量B为滤波器的系数 单位冲激响应序列 其长度为N 1 4 9 1窗函数设计法窗函数设计的基本原理 h n w n hd n w n 为窗函数 hd n 理想数字滤波器的单位冲激响应 在MATLAB信号处理工具箱中为用户提供了Boxcar 矩形 Bartlet 巴特利特 Hanning 汉宁 等窗函数 这些窗函数的调用格式相同 FIR数字滤波器的单位冲激响应h n 满足偶 奇 对称h n h N n 1 或h n h N n 1 FIR数字滤波器具有线性相位 或 例4 44 用矩形窗设计线性相位FIR低通滤波器 该滤波器的通带截止频率wc pi 4 单位脉冲响h n 的长度M 21 并绘出h n 及其幅度响应特性曲线 MATLAB源程序为 M 21 wc pi 4 理想低通滤波器参数n 0 M 1 r M 1 2 nr n r eps n r 0 hdn sin wc nr pi nr