《信号与系统》- 第七章系统函数

啊阿狸不会拉杆

发布于 2026-01-21 11:52:55

1200

7.1 系统函数与系统特性

系统函数是分析线性时不变系统的重要工具，它将系统的时域特性和频域特性联系起来，通过系统函数可以全面了解系统的行为。

一、系统函数的零点与极点

系统函数 H (s) 或 H (z) 的零点和极点分布决定了系统的基本特性。对于连续时间系统，系统函数定义为输出信号与输入信号的拉普拉斯变换之比；对于离散时间系统，则为 Z 变换之比。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from control import tf, impulse_response, step_response

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 案例：连续时间系统的零极点分析
# 定义系统函数分子和分母多项式系数
num = [1, 4, 3]  # 分子多项式 s^2 + 4s + 3
den = [1, 5, 6, 0]  # 分母多项式 s^3 + 5s^2 + 6s

# 创建传递函数对象
sys = tf(num, den)

# 计算系统的零点和极点
zeros = sys.zero()
poles = sys.pole()

# 打印零极点
print("系统零点：", zeros)
print("系统极点：", poles)

# 绘制零极点图
plt.figure(figsize=(8, 6))
# 修改 color 参数为有效的颜色值，这里使用 'b' 表示蓝色
plt.scatter(zeros.real, zeros.imag, color='b', marker='o', s=200, facecolors='none', edgecolors='b', label='零点')
# 对于 'x' 标记，不需要设置 edgecolors，因为它没有填充区域
plt.scatter(poles.real, poles.imag, color='r', marker='x', s=200, label='极点')
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)
plt.grid(True, alpha=0.3)
plt.title('系统零极点图')
plt.xlabel('实部')
plt.ylabel('虚部')
plt.legend()
plt.axis('equal')
plt.show()

# 离散时间系统零极点分析
num_z = [1, -1.5, 0.5]  # 分子多项式 z^2 - 1.5z + 0.5
den_z = [1, -0.5, 0.25]  # 分母多项式 z^2 - 0.5z + 0.25

# 创建离散时间传递函数
sys_z = tf(num_z, den_z, dt=1)

# 计算离散系统零极点
zeros_z = sys_z.zero()
poles_z = sys_z.pole()

print("离散系统零点：", zeros_z)
print("离散系统极点：", poles_z)

图 1：系统零极点图放置位置

二、系统函数与时域响应

系统函数的极点决定了系统时域响应的模态，不同的极点位置对应不同的响应特性（如指数衰减、等幅振荡等）。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from control import tf, impulse_response, step_response

# 设置中文字体显示
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 定义系统函数
# 连续时间系统
num = [1, 4, 3]  # 分子多项式 s^2 + 4s + 3
den = [1, 5, 6, 0]  # 分母多项式 s^3 + 5s^2 + 6s
sys = tf(num, den)

# 离散时间系统
num_z = [1, -1.5, 0.5]  # 分子多项式 z^2 - 1.5z + 0.5
den_z = [1, -0.5, 0.25]  # 分母多项式 z^2 - 0.5z + 0.25
sys_z = tf(num_z, den_z, dt=1)

# 案例：系统函数与时域响应
# 连续时间系统冲激响应
try:
    t = np.linspace(0, 10, 1000)
    t, y = impulse_response(sys, T=t)

    plt.figure(figsize=(10, 6))
    plt.plot(t, y[0])
    plt.title('系统冲激响应')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅值')
    plt.grid(True, alpha=0.3)
    plt.show()
except Exception as e:
    print(f"计算冲激响应时出错: {e}")
    # 备选方法：使用scipy.signal.impulse
    try:
        t, y = signal.impulse((num, den), T=t)
        plt.figure(figsize=(10, 6))
        plt.plot(t, y)
        plt.title('系统冲激响应 (备选方法)')
        plt.xlabel('时间 (s)')
        plt.ylabel('幅值')
        plt.grid(True, alpha=0.3)
        plt.show()
    except:
        print("无法计算冲激响应")

# 连续时间系统阶跃响应
try:
    t, y_step = step_response(sys, T=t)

    plt.figure(figsize=(10, 6))
    plt.plot(t, y_step[0])
    plt.title('系统阶跃响应')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅值')
    plt.grid(True, alpha=0.3)
    plt.show()
except Exception as e:
    print(f"计算阶跃响应时出错: {e}")
    # 备选方法：使用scipy.signal.step
    try:
        t, y_step = signal.step((num, den), T=t)
        plt.figure(figsize=(10, 6))
        plt.plot(t, y_step)
        plt.title('系统阶跃响应 (备选方法)')
        plt.xlabel('时间 (s)')
        plt.ylabel('幅值')
        plt.grid(True, alpha=0.3)
        plt.show()
    except:
        print("无法计算阶跃响应")

# 离散时间系统单位脉冲响应
try:
    t_z = np.arange(0, 20)
    y_impulse = signal.lfilter(num_z, den_z, [1] + [0] * (len(t_z) - 1))

    plt.figure(figsize=(10, 6))
    plt.stem(t_z, y_impulse)
    plt.title('离散系统单位脉冲响应')
    plt.xlabel('时间 (n)')
    plt.ylabel('幅值')
    plt.grid(True, alpha=0.3)
    plt.show()
except Exception as e:
    print(f"计算离散系统脉冲响应时出错: {e}")

# 离散时间系统单位阶跃响应
try:
    t_z = np.arange(0, 20)
    y_step_z = signal.lfilter(num_z, den_z, np.ones_like(t_z))

    plt.figure(figsize=(10, 6))
    plt.stem(t_z, y_step_z)
    plt.title('离散系统单位阶跃响应')
    plt.xlabel('时间 (n)')
    plt.ylabel('幅值')
    plt.grid(True, alpha=0.3)
    plt.show()
except Exception as e:
    print(f"计算离散系统阶跃响应时出错: {e}")

图 2：系统冲激响应图放置位置

图 3：系统阶跃响应图放置位置

三、系统函数与频域响应

系统函数在 jω 轴上的取值即为频率响应，包括幅频响应和相频响应，反映了系统对不同频率信号的处理特性。

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

# 设置中文字体显示
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC", "Arial Unicode MS"]
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 定义系统函数
# 连续时间系统
num = [1, 4, 3]  # 分子多项式 s^2 + 4s + 3
den = [1, 5, 6, 0]  # 分母多项式 s^3 + 5s^2 + 6s

# 离散时间系统
num_z = [1, -1.5, 0.5]  # 分子多项式 z^2 - 1.5z + 0.5
den_z = [1, -0.5, 0.25]  # 分母多项式 z^2 - 0.5z + 0.25

# 案例：系统函数与频域响应
# 连续时间系统频域响应 - Bode图
try:
    w, mag, phase = signal.bode((num, den))

    plt.figure(figsize=(12, 8))
    plt.subplot(2, 1, 1)
    plt.semilogx(w, mag)  # 幅频响应
    plt.title('连续系统频域响应')
    plt.ylabel('幅值 (dB)')
    plt.grid(True, which='both', alpha=0.3)
    plt.xticks(rotation=45)  # 旋转x轴标签，避免重叠

    plt.subplot(2, 1, 2)
    plt.semilogx(w, phase)  # 相频响应
    plt.xlabel('角频率 (rad/s)')
    plt.ylabel('相位 (度)')
    plt.grid(True, which='both', alpha=0.3)
    plt.xticks(rotation=45)

    plt.tight_layout()
    plt.show()
except Exception as e:
    print(f"计算连续系统Bode图时出错: {e}")

# 离散时间系统频域响应
try:
    # 计算8000个频率点，提高分辨率
    w_z, h_z = signal.freqz(num_z, den_z, worN=8000)
    mag_z = np.abs(h_z)
    phase_z = np.angle(h_z, deg=True)

    # 处理幅值为零的情况，避免对数计算警告
    eps = 1e-10  # 一个很小的数，避免对零取对数
    mag_z_db = 20 * np.log10(np.maximum(mag_z, eps))

    plt.figure(figsize=(12, 8))
    plt.subplot(2, 1, 1)
    plt.plot(w_z / np.pi, mag_z_db)  # 幅频响应(dB)，使用处理后的mag_z_db
    plt.title('离散系统频域响应')
    plt.ylabel('幅值 (dB)')
    plt.grid(True, alpha=0.3)

    # 设置y轴下限，避免显示过小的值
    plt.ylim(bottom=-100)

    plt.subplot(2, 1, 2)
    plt.plot(w_z / np.pi, phase_z)  # 相频响应(度)
    plt.xlabel('归一化频率 (π为单位)')
    plt.ylabel('相位 (度)')
    plt.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()
except Exception as e:
    print(f"计算离散系统频域响应时出错: {e}")

# 频域响应的替代计算方法（使用control库）
try:
    from control import tf, bode_plot

    # 连续系统
    sys = tf(num, den)
    plt.figure(figsize=(12, 8))
    bode_plot(sys, dB=True, deg=True, omega_limits=(0.01, 100))
    plt.suptitle('连续系统频域响应 (control库)')
    plt.tight_layout(rect=[0, 0, 1, 0.96])  # 为suptitle留出空间
    plt.show()
except Exception as e:
    print(f"使用control库计算频域响应时出错: {e}")

图 4：连续系统频域响应图放置位置

图 5：离散系统频域响应图放置位置

7.1 思维导图

图 6：7.1 节思维导图

7.2 系统的因果性与稳定性

因果性和稳定性是系统的重要特性，决定了系统在实际应用中的可行性。

一、系统的因果性

因果系统的响应只取决于当前和过去的输入，因果性在频域和时域都有明确的条件。

# 导入必要的库
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from control import tf

# 设置中文字体显示
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 定义系统函数
# 连续时间系统
num = [1, 4, 3]  # 分子多项式 s^2 + 4s + 3
den = [1, 5, 6, 0]  # 分母多项式 s^3 + 5s^2 + 6s

# 创建连续时间系统
sys = tf(num, den)

# 计算系统的极点
poles = sys.pole()  # 确保 poles 变量在此处被定义

# 案例：因果性判断
# 连续时间系统因果性判断
# 因果系统要求极点全部在左半平面
print("连续系统因果性判断：")
if all(poles.real < 0):
    print("系统是因果稳定的")
else:
    print("系统不是因果稳定的")

# 离散时间系统
num_z = [1, -1.5, 0.5]  # 分子多项式 z^2 - 1.5z + 0.5
den_z = [1, -0.5, 0.25]  # 分母多项式 z^2 - 0.5z + 0.25

# 创建离散时间系统
sys_z = tf(num_z, den_z, dt=1)

# 计算离散系统极点
poles_z = sys_z.pole()  # 确保 poles_z 变量在此处被定义

# 离散时间系统因果性判断
# 因果离散系统要求极点全部在单位圆内
print("\n离散系统因果性判断：")
if all(np.abs(poles_z) < 1):
    print("离散系统是因果稳定的")
else:
    print("离散系统不是因果稳定的")

# 非因果系统示例
num_noncausal = [1, 0, -1]
den_noncausal = [1, 0.5, 0.25]
sys_noncausal = tf(num_noncausal, den_noncausal, dt=1)
poles_noncausal = sys_noncausal.pole()

print("\n非因果系统极点：", poles_noncausal)

# 正确判断非因果系统的因果性和稳定性
# 1. 先判断因果性（极点是否都在单位圆内）
is_causal = all(np.abs(poles_noncausal) < 1)

# 2. 判断稳定性（收敛域是否包含单位圆）
# 对于因果系统，稳定性等价于因果性
# 对于非因果系统，需要额外分析收敛域
if is_causal:
    print("系统是因果的")
    if is_causal:  # 对于因果系统，稳定性由因果性决定
        print("系统是稳定的")
    else:
        print("系统是不稳定的")
else:
    print("系统是非因果的")
    # 非因果系统的稳定性需要额外分析
    # 这里简化为：如果所有极点模都大于1，则系统不稳定
    if all(np.abs(poles_noncausal) > 1):
        print("系统是不稳定的")
    else:
        # 部分极点在单位圆内，部分在单位圆外
        # 此时稳定性取决于收敛域的选择
        print("系统稳定性取决于收敛域选择")

二、系统的稳定性

稳定系统对有界输入产生有界输出 (BIBO 稳定)，连续系统和离散系统的稳定性条件与极点位置密切相关。

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

# 设置中文字体
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams['axes.unicode_minus'] = False

# 连续系统
num = [1, 4, 3]
den = [1, 5, 6, 0]
sys = tf(num, den)
poles = sys.pole()  # 使用pole()方法获取极点

# 离散系统
num_z = [1, -1.5, 0.5]
den_z = [1, -0.5, 0.25]
sys_z = tf(num_z, den_z, dt=1)
poles_z = sys_z.pole()  # 使用pole()方法获取极点

# 案例：稳定性分析
# 连续系统稳定性判断
print("连续系统稳定性判断：")
if all(poles.real < 0):
    print("系统是稳定的")
else:
    print("系统是不稳定的")

# 离散系统稳定性判断
print("\n离散系统稳定性判断：")
if all(np.abs(poles_z) < 1):
    print("离散系统是稳定的")
else:
    print("离散系统是不稳定的")

# 不稳定系统示例
num_unstable = [1, 2]
den_unstable = [1, -1.5, 0.5]
sys_unstable = tf(num_unstable, den_unstable, dt=1)
poles_unstable = sys_unstable.pole()  # 统一使用pole()方法

print("\n不稳定系统极点：", poles_unstable)
print("不稳定系统稳定性判断：")
if all(np.abs(poles_unstable) < 1):
    print("系统是稳定的")
else:
    print("系统是不稳定的")

# 绘制不稳定系统的响应
t_unstable = np.arange(0, 20)
y_unstable = signal.lfilter(num_unstable, den_unstable, [1] + [0]*(len(t_unstable)-1))

plt.figure(figsize=(10, 6))
plt.stem(t_unstable, y_unstable)
plt.title('不稳定系统的单位脉冲响应')
plt.xlabel('时间 (n)')
plt.ylabel('幅值')
plt.grid(True, alpha=0.3)
plt.show()

图 7：不稳定系统脉冲响应图放置位置

7.2 流程图

图 8：7.2 节流程图

7.3 信号流图

信号流图是一种用节点和支路表示系统信号流动和运算关系的图形表示方法，便于系统分析和结构设计。

一、信号流图

信号流图由节点和支路组成，节点表示信号变量，支路表示信号的传输关系。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from control import tf
from sympy import symbols  # 导入sympy库中的symbols函数

# 设置中文字体
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams['axes.unicode_minus'] = False

# 定义符号变量s
s = symbols('s')

# 连续系统
num = [1, 4, 3]
den = [1, 5, 6, 0]
sys = tf(num, den)
poles = sys.pole()

# 离散系统
num_z = [1, -1.5, 0.5]
den_z = [1, -0.5, 0.25]
sys_z = tf(num_z, den_z, dt=1)
poles_z = sys_z.pole()

# 案例：稳定性分析
# 连续系统稳定性判断
print("连续系统稳定性判断：")
if all(poles.real < 0):
    print("系统是稳定的")
else:
    print("系统是不稳定的")

# 离散系统稳定性判断
print("\n离散系统稳定性判断：")
if all(np.abs(poles_z) < 1):
    print("离散系统是稳定的")
else:
    print("离散系统是不稳定的")

# 不稳定系统示例
num_unstable = [1, 2]
den_unstable = [1, -1.5, 0.5]
sys_unstable = tf(num_unstable, den_unstable, dt=1)
poles_unstable = sys_unstable.pole()

print("\n不稳定系统极点：", poles_unstable)
print("不稳定系统稳定性判断：")
if all(np.abs(poles_unstable) < 1):
    print("系统是稳定的")
else:
    print("系统是不稳定的")

# 绘制不稳定系统的响应
t_unstable = np.arange(0, 20)
y_unstable = signal.lfilter(num_unstable, den_unstable, [1] + [0]*(len(t_unstable)-1))

plt.figure(figsize=(10, 6))
plt.stem(t_unstable, y_unstable)
plt.title('不稳定系统的单位脉冲响应')
plt.xlabel('时间 (n)')
plt.ylabel('幅值')
plt.grid(True, alpha=0.3)
plt.show()

# 案例：信号流图表示与分析
# 以二阶系统为例：H(s) = (s+2)/(s^2+3s+2)
# 绘制信号流图的Mermaid表示
print("二阶系统信号流图Mermaid表示：")
print("```mermaid")
print("flowchart TD")
print("    A[X(s)] -->|1| W")
print("    W -->|s| Y")
print("    Y -->|s| Z")
print("    Z -->|-3| W")
print("    Z -->|-2| W")
print("    W -->|1| V")
print("    Y -->|1| V")
print("    V -->|2| U")
print("    U -->|1| Y_out[Y(s)]")
print("```")

# 信号流图到系统函数的转换
# 通过Mason公式计算系统函数
print("\n通过Mason公式计算系统函数：")

# 前向通路增益
forward_paths = [
    {'gain': 1*s*1, 'non_touching_loops': []},  # 路径1: X->W->Y->V->Y_out
    {'gain': 1*s*s*1*2*1, 'non_touching_loops': []}  # 路径2: X->W->Y->Z->W->V->U->Y_out
]

# 单独回路增益
loops = [
    {'gain': s*s*(-3), 'name': 'L1'},  # 回路1: Y->Z->W->Y
    {'gain': s*s*(-2), 'name': 'L2'}   # 回路2: Y->Z->W->Y(另一条)
]

# 计算特征式Δ
delta = 1
# 减去所有单独回路增益的和
loop_sum = sum(loop['gain'] for loop in loops)
delta -= loop_sum

# 计算两两不接触回路的乘积和
non_touching_pairs = []  # 假设没有不接触的回路
for pair in non_touching_pairs:
    loop1, loop2 = pair
    delta += loop1['gain'] * loop2['gain']

# 假设没有三个及以上不接触的回路

print(f"特征式 Δ = {delta}")

# 计算每个前向通路的余因子
for i, path in enumerate(forward_paths):
    path_delta = 1
    # 减去与该路径不接触的回路增益
    for loop in path['non_touching_loops']:
        path_delta -= loop['gain']
    # 加上与该路径不接触的回路对的乘积
    for pair in non_touching_pairs:
        if all(loop in path['non_touching_loops'] for loop in pair):
            path_delta += pair[0]['gain'] * pair[1]['gain']
    forward_paths[i]['delta'] = path_delta
    print(f"前向通路 {i+1} 的余因子 Δ_{i+1} = {path_delta}")

# 计算总系统函数
H = sum(path['gain'] * path['delta'] for path in forward_paths) / delta
print(f"系统函数 H(s) = {H}")

# 验证：与已知的系统函数比较
print("\n验证：")
print(f"已知系统函数: H(s) = (s+2)/(s^2+3s+2)")
print(f"Mason公式计算结果: H(s) = {H}")

二、梅森公式

梅森公式用于从信号流图直接计算系统函数，避免复杂的代数运算。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from control import tf
from sympy import symbols

# 设置中文字体
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC", "Arial Unicode MS"]
plt.rcParams['axes.unicode_minus'] = False

# 案例：梅森公式应用
# 考虑一个包含前向通路和回路的信号流图
# 前向通路增益
P1 = 5  # 前向通路1增益
P2 = 3  # 前向通路2增益

# 单独回路增益
L1 = -2  # 回路1增益
L2 = -1  # 回路2增益
L3 = -0.5  # 回路3增益

# 互不接触回路组合
L1L2 = L1 * L2  # 回路1和回路2互不接触
L1L3 = L1 * L3  # 回路1和回路3互不接触
L2L3 = L2 * L3  # 回路2和回路3互不接触

# 三个回路都互不接触的情况（如果存在）
L1L2L3 = L1 * L2 * L3  # 假设三个回路都互不接触

# 特征式Δ
Delta = 1 - (L1 + L2 + L3) + (L1L2 + L1L3 + L2L3) - L1L2L3

# 前向通路P1的余因子Δ1（假设P1与所有回路接触）
Delta1 = 1

# 前向通路P2的余因子Δ2（假设P2与回路3不接触）
Delta2 = 1 - L3

# 系统函数H = (P1*Delta1 + P2*Delta2)/Delta
H = (P1 * Delta1 + P2 * Delta2) / Delta

print("梅森公式计算系统函数案例：")
print(f"特征式Δ = {Delta:.4f}")
print(f"系统函数H = {H:.4f}")

# 验证：将信号流图转换为系统函数
# 修正：分子和分母应该是多项式系数列表
# 对于SISO系统，分子和分母需要至少两个元素，第一个元素是最高阶系数
num_mason = [0, P1*Delta1 + P2*Delta2]  # 分子为常数项，对应s的0次幂
den_mason = [Delta, 0]  # 分母为Delta，对应s的0次幂

try:
    sys_mason = tf(num_mason, den_mason)
    print("通过梅森公式得到的系统函数：", sys_mason)
except Exception as e:
    print(f"创建系统函数时出错: {e}")
    print("修正后的分子:", num_mason)
    print("修正后的分母:", den_mason)

7.3 信号流图示例

图 9：7.3 节信号流图示例

7.4 系统的结构

系统可以用不同的结构实现，不同结构在运算效率、存储需求和对有限字长效应的敏感性等方面有不同特点。

一、直接实现

直接按照系统函数的分子和分母多项式系数直接构建系统结构。

import numpy as np
import matplotlib.pyplot as plt
from control import tf
from control.matlab import step

# 设置中文字体
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams['axes.unicode_minus'] = False

# 案例：系统的直接实现
# 以三阶系统为例：H(s) = (s^2+4s+3)/(s^3+5s^2+6s)
# 直接实现结构的模拟

# 定义系统传递函数
num = [1, 4, 3]
den = [1, 5, 6, 0]
sys_direct = tf(num, den)

# 生成输入信号
t = np.linspace(0, 10, 1000)

# 计算直接实现的系统响应（使用step函数的正确方式）
try:
    t, y_direct = step(sys_direct, T=t)

    plt.figure(figsize=(10, 6))
    plt.plot(t, y_direct)  # 修正：直接使用y_direct，不需要索引[0]
    plt.title('直接实现的系统阶跃响应')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅值')
    plt.grid(True, alpha=0.3)
    plt.show()
except Exception as e:
    print(f"计算系统响应时出错: {e}")

# 直接实现的结构示意图（Mermaid表示）
print("\n直接实现结构示意图Mermaid表示：")
print("```mermaid")
print("flowchart TD")
print("    A[输入x(t)] --> B[加法器1]")
print("    B --> C1[积分器1]")
print("    C1 --> D1[积分器2]")
print("    D1 --> E1[积分器3]")
print("    E1 --> F[输出y(t)]")
print("    C1 -.->|增益-5| B")  # 修正：使用正确的Mermaid连线语法
print("    D1 -.->|增益-6| B")  # 修正：使用正确的Mermaid连线语法
print("    D1 -.->|增益0| B")  # 修正：从D1引出，而不是E1
print("    A -->|增益1| B")  # 修正：添加箭头标签语法
print("    C1 -->|增益1| C2[增益1]")
print("    D1 -->|增益1| D2[增益4]")
print("    E1 -->|增益1| E2[增益3]")
print("    C2 --> D2 --> E2 --> F")  # 合并连线
print("```")

二、级联和并联实现

将系统分解为多个子系统的级联或并联组合，便于分析和设计。

import numpy as np
import matplotlib.pyplot as plt
from control import tf, series, parallel
from control.matlab import step
from scipy import signal

# 设置中文字体
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams['axes.unicode_minus'] = False

# 案例：系统的级联和并联实现
# 级联实现
# 将系统分解为一阶或二阶子系统的级联
num1 = [1, 1]  # 分解为(s+1)
den1 = [1, 2]  # (s+2)
num2 = [1, 3]  # (s+3)
den2 = [1, 3, 0]  # s(s+3)

try:
    sys_cascade = series(tf(num1, den1), tf(num2, den2))
    print("级联实现的系统函数：", sys_cascade)
except Exception as e:
    print(f"创建级联系统时出错: {e}")

# 并联实现
# 将系统分解为部分分式之和
num_par = [1, 4, 3]
den_par = [1, 5, 6, 0]

# 部分分式展开
try:
    r, p, k = signal.residue(num_par, den_par)
    print("\n部分分式展开系数：", r)
    print("极点：", p)
    print("余项：", k)

    # 构建并联子系统
    sys_par1 = tf([r[0]], [1, -p[0]])
    sys_par2 = tf([r[1]], [1, -p[1]])

    # 处理常数项（如果有余项）
    if len(k) > 0:
        sys_par3 = tf([k[0]], [1])  # 常数项作为增益
    else:
        sys_par3 = tf([0], [1])  # 如果没有余项，创建零系统

    sys_parallel = parallel(sys_par1, sys_par2, sys_par3)
    print("并联实现的系统函数：", sys_parallel)
except Exception as e:
    print(f"创建并联系统时出错: {e}")

# 比较不同实现的响应
try:
    # 使用np.linspace生成时间点，然后验证并调整
    t = np.linspace(0, 10, 1001)  # 使用1001确保包含端点且等间距

    # 验证时间数组是否等间距
    dt = t[1] - t[0]
    if not np.allclose(np.diff(t), dt):
        raise ValueError("时间数组不是等间距的")

    # 计算直接实现的响应（如果之前已定义）
    try:
        t, y_direct = step(sys_direct, T=t)
    except NameError:
        # 如果sys_direct未定义，重新定义
        num_direct = [1, 4, 3]
        den_direct = [1, 5, 6, 0]
        sys_direct = tf(num_direct, den_direct)
        t, y_direct = step(sys_direct, T=t)

    # 计算级联和并联实现的响应
    t, y_cascade = step(sys_cascade, T=t)
    t, y_parallel = step(sys_parallel, T=t)

    plt.figure(figsize=(12, 8))
    plt.subplot(2, 1, 1)
    plt.plot(t, y_direct, 'b-', label='直接实现')
    plt.plot(t, y_cascade, 'r--', label='级联实现')
    plt.title('不同实现方式的阶跃响应比较')
    plt.ylabel('幅值')
    plt.legend()
    plt.grid(True, alpha=0.3)

    plt.subplot(2, 1, 2)
    plt.plot(t, y_direct, 'b-', label='直接实现')
    plt.plot(t, y_parallel, 'g-.', label='并联实现')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅值')
    plt.legend()
    plt.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()
except Exception as e:
    print(f"绘制响应比较图时出错: {e}")
    # 备选方案：不指定时间数组，让step函数自动生成
    try:
        print("尝试备选方案：使用自动生成的时间数组...")
        y_direct = step(sys_direct)[1]
        y_cascade = step(sys_cascade)[1]
        y_parallel = step(sys_parallel)[1]

        # 获取自动生成的时间数组
        t_auto = np.linspace(0, 10, len(y_direct))

        plt.figure(figsize=(12, 8))
        plt.subplot(2, 1, 1)
        plt.plot(t_auto, y_direct, 'b-', label='直接实现')
        plt.plot(t_auto, y_cascade, 'r--', label='级联实现')
        plt.title('不同实现方式的阶跃响应比较')
        plt.ylabel('幅值')
        plt.legend()
        plt.grid(True, alpha=0.3)

        plt.subplot(2, 1, 2)
        plt.plot(t_auto, y_direct, 'b-', label='直接实现')
        plt.plot(t_auto, y_parallel, 'g-.', label='并联实现')
        plt.xlabel('时间 (s)')
        plt.ylabel('幅值')
        plt.legend()
        plt.grid(True, alpha=0.3)

        plt.tight_layout()
        plt.show()
    except Exception as e2:
        print(f"备选方案也失败了: {e2}")