Matlab滤波器设计：Z变换与Z逆变换原理及Matlab实现代码

用户1143655

发布于 2023-02-23 11:52:45

2.9K0

发布于 2023-02-23 11:52:45

文章被收录于专栏：全栈之殇

Z变换在离散时间信号与系统中的地位相当于拉普拉斯变换在连续时间信号与系统中的地位。它可以求解常系数差分方程，进而估算一个线性时不变系统的响应及线性滤波器的设计。

一、Z变换的数学表述方法及Matlab实现代码

1、Z变换的数学表述方法

对于连续时间信号

x(t)

，其拉普拉斯变换为

X(s)

，即

x(t) \leftrightarrow X(s)

。经过离散化处理，每隔一个采样周期

T_s

取得一个采样点，进而根据拉普拉斯变换延迟定理，可以得到

x(t-nT_s) \leftrightarrow e^{-snT_s}X(s)

。

因此，对于离散时间序列

x(k) = {x(0), x(1), x(2), ...}

，采用拉普拉斯变换的延迟定理可以得到

X_A(s) = X_A(0) + X_A(1)e^{-sT_s} + x_A(2)e^{-2sT_s} + ...

其中，

x_A(n)

与

X_A(s)

表示采样后的时间信号与对应的拉普拉斯变换，进而

X_A(s)

可以重新表示为：

X_A(s) = \sum_{n=-\infty}^{\infty}x_A(nT_s)e^{-nsT_s} \tag{1}

其中，

x_A(nT_s)

：为模拟信号

x(t)

在各个

nT_s

时刻的采样值。

为了进一步简化数学表示，得到如下的数学变换，

z=e^{sT_s} \ \ 或者 \ \ z^{-1}=e^{-sT_s} \tag{2}

进而得到Z变换的数学表达式：

X(z) = \sum_{-\infty}^\infty x(n)z^{-n} \tag{3}

其中具有如下关系：

X(z)|_{z=e^{sT}} = X(e^{sT}) = X_A(S) \tag{4}

此时就从连续系统的拉普拉斯变换转换为离散系统的Z变换，因此，Z变换就是离散系统的拉普拉斯变换。

其中

(3)

式是双边Z变换，但是对于因果关系的时间序列，如下所示，单边Z变换是用单边求和的形式来描述：

X(z) = \sum_{n=0}^\infty x(x)z^{-n} \tag{5}

而对于有限时间序列，如果

x(n)

序列长度为

，Z变换的表达式可以用有限项求和形式来表示：

X(z) = \sum_{n=0}^{N-1} x(x)z^{-n} \tag{6}

S平面与Z平面的映射关系如下所示：在S平面上通常使用直角坐标系统，因此复变量s可以表示为

s=\sigma+{\rm j} \Omega

；而在Z平面上通常使用极坐标系统，即

z=re^{{\rm j} \omega}

，将它们代入

(2)

得到，并将采样周期简写为

：

re^{{\rm j}\omega} = e^{(\sigma + {\rm j}\Omega)T} = e^{\sigma T}e^{{\rm j}\Omega T} \tag{7}

即，

\begin{cases} r=e^{\sigma T} \newline \omega = \Omega T \end{cases} \tag{8}

上式表明，z的模值

与

的实部

\sigma

相关；而z的幅角

\omega

与

的虚部

\Omega

相关，其具体的关系如下所示：

（1）对于

与

\sigma

的关系，

表示z的模值，即z变量到原点的距离；而

\sigma

表示

变量的实部，由

(8)

式可得如下所示两者对应关系：

\begin{cases} \sigma = 0 \rightarrow r=1 \newline \sigma < 0 \rightarrow r < 1 \newline \sigma >0 \rightarrow r > 1 \end{cases} \tag{9}

如下图所示，上式表明S平面的虚轴（即

\sigma=0

）时，映射到Z平面上半径为

（

r=1

）的圆，即单位圆；而S的左半平面（

\sigma < 0

）映射到Z平面上的单位圆内区域（

r<1

）；S平面右半平面（

\sigma > 0

）映射到Z平面上是单位圆外（

r>1

）。

（2）S平面的虚变量

\Omega

与Z平面的幅角

\omega

之间的关系，显然

\omega=\Omega T

是一种线性关系。当

\Omega=0

时，

\omega=0

，表明S平面的实轴映射到Z平面上是正实轴。根据公式

(9)

可得，S平面上的原点

s=0

正好映射到Z平面的

z=1

处。

\Omega \ne 0

时，

\Omega

与

\omega

的映射关系就稍微复杂。当

\Omega

由

-\pi/T

增长到

\pi/T

时，即从

-\Omega_s/2

增长到

\Omega_s/2

时（

\Omega_s = 2\pi /T

为抽样频率），

\omega

从

-\pi

增长到

\pi

，对应Z平面上的幅角旋转了一周。也就是说S平面的虚轴仅从

-\pi/T

到

\pi/T

的区域，就已经映射到整个Z平面。

如果

\sigma < 0

，则对应于Z平面单位圆内旋转一周；

如果

\sigma > 0

，则对应于Z平面单位圆外旋转一周；

如果

\sigma = 0

，则对应于Z平面单位圆上旋转一周。

进而，当S的虚部

\Omega

由

\pi/T

增长到

3\pi/T

（即

\Omega

由

\Omega_s/2

增长到

3\Omega_s/2

）时，

\omega

由

\pi

增长到

3\pi

，由于

z=re^{{\rm j}\omega}

是

\omega

的周期函数，所以此时仍然映射到Z平面上的同样位置，只不过是在旋转一周的基础上再旋转一周，进而在Z平面上重叠一次。此时，这种多值函数的映射关系可以直观看作将S平面分割成一条条宽度为

\Omega_s

的横条。这些横条中的每一条都同样地映射到整个Z平面上且相互重叠在一起，可以表示为如下所示的数学公式：

H(z)|_{z=e^{sT}}=\frac{1}{T}\sum_{k=-\infty}^\infty H_a(s-{\rm j}k\Omega_s)=\frac{1}{T}\sum_{k=-\infty}^\infty H_a\bigg(s - {\rm j}k\frac{2\pi}{T} \bigg) \tag{10}

2、Z变换的Matlab实现代码

在Matlab的符号运算中Z变换的函数为ztrans。

（1）当

|a|<1

时，对

x(n)=a^nu(n)

进行Z变换，Matlab实现代码如下所示：

syms a n  % 声明符号变量

x = a^n;
X = ztrans(x);
X

运行代码，执行结果为X=-z/(a-z)，其数学公式如下所示：

X(z)=\frac{1}{1-az^{-1}}

（2）将正弦波

x(n)=sin(an)u(n)

进行Z变换，其Matlab命令如下所示：

syms a n

x = sin(a * n);
X = ztrans(x);
X

运行代码，执行结果为X=(z*(sin a))/(z^2-2*cos a*z + 1)，其数学公式如下所示：

X(z)=(sin a)z/(z^2-2z\cos a + 1)

二、Z变换的收敛域

通常，序列的Z变换

X(z)=\sum_{n=-\infty}^\infty x(n)z^{-n}

并不一定对任何

值都收敛，在Z平面上满足级数收敛的区域成为收敛域（Region of Convergence, ROC）。根据级数的知识，级数一致收敛条件是绝对可积。也就是说，如果在Z平面上的某处（即某个z值点）级数一致收敛，则对此z值满足：

\sum_{n=-\infty}^{\infty} |x(n)z^{-n}| < \infty \tag{11}

即其各项模值的和必须有界。由此可以想到，Z平面上的收敛域总是环形，因此

(11)

式可以表示为：

\sum_{n=-\infty}^\infty |x(n)z^{-n}| = \sum_{n=-\infty}^\infty |x(n)| \cdot |z^{-n}| < \infty \tag{12}

因此，式

X(z)=\sum_{n=0}^\infty x(n)z^{-n}

幂级数的收敛域由满足

(12)

式的全部

值所组成。因此，如果某个

z=z_1

值是在ROC内，即全部由

|z|=|z_1|

确定的圆上的z值也一定在ROC内。如下图所示，结果收敛域一定由在Z平面内以原点为中心的圆环所组成。收敛域的外边界是一个圆（或者可能外向延伸到无穷大），而内边界也是一个圆（或者ROC向内可包括原点）。

通常，级数在Z平面上的收敛域范围可以表示为：

R_1 < |Z| < R_2 \tag{13}

上式表明收敛域是一个以

R_1

和

R_2

为半径的两个圆所围城的环带区域，其中

R_1

为内圆半径，

R_2

为外圆半径，同时

R_1

与

R_2

也称为收敛半径。

三、基本Z变换对

通常，简单信号的Z变换我们可以通过下表查到其变换结果及其收敛域ROC：

四、线性系统的Z变换

在Z平面上对数字线性系统进行建模与分析时，通常的方法是用

\delta

函数作为输入激励序列，通过系统输出（脉冲响应）来分析线性系统。根据这些研究，可以推导出任意输入信号的线性系统响应。通常，松弛型（初始条件为零）线性常系数系统或滤波器的输入-输出关系可由差分方程表示为：

\sum_{m=0}^N a_m y(n-m) = \sum_{m=0}^M b_m x(n-m) \tag{14}

其中，

y(n)

为第

次的输出采样值；

x(n)

为第

次的输入采样值。

设

y(n)

的Z变换为

Y(z)

，

x(n)

的Z变换为

X(z)

，则根据延迟定理可得

y(n-m) \leftrightarrow z^{-m}

以及

x(n-m) \leftrightarrow z^{-m}X(z)

，代入

(14)

式得到：

\bigg(\sum_{m=0}^N a_mz^{-m} \bigg)Y(z) = \bigg(\sum_{m=0}^M b_mz^{-m} \bigg)X(z) \tag{15}

其中

Y(z)

与

X(z)

的比值称为传递函数，记为

H(z)

。对于

(14)

式所描述的线性系统，其传递函数为：

H(z)=\frac{Y(z)}{X(z)}=\frac{\sum_{m=0}^M b_mz^{-m}}{\sum_{m=0}^N a_m z^{-m}} \tag{16}

(15)

式与

(16)

式描述了输入信号Z变换与输出信号Z变换之间的变换关系，该转换关系涵盖了系统的大量信息。

五、Z变换的特性

通常，大部分重要信号可由1.3部分表格中的一个或多个初等函数表示，然而这些原始信号可能需要经过适当的调整、修改与合并。结合下表所示的Z变换特性可以得到如何通过初等信号的变换与组合来构建复杂信号：

六、Z逆变换

1、Z逆变换的数学原理

已知函数

X(z)

及其收敛域，反过来求序列的变换称为Z逆变换，表示为：

x(n)=Z^{-1}[X(z)] \tag{17}

其数学表述方法通常表示为如下数学公式：

x(n)=\frac{1}{2\pi{\rm j}} \oint_c X(z)z^{n-1}{\rm d}z, \ \ \ \ c \in (R_{x1}, R_{x2}) \tag{18}

其中，如下图所示，

是包围

X(z)z^{n-1}

所有极点的逆时针闭合积分路线，通常选Z平面收敛域内以原点为中心的圆。

以下给出由Z变换定义表达式导出逆变换式的过程，已知

X(z)=\sum_{n=-\infty}^\infty x(n)z^{-n}

，对其两端分别乘以

z^{m-1}

，如上图所示然后沿围线

进行积分，得到：

\oint_c z^{m-1}X(z)dz = \sum_{n=-\infty}^{\infty}x(n)\oint_c z^{m-n-1}{\rm d}z \tag{19}

根据复变函数中的柯西定理：

\oint_c z^{k-1}{\rm d}z = \begin{cases} 2\pi{\rm j}, \ \ k=0 \newline 0, \ \ \ \ \ \ k \ne 0 \end{cases}

由于

(19)

式右边只有

m=n

项，其余均为

，于是

(19)

式可写为：

\oint_c X(z)z^{n-1}dz = 2\pi{\rm j}x(n)

最终可以得到Z逆变换的数学表达式

(18)

。

直接计算围线积分是比较麻烦的，实际上求Z逆变换时的方法包括：

围线积分法（留数法）；
部分分式展开法；
幂级数展开法（长除法）。

为了更好的理解如何使用Matlab现成的函数求Z逆变换，下面以部分分式展开法为例，介绍Z逆变换的求解过程：

在数字信号处理中，

X(z)

通常是

z^{-1}

的有理函数，通常可采用部分分式分解将其变换为简单因式的和。假设

X(z)

为

X(z)=\frac{b_0+b_1z^{-1}+...+b_Mz^{-M}}{1+a_1z^{-1}+...+a_Nz^{-N}} \tag{20}

如果

M \ge N

，上式可表示为：

X(z) = \underbrace{\frac{\overline{b}_0 + \overline{b}_1z^{-1} + ... + \overline{b}_Mz^{-M}}{1+a_1z^{-1}+...+a_Nz^{-N}}}_{\text{真有理分式}} + \underbrace{\sum_{k=0}^{M-N}C_kz^{-k}}_{\text{如果}M \ge N \text{，则为直接多项式}} \tag{21}

将(21)式中的真有理式部分进行部分分式展开，可得：

X_1(z)=\sum_{k=1}^N \frac{R_k}{1-p_kz^{-1}} \tag{22}

根据本文第三部分的基本Z变换对表格及Z变换的可叠加性，可以得到

(22)

式的Z逆变换：

x(n)=\sum_{k=1}^NR_kp_k^nu(n)+\sum_{k=0}^{M-N}C_k\delta(n-k) \tag{23}

2、Z逆变换的Matlab留数函数实现方法

在实际应用中我们不必手算，可以使用计算机代替。Matlab提供了极点留数计算的两个函数：

(1) residue：针对拉普拉斯算子

的极点留数计算函数,适用于连续系统；

(2) residuez：针对Z变换算子的，适用于离散系统。

residuez函数的调用格式为：

[r,p,C] = residuez(b,a)

其中，

b和a为按照z^{-1}升幂排列的多项式

(20)

的分子和分布系数向量；

p为分母的根向量，即

X(z)

的极点向量；

r为对应于分母根向量中各个根的留数向量；
C为无穷项多项式系数向量，仅在

M \ge N

时存在。

由residuez函数计算得到r、p、C向量后，就可以将多项式

(20)

分解成如下形式：

\frac{B(z)}{A(z)}=\frac{r(1)}{1-p(1)z^{-1}}+\frac{r(2)}{1-p(2)z^{-1}}+...+\frac{r(N)}{1-p(N)z^{-1}}+C(1)+C(2)z^{-1}+... \tag{21}

根据

(23)

式，对上式作逆变换，得到其时域信号表达式为：

y(n)=r(1)p(1)^nu(n)+...+r(N)p(N)^nu(n)+C(1)\delta(n)+C(2)\delta(n-1)+...

其中，

u(n)

为阶跃函数，

\delta(n)

为脉冲函数。

3、Z逆变换Matlab符号函数实现方法

🚀 除了上面的residuez函数，Matlab的符号运算中有Z逆变换的iztrans函数。

比如我们已经知道

X(z)=z/(z-0.5)

，则可以很容易地通过下面函数实现：

syms z n;
X = z/(z-0.5);
x = iztrans(X,z,n);

但是函数iztrans处理高阶

X(z)

的Z逆变换时，有时候很难解释Z逆变换符号运算的结果。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-12-15，如有侵权请联系 cloudcommunity@tencent.com 删除

matlab

编程算法

本文分享自人工智能技术栈微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

matlab

编程算法

登录后参与评论

0 条评论

热度