信号与系统仿真：离散时间系统仿真_（2）.离散时间系统的数学模型

离散时间系统的数学模型

在信号与系统仿真中，离散时间系统的数学模型是描述系统行为的基础。离散时间系统是指在离散时间点上处理信号的系统，通常通过采样连续时间信号得到离散时间信号。本节将详细介绍离散时间系统的几种常见数学模型，包括差分方程、状态空间表示和Z变换。

差分方程

差分方程是描述离散时间系统最直接的方式之一。差分方程通过表示系统的当前输出与过去输出和当前及过去输入之间的关系来描述系统的行为。一般形式的差分方程可以表示为：

$$y[n]=\sum _{k=0}^{N-1}{a}_{k}y[n-k]+\sum _{k=0}^{M-1}{b}_{k}x[n-k]$$

其中，$y[n]$ 是系统的输出信号，$x[n]$ 是系统的输入信号，$a_k$ 和 $b_k$ 是系统的系数，$N$ 和 $M$ 分别是输出和输入的最大延迟。

一阶差分方程

考虑一个一阶差分方程：

$$y[n]=ay[n-1]+bx[n]$$

这个方程描述了一个简单的离散时间系统，其中输出 $y[n]$ 依赖于前一个输出 $y[n-1]$ 和当前输入 $x[n]$。

例子：一阶差分方程的实现

假设我们有一个简单的离散时间系统，其差分方程为：

$$y[n]=0.5y[n-1]+1.0x[n]$$

我们可以通过Python代码来实现这个系统。首先，我们需要定义输入信号和初始条件，然后通过迭代计算输出信号。

import numpy as np
import matplotlib.pyplot as plt

# 定义系统参数
a = 0.5
b = 1.0

# 定义输入信号
x = np.zeros(100)
x[10:20] = 1  # 在第10到19个时间点上输入为1

# 初始化输出信号
y = np.zeros(100)
y[0] = 0  # 初始条件

# 计算输出信号
for n in range(1, 100):
    y[n] = a * y[n-1] + b * x[n]

# 绘制输入信号和输出信号
plt.figure(figsize=(10, 5))
plt.subplot(2, 1, 1)
plt.plot(x, label='输入信号 x[n]')
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(y, label='输出信号 y[n]')
plt.legend()
plt.show()

高阶差分方程

对于高阶差分方程，形式更为复杂，但原理相同。考虑一个二阶差分方程：

$$y[n]=a_{1}y[n-1]+a_{2}y[n-2]+b_{0}x[n]+b_{1}x[n-1]$$

这个方程描述了一个二阶离散时间系统，其中输出 $y[n]$ 依赖于前两个输出 $y[n-1]$ 和 $y[n-2]$，以及当前和前一个输入 $x[n]$ 和 $x[n-1]$。

例子：二阶差分方程的实现

假设我们有一个二阶差分方程系统：

$$y[n]=0.5y[n-1]+0.2y[n-2]+1.0x[n]+0.8x[n-1]$$

我们可以通过Python代码来实现这个系统。首先，定义输入信号和初始条件，然后通过迭代计算输出信号。

import numpy as np
import matplotlib.pyplot as plt

# 定义系统参数
a1 = 0.5
a2 = 0.2
b0 = 1.0
b1 = 0.8

# 定义输入信号
x = np.zeros(100)
x[10:20] = 1  # 在第10到19个时间点上输入为1

# 初始化输出信号
y = np.zeros(100)
y[0] = 0  # 初始条件
y[1] = 0  # 初始条件

# 计算输出信号
for n in range(2, 100):
    y[n] = a1 * y[n-1] + a2 * y[n-2] + b0 * x[n] + b1 * x[n-1]

# 绘制输入信号和输出信号
plt.figure(figsize=(10, 5))
plt.subplot(2, 1, 1)
plt.plot(x, label='输入信号 x[n]')
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(y, label='输出信号 y[n]')
plt.legend()
plt.show()

状态空间表示

状态空间表示是描述离散时间系统的一种更通用的方法。通过定义系统的状态变量，可以将系统的动态行为描述为状态变量的演化过程。状态空间模型通常包括状态方程和输出方程。

状态方程

状态方程描述了系统状态变量的演化过程：

$$\mathbf{x}[n+1]=\mathbf{Ax}[n]+\mathbf{Bu}[n]$$

其中，$\mathbf{x}[n]$ 是状态向量，$\mathbf{u}[n]$ 是输入向量，$\mathbf{A}$ 是状态矩阵，$\mathbf{B}$ 是输入矩阵。

输出方程

输出方程描述了系统输出与状态变量和输入之间的关系：

$$\mathbf{y}[n]=\mathbf{Cx}[n]+\mathbf{Du}[n]$$

其中，$\mathbf{y}[n]$ 是输出向量，$\mathbf{C}$ 是输出矩阵，$\mathbf{D}$ 是直接传输矩阵。

例子：状态空间模型的实现

假设我们有一个二阶离散时间系统，其状态空间表示如下：

$$\mathbf{x}[n+1]=\left[\begin{array}{cc}0.5& 0.2\\ 1& 0\end{array}\right]\mathbf{x}[n]+\left[\begin{array}{c}1.0\\ 0.8\end{array}\right]\mathbf{u}[n]$$

$$\mathbf{y}[n]=\left[\begin{array}{cc}1& 0\end{array}\right]\mathbf{x}[n]+0\mathbf{u}[n]$$

我们可以通过Python代码来实现这个系统。首先，定义输入信号和初始状态，然后通过迭代计算状态变量和输出信号。

import numpy as np
import matplotlib.pyplot as plt

# 定义系统参数
A = np.array([[0.5, 0.2], [1, 0]])
B = np.array([[1.0], [0.8]])
C = np.array([[1, 0]])
D = 0

# 定义输入信号
u = np.zeros(100)
u[10:20] = 1  # 在第10到19个时间点上输入为1

# 初始化状态向量
x = np.zeros((2, 100))
x[:, 0] = np.array([0, 0])  # 初始状态

# 初始化输出信号
y = np.zeros(100)

# 计算状态变量和输出信号
for n in range(99):
    x[:, n+1] = A @ x[:, n] + B * u[n]
    y[n] = C @ x[:, n] + D * u[n]

# 绘制输入信号和输出信号
plt.figure(figsize=(10, 5))
plt.subplot(2, 1, 1)
plt.plot(u, label='输入信号 u[n]')
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(y, label='输出信号 y[n]')
plt.legend()
plt.show()

Z变换

Z变换是分析离散时间系统的一种重要工具。Z变换将离散时间信号从时域变换到复频域，使得系统的分析和设计更为方便。Z变换的定义为：

$$X(z)=\sum _{n=-\infty }^{\infty }x[n]z^{-n}$$

其中，$X(z)$ 是 $x[n]$ 的Z变换，$z$ 是复变量。

传递函数

通过Z变换，可以得到系统的传递函数。传递函数表示了系统输出的Z变换与输入的Z变换之间的关系：

$$H(z)=\frac{Y(z)}{X(z)}$$

对于一个差分方程：

$$y[n]=ay[n-1]+bx[n]$$

其传递函数为：

$$H(z)=\frac{b}{1-a{z}^{-1}}$$

例子：离散时间系统的传递函数

假设我们有一个简单的离散时间系统，其差分方程为：

$$y[n]=0.5y[n-1]+1.0x[n]$$

我们可以通过Z变换得到其传递函数：

$$H(z)=\frac{1.0}{1-0.5z^{-1}}$$

然后，我们可以使用Python中的 scipy 库来分析这个系统的频率响应。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 定义系统参数
b = [1.0]
a = [1, -0.5]

# 创建离散时间系统
system = signal.dlti(b, a)

# 计算频率响应
w, h = signal.dfreqresp(system, w=np.logspace(-2, 2, 1000))

# 绘制频率响应
plt.figure(figsize=(10, 5))
plt.subplot(2, 1, 1)
plt.semilogx(w, 20 * np.log10(abs(h)), label='幅度响应')
plt.ylabel('幅度 [dB]')
plt.legend()
plt.subplot(2, 1, 2)
plt.semilogx(w, np.angle(h), label='相位响应')
plt.ylabel('相位 [rad]')
plt.xlabel('频率 [rad/sample]')
plt.legend()
plt.show()

逆Z变换

逆Z变换将系统的传递函数从复频域变换回时域。逆Z变换的定义为：

$$x[n]=\frac{1}{2\pi j}{\oint }_{C}X(z)z^{n-1}dz$$

其中，$C$ 是复平面中的闭合路径。

例子：逆Z变换的实现

假设我们有一个传递函数：

$$H(z)=\frac{1.0}{1-0.5z^{-1}}$$

我们可以通过逆Z变换得到其时域表示。在Python中，可以使用 scipy 库中的 signal.dimpulse 函数来计算系统的冲激响应。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 定义系统参数
b = [1.0]
a = [1, -0.5]

# 创建离散时间系统
system = signal.dlti(b, a)

# 计算冲激响应
t, y = signal.dimpulse(system, n=100)

# 绘制冲激响应
plt.figure(figsize=(10, 5))
plt.plot(t, y[0], label='冲激响应 y[n]')
plt.xlabel('时间 [n]')
plt.ylabel('幅度')
plt.legend()
plt.show()

系统的稳定性

系统的稳定性是离散时间系统分析中一个重要的概念。一个离散时间系统是稳定的，如果其冲激响应是绝对可和的，即：

$$\sum _{n=-\infty }^{\infty }|h[n]|<\infty$$

例子：系统稳定性的分析

假设我们有一个离散时间系统的传递函数：

$$H(z)=\frac{1.0}{1-0.5z^{-1}}$$

我们可以通过计算系统的极点来分析其稳定性。如果所有极点都在单位圆内，则系统是稳定的。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 定义系统参数
b = [1.0]
a = [1, -0.5]

# 创建离散时间系统
system = signal.dlti(b, a)

# 计算系统的极点
poles = signal.dlti(b, a).zeros

# 绘制极点
plt.figure(figsize=(6, 6))
plt.plot(np.real(poles), np.imag(poles), 'x', label='极点')
plt.axhline(0, color='black', linewidth=0.5)
plt.axvline(0, color='black', linewidth=0.5)
plt.legend()
plt.title('极点图')
plt.xlabel('实部')
plt.ylabel('虚部')
plt.grid(True)
plt.show()

系统的频率响应

系统的频率响应是描述系统在不同频率下输出与输入之间关系的重要指标。可以通过Z变换或离散傅里叶变换（DFT）来计算系统的频率响应。

例子：系统的频率响应分析

假设我们有一个离散时间系统的传递函数：

$$H(z)=\frac{1.0}{1-0.5z^{-1}}$$

我们可以通过 scipy 库中的 signal.dfreqresp 函数来计算系统的频率响应。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 定义系统参数
b = [1.0]
a = [1, -0.5]

# 创建离散时间系统
system = signal.dlti(b, a)

# 计算频率响应
w, h = signal.dfreqresp(system, w=np.logspace(-2, 2, 1000))

# 绘制频率响应
plt.figure(figsize=(10, 5))
plt.subplot(2, 1, 1)
plt.semilogx(w, 20 * np.log10(abs(h)), label='幅度响应')
plt.ylabel('幅度 [dB]')
plt.legend()
plt.subplot(2, 1, 2)
plt.semilogx(w, np.angle(h), label='相位响应')
plt.ylabel('相位 [rad]')
plt.xlabel('频率 [rad/sample]')
plt.legend()
plt.show()

离散时间系统的仿真

离散时间系统的仿真是在计算机上模拟系统的行为，以便验证和分析系统的性能。常用的仿真工具有MATLAB、Python等。

例子：使用Python进行离散时间系统的仿真

假设我们有一个离散时间系统的差分方程：

$$y[n]=0.5y[n-1]+1.0x[n]$$

我们可以通过Python代码来实现系统的仿真，包括计算输出信号和绘制结果。

import numpy as np
import matplotlib.pyplot as plt

# 定义系统参数
a = 0.5
b = 1.0

# 定义输入信号
x = np.zeros(100)
x[10:20] = 1  # 在第10到19个时间点上输入为1

# 初始化输出信号
y = np.zeros(100)
y[0] = 0  # 初始条件

# 计算输出信号
for n in range(1, 100):
    y[n] = a * y[n-1] + b * x[n]

# 绘制输入信号和输出信号
plt.figure(figsize=(10, 5))
plt.subplot(2, 1, 1)
plt.plot(x, label='输入信号 x[n]')
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(y, label='输出信号 y[n]')
plt.legend()
plt.show()

例子：使用MATLAB进行离散时间系统的仿真

假设我们有一个离散时间系统的差分方程：

$$y[n]=0.5y[n-1]+1.0x[n]$$

我们可以通过MATLAB代码来实现系统的仿真，包括计算输出信号和绘制结果。

% 定义系统参数
a = 0.5;
b = 1.0;

% 定义输入信号
x = zeros(1, 100);
x(11:20) = 1;  % 在第10到19个时间点上输入为1

% 初始化输出信号
y = zeros(1, 100);
y(1) = 0;  % 初始条件

% 计算输出信号
for n = 2:100
    y(n) = a * y(n-1) + b * x(n);
end

% 绘制输入信号和输出信号
figure;
subplot(2, 1, 1);
plot(x, 'DisplayName', '输入信号 x[n]');
legend;
subplot(2, 1, 2);
plot(y, 'DisplayName', '输出信号 y[n]');
legend;

系统的辨识

系统的辨识是指通过实验数据来确定系统的数学模型。常见的辨识方法包括最小二乘法、递推最小二乘法等。

最小二乘法

最小二乘法通过最小化输出误差的平方和来估计系统的参数。假设我们有一个线性离散时间系统：

$$y[n]=a_{1}y[n-1]+a_{2}y[n-2]+b_{0}x[n]+b_{1}x[n-1]$$

可以通过最小二乘法来估计参数 $a_1,a_2,b_0,b_1$。

例子：最小二乘法的实现

假设我们有一组实验数据，输入信号 $x[n]$ 和输出信号 $y[n]$，我们可以通过最小二乘法来估计系统的参数。

import numpy as np
import matplotlib.pyplot as plt

# 定义实验数据
x = np.random.randn(100)  # 随机生成的输入信号
a1_true = 0.5
a2_true = 0.2
b0_true = 1.0
b1_true = 0.8

# 生成输出信号（假设已知系统的参数）
y = np.zeros(100)
y[0] = 0  # 初始条件
y[1] = 0  # 初始条件
for n in range(2, 100):
    y[n] = a1_true * y[n-1] + a2_true * y[n-2] + b0_true * x[n] + b1_true * x[n-1]

# 构建数据矩阵
X = np.column_stack((y[:-2], y[:-1], x[:-2], x[:-1]))
Y = y[2:]

# 使用最小二乘法估计参数
params = np.linalg.lstsq(X, Y, rcond=None)[0]

# 提取估计的参数
a1_est = params[1]
a2_est = params[0]
b0_est = params[2]
b1_est = params[3]

print(f"估计的参数: a1 = {a1_est}, a2 = {a2_est}, b0 = {b0_est}, b1 = {b1_est}")

# 使用估计的参数重新生成输出信号
y_est = np.zeros(100)
y_est[0] = 0  # 初始条件
y_est[1] = 0  # 初始条件
for n in range(2, 100):
    y_est[n] = a1_est * y_est[n-1] + a2_est * y_est[n-2] + b0_est * x[n] + b1_est * x[n-1]

# 绘制实际输出信号和估计输出信号
plt.figure(figsize=(10, 5))
plt.plot(y, label='实际输出信号 y[n]')
plt.plot(y_est, label='估计输出信号 y_est[n]', linestyle='--')
plt.legend()
plt.xlabel('时间 [n]')
plt.ylabel('幅度')
plt.title('实际输出信号与估计输出信号的比较')
plt.grid(True)
plt.show()

递推最小二乘法

递推最小二乘法是一种在线辨识方法，可以在系统运行过程中不断更新参数估计。递推最小二乘法通过递推的方式逐步减小输出误差的平方和。

例子：递推最小二乘法的实现

假设我们有一组实验数据，输入信号 $x[n]$ 和输出信号 $y[n]$，我们可以通过递推最小二乘法来估计系统的参数。

import numpy as np
import matplotlib.pyplot as plt

# 定义实验数据
x = np.random.randn(100)  # 随机生成的输入信号
a1_true = 0.5
a2_true = 0.2
b0_true = 1.0
b1_true = 0.8

# 生成输出信号（假设已知系统的参数）
y = np.zeros(100)
y[0] = 0  # 初始条件
y[1] = 0  # 初始条件
for n in range(2, 100):
    y[n] = a1_true * y[n-1] + a2_true * y[n-2] + b0_true * x[n] + b1_true * x[n-1]

# 初始化参数估计
params = np.zeros(4)
P = np.eye(4) * 1000  # 初始化协方差矩阵
lambda_ = 0.98  # 遗忘因子

# 估计参数
for n in range(2, 100):
    phi = np.array([y[n-2], y[n-1], x[n-2], x[n-1]])
    K = (P @ phi) / (lambda_ + phi.T @ P @ phi)
    params = params + K * (y[n] - phi.T @ params)
    P = (P - K @ phi.T @ P) / lambda_

# 提取估计的参数
a1_est = params[1]
a2_est = params[0]
b0_est = params[2]
b1_est = params[3]

print(f"估计的参数: a1 = {a1_est}, a2 = {a2_est}, b0 = {b0_est}, b1 = {b1_est}")

# 使用估计的参数重新生成输出信号
y_est = np.zeros(100)
y_est[0] = 0  # 初始条件
y_est[1] = 0  # 初始条件
for n in range(2, 100):
    y_est[n] = a1_est * y_est[n-1] + a2_est * y_est[n-2] + b0_est * x[n] + b1_est * x[n-1]

# 绘制实际输出信号和估计输出信号
plt.figure(figsize=(10, 5))
plt.plot(y, label='实际输出信号 y[n]')
plt.plot(y_est, label='估计输出信号 y_est[n]', linestyle='--')
plt.legend()
plt.xlabel('时间 [n]')
plt.ylabel('幅度')
plt.title('实际输出信号与估计输出信号的比较')
plt.grid(True)
plt.show()

系统的控制

在信号与系统分析中，系统的控制是一个重要的应用领域。通过设计控制器，可以使系统在各种条件下保持稳定和响应性能。常见的控制器设计方法包括PID控制、状态反馈控制等。

PID控制

PID控制是一种常用的控制方法，通过比例、积分和微分三个部分来调节系统的输出。PID控制器的输出可以表示为：

$$u[n]=K_{p}e[n]+K_i\sum _{k=0}^{n}e[k]+K_d(e[n]-e[n-1])$$

其中，$e[n]=r[n]-y[n]$ 是误差信号，$r[n]$ 是参考信号，$K_p,K_i,K_d$ 是控制器的增益。

例子：PID控制器的实现

假设我们有一个离散时间系统，其差分方程为：

$$y[n]=0.5y[n-1]+1.0x[n]$$

我们可以通过设计一个PID控制器来调节系统的输出。

import numpy as np
import matplotlib.pyplot as plt

# 定义系统参数
a = 0.5
b = 1.0

# 定义输入信号
x = np.zeros(100)
x[10:20] = 1  # 在第10到19个时间点上输入为1

# 初始化输出信号
y = np.zeros(100)
y[0] = 0  # 初始条件

# 定义参考信号
r = np.ones(100) * 1.0  # 参考信号为1

# 定义PID控制器参数
Kp = 0.8
Ki = 0.1
Kd = 0.05

# 初始化PID控制器的变量
integral = 0
previous_error = 0

# 计算输出信号
for n in range(1, 100):
    error = r[n] - y[n-1]
    integral += error
    derivative = error - previous_error
    u = Kp * error + Ki * integral + Kd * derivative
    y[n] = a * y[n-1] + b * u
    previous_error = error

# 绘制参考信号、输入信号和输出信号
plt.figure(figsize=(10, 5))
plt.subplot(3, 1, 1)
plt.plot(r, label='参考信号 r[n]')
plt.legend()
plt.subplot(3, 1, 2)
plt.plot(x, label='输入信号 x[n]')
plt.legend()
plt.subplot(3, 1, 3)
plt.plot(y, label='输出信号 y[n]')
plt.legend()
plt.xlabel('时间 [n]')
plt.ylabel('幅度')
plt.title('PID控制器的输出信号')
plt.grid(True)
plt.show()

状态反馈控制

状态反馈控制是一种基于系统状态变量的控制方法。通过选择合适的反馈增益矩阵 $\mathbf{K}$，可以使系统的闭环极点在期望的位置。

例子：状态反馈控制器的实现

假设我们有一个二阶离散时间系统，其状态空间表示如下：

$$\mathbf{x}[n+1]=\left[\begin{array}{cc}0.5& 0.2\\ 1& 0\end{array}\right]\mathbf{x}[n]+\left[\begin{array}{c}1.0\\ 0.8\end{array}\right]\mathbf{u}[n]$$

$$\mathbf{y}[n]=\left[\begin{array}{cc}1& 0\end{array}\right]\mathbf{x}[n]+0\mathbf{u}[n]$$

我们可以通过设计一个状态反馈控制器来调节系统的输出。

import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg

# 定义系统参数
A = np.array([[0.5, 0.2], [1, 0]])
B = np.array([[1.0], [0.8]])
C = np.array([[1, 0]])
D = 0

# 定义参考信号
r = np.ones(100) * 1.0  # 参考信号为1

# 定义状态反馈控制器参数
K = np.array([1.0, 0.5])  # 反馈增益矩阵

# 初始化状态向量
x = np.zeros((2, 100))
x[:, 0] = np.array([0, 0])  # 初始状态

# 初始化输出信号
y = np.zeros(100)

# 计算状态变量和输出信号
for n in range(99):
    error = r[n] - y[n]
    u = -K @ x[:, n] + error
    x[:, n+1] = A @ x[:, n] + B * u
    y[n] = C @ x[:, n] + D * u

# 绘制参考信号、输入信号和输出信号
plt.figure(figsize=(10, 5))
plt.subplot(3, 1, 1)
plt.plot(r, label='参考信号 r[n]')
plt.legend()
plt.subplot(3, 1, 2)
plt.plot(x[0, :], label='状态 x1[n]')
plt.plot(x[1, :], label='状态 x2[n]')
plt.legend()
plt.subplot(3, 1, 3)
plt.plot(y, label='输出信号 y[n]')
plt.legend()
plt.xlabel('时间 [n]')
plt.ylabel('幅度')
plt.title('状态反馈控制器的输出信号')
plt.grid(True)
plt.show()

总结

离散时间系统的数学模型是描述系统行为的基础。差分方程、状态空间表示和Z变换是常用的数学模型，分别适用于不同的分析和设计需求。系统的稳定性分析和频率响应计算是确保系统性能的重要手段。此外，通过系统辨识可以确定系统的参数，而通过控制系统设计可以实现期望的系统行为。这些工具和方法在信号与系统仿真中发挥着重要作用，帮助工程师和研究人员更好地理解系统的动态特性。