序列逆Z变换的Matlab实现

一 序列逆Z变换的Matlab实现 函数residuez 适合计算离散系统有理函数的留数和极点 可以用于求解序列的逆Z变换 函数residuez基本调用方式 r p c residuez b a 输入参数 b b0 b1 bM 为分子多项式的系数 a a0 a1 aN 为分母多项式的系数 这些多项式都按z的降幂排列输出参数 r是极点的留数 p是极点 c是无穷项多项式的系数项 仅当M N时存在 例 计算逆Z变换 例 计算的逆Z变换 解 有理分式X z 分子和分母多项式都按z的降幂排列 b 0 1 a 2 3 1 多项式的系数 r p c residuez b a 求留数 极点和系数项disp 留数 disp r 显示输出参数disp 极点 disp p disp 系数项 disp c 程序运行结果为留数 1 1极点 1 00000 5000系数项 X z 的部分分式形式为 逆Z变换为 二 DFT物理意义的Matlab实现 序列的N点DFT的物理意义 对X ej 在 0 2 上进行N点的等间隔取样 函数fft用于快速计算离散傅里叶变换 调用方式为 y fft x y fft x N y fft x 利用FFT算法计算序列x的离散傅里叶变换 当x为矩阵时 y为矩阵x每一列的FFT 当x长度为2的整数次幂时 函数fft采用基 2的FFT算法 否则采用混合基算法 y fft x N 采用N点FFT 当序列x长度小于N时 函数fft自动对序列尾部补零 构成N点数据 当x长度大于N时 函数fft自动截取序列前面N点数据进行FFT 函数ifft用于快速计算向量或矩阵的离散傅里叶逆变换 与函数fft的调用规则基本相同 调用方式为 y ifft x y ifft x N 例 利用FFT实现线性卷积 例 利用FFT实现线性卷积 已知序列x n R4 n 求 1 用conv函数求x n 与x n 的线性卷积y n 并绘出图形 2 用FFT求x n 与x n 的4点循环卷积y1 n 并绘出图形 3 用FFT求x n 与x n 的8点循环卷积y2 n 并将结果与 1 比较 说明线性卷积与循环卷积之间的关系 解程序如下 N1 4 N2 8 n1 0 1 N1 1 n2 0 1 N2 1 x 1 1 1 1 构造序列x n x1 1 1 1 1 0 0 0 0 在序列x n 后补4个零figure 1 subplot 2 2 1 stem n1 x gridon title 序列x n y1 conv x x y1为x n 与x n 的线性卷积subplot 2 2 2 stem 0 1 length y1 1 y1 gridon title x n 与x n 线性卷积 X2 fft x 计算x n 与x n 的4点循环卷积Y2 X2 X2 y2 ifft Y2 subplot 2 2 3 stem n1 y2 gridon title x n 与x n 的4点循环卷积 X3 fft x1 计算x n 与x n 的8点循环卷积Y3 X3 X3 y3 ifft Y3 subplot 2 2 4 stem n2 y3 gridon title x n 与x n 的8点循环卷积 程序运行结果图 三 频域取样定理的Matlab实现 例 设x n 0 7 n u n 在单位圆上以M 5和M 20 对其Z变换取样 研究时域信号受M变化的影响 1 对x n 进行Z变换 2 对X z 进行等角取样 取样点数为M 求X k 3 对X k 进行IDFT变化 得到M点序列 请比较几个序列 并作分析 解x n 0 7 n u n 的Z变换为 程序清单 n 0 19 x 0 7 n na 0 4 za exp j 2 pi na 5 在z平面的单位圆上对其进行5点的等角距取样Xa za za 0 7 xa abs ifft Xa nb 0 19 zb exp j 2 nb pi 20 在z平面的单位圆上对其进行20点的等角距取样Xb zb zb 0 7 xb abs ifft Xb figure 1 subplot 2 2 1 画出原始时域信号stem n x title 时域信号x n subplot 2 2 2 xa xa xa xa xa stem n xa title 5点取样恢复的序列 subplot 2 2 3 stem n xb title 20点取样恢复的序列 程序运行结果 四 用FFT进行谱分析的Matlab实现 设模拟信号 以t 0 01n n 0 N 1 进行取样 试用fft函数对其做频谱分析 N分别为 1 N 45 2 N 50 3 N 55 2 N 60 程序清单如下 计算N 45的FFT并绘出其幅频曲线N 45 n 0 N 1 t 0 01 n q n 2 pi N x 2 sin 4 pi t 5 cos 8 pi t y fft x N figure 1 subplot 2 2 1 plot q abs y title FFTN 45 程序清单 计算N 50的FFT并绘出其幅频曲线N 50 n 0 N 1 t 0 01 n q n 2 pi N x 2 sin 4 pi t 5 cos 8 pi t y fft x N figure 1 subplot 2 2 2 plot q abs y title FFTN 50 计算N 55的FFT并绘出其幅频曲线N 55 n 0 N 1 t 0 01 n q n 2 pi N x 2 sin 4 pi t 5 cos 8 pi t y fft x N figure 1 subplot 2 2 3 plot q abs y title FFTN 55 计算N 60的FFT并绘出其幅频曲线N 60 n 0 N 1 t 0 01 n q n 2 pi N x 2 sin 4 pi t 5 cos 8 pi t y fft x N figure 1 subplot 2 2 4 plot q abs y title FFTN 60 程序运行结果 从图中可以看出 这几种情况下均有较好的精度 程序运行结果分析 分析 由t 0 01n进行取样可得 采样频率fs 100Hz 而连续信号的最高模拟角频率为 8 由 2 f可得 最高频率为8 2 4Hz 因此 满足采样定理的要求 采样序列为 即 为周期序列 周期N 50 将程序中plot改为stem函数 则可以更清楚地看出频谱 修改程序运行结果 6 5IIR数字滤波器的Matlab仿真实现 IIR数字滤波器设计模拟滤波器到数字滤波器的转换 五 IIR数字滤波器设计 设数字滤波器系统函数为模拟滤波器的系统函数为函数butter和cheby1可以确定Butterworth和ChebyshevI型滤波器的系统函数 函数butter的调用格式 函数butter的调用格式为 b a butter n Wc 设计数字Butterworth滤波器 b a butter n Wc ftype 设计模拟Butterworth滤波器其中 n为滤波器阶数 Wc为截止频率 函数cheby1的调用格式 函数cheby1的调用格式为 b a cheby1 n Rp Wc 设计数字Chebyshev滤波器 b a cheby1 n Rp Wc ftype 设计模拟Chebyshev滤波器其中 n为滤波器阶数 Rp为通带内的纹波系数 Wc为截止频率 例 设计butterworth低通滤波器 设计一模拟butterworth低通滤波器 通带截止频率300Hz 通带最大衰减2dB 阻带截止频率800Hz 阻带最小衰减30dB 解滤波器的阶数和截止频率可由式和确定 程序段为 Wp 2 pi 300 Ws 2 pi 800 Rp 2 Rs 30 N ceil log10 10 0 1 Rs 1 10 0 1 Rp 1 2 log10 Ws Wp Wc Wp 10 Rp 10 1 1 2 N b a butter N Wc s freqs b a 程序运行结果 运行程序 得到N 4 Wc 2 0157e 003 幅频特性和相频特性如图7 16所示 模拟滤波器到数字滤波器的转换 设模拟滤波器系统函数为数字滤波器的系统函数为从模拟滤波器到数字滤波器的转换有两种方法 即脉冲响应不变法和双线性变换法 b a分别为模拟滤波器的系统函数的系数 bz az分别为数字滤波器的系统函数的系数 脉冲响应不变法 脉冲响应不变法 用代换Ha s 中的 s sk 即可得到H z 从而将模拟滤波器转换为数字滤波器格式 可用函数impinvar实现 调用格式为 bz az impinvar b a fs 其中 fs为取样频率 双线性变换法 双线性变换法 用代换Ha s 中的s即可得到H z 从而将模拟滤波器转换为数字滤波器格式 可用函数bilinear实现 调用格式为 zd pd kd bilinear z p k fs 其中 z p k和zd pd kd分别为s域和z域系统函数的零点 极点和增益 例 模拟滤波器转换数字滤波器 利用impinvar将一模拟低通滤波器变换成数字滤波器 取样频率为10Hz 程序段为 b a butter 4 3 s bz az impinvar b a 10 程序运行结果为bz 1 0e 006 0 00000 13240 51920 12730az 1 0000 3 92165 7679 3 77090 9246 六 FIR数字滤波器的Matlab仿真实现 窗函数法设计FIR滤波器FIR滤波器的优化设计 窗函数法设计FIR滤波器 窗函数法是通过对理想滤波器的单位取样响应加窗来逼近理想滤波器的 函数fir1用于设计标准的低通 带通 高通和带阻滤波器 函数fir1的调用格式 函数fir1的调用格式为 b fir1 n Wc ftype Windows 其中 n为滤波器阶数 Wc为截止频率ftype决定滤波器类型 ftype high 设计高通FIR滤波器 ftype stop 设计带阻FIR滤波器 Windows指定窗函数类型 默认为Hamming窗 可选Hanning Hamming Blackman triangle bartlett和boxcar窗 每种窗都可以由Matlab的相应函数生成 例 设计butterworth低通滤波器 设计一个15阶FIR低通滤波器 截止频率为0 2 b fir1 15 0 2 freqz b 1 512 函数freqz b a N 用于计算由a和b构成的数字滤波器的频率响应 并用图形方式分别表示其幅度响应和相位响应 程序运行结果 程序运行结果如图所示 FIR滤波器的优化设计 前面介绍了FIR滤波器的优化设计方法 通过迭代的方法求解FIR滤波器 过程十分复杂 在Matlab中 可以调用函数remez实现滤波器的设计 b