《信号与系统》- 第六章 离散系统的 z 域分析

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, Sum, oo, simplify, KroneckerDelta

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义符号变量
n, z = symbols('n z')

# 单位脉冲序列 δ[n] 的符号定义
delta_expr = KroneckerDelta(n, 0)

# 计算 z 变换
X_z = Sum(delta_expr * z**(-n), (n, -oo, oo))
X_z_simplify = simplify(X_z.doit())

print(f"单位脉冲序列 δ[n] 的 z 变换: {X_z_simplify}")

# 数值计算部分 - 绘制单位脉冲序列
n_vals = np.arange(-5, 6)
x_n = np.zeros_like(n_vals, dtype=int)
x_n[n_vals == 0] = 1  # 在n=0处值为1

plt.figure(figsize=(10, 4))
plt.subplot(121)
# 使用新的stem绘制方式
markerline, stemlines, baseline = plt.stem(n_vals, x_n)
plt.setp(markerline, markersize=8, markeredgewidth=2, color='blue')
plt.setp(stemlines, linewidth=1.5, color='blue')
plt.title('单位脉冲序列 δ[n]')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.grid(True)

# 添加z变换结果可视化
plt.subplot(122)
# 修复转义序列警告 - 使用原始字符串
plt.text(0.1, 0.5, r"$\mathcal{Z}\{\delta[n]\} = 1$",
         fontsize=16, verticalalignment='center')
plt.title('z变换结果')
plt.axis('off')
plt.tight_layout()
plt.savefig('delta_z_transform.png', dpi=300)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, Sum, oo, simplify

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

# 创建符号变量
n, z = symbols('n z')

# 右边序列 x[n] = a^n u[n], |a| < 1
a = 0.5
n_vals = np.arange(0, 11)
x_n = [a**n for n in n_vals]

# 计算 z 变换及收敛域
# 使用 SymPy 的 Sum 函数表示无限求和
X_z = Sum((a**n) * z**(-n), (n, 0, oo))
X_z_simplify = simplify(X_z.doit())  # doit() 执行求和计算

print(f"右边序列 a^n u[n] 的 z 变换: {X_z_simplify}")
print(f"收敛域: |z| > {a}")

# 绘制序列
plt.figure(figsize=(8, 6))
plt.stem(n_vals, x_n)
plt.title(f'右边序列 a^n u[n] (a={a})')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.grid(True)
plt.tight_layout()
plt.savefig('right_sided_sequence.png', dpi=300)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, Sum, oo, simplify, Heaviside

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

# 创建符号变量
n, z = symbols('n z', integer=True)  # 指定 n 为整数

# 定义两个序列（使用 SymPy 表达式）
def x1(n):
    return 1 if n == 0 else 0  # δ[n]

def x2(n):
    return (0.5**n) * Heaviside(n, 1)  # (0.5)^n u[n]，Heaviside(n, 1) 表示 u[n]

# 计算线性组合的 z 变换
a, b = 2, 3
n_vals = np.arange(0, 10)
x_combined = [a*x1(n) + b*x2(n) for n in n_vals]

# 使用 SymPy 的 Sum 函数进行符号求和
X1_z = Sum(x1(n)*z**(-n), (n, 0, oo)).doit()
X2_z = Sum(x2(n)*z**(-n), (n, 0, oo)).doit()
X_combined_z = a*X1_z + b*X2_z

# 简化结果
X1_z_simplified = simplify(X1_z)
X2_z_simplified = simplify(X2_z)
X_combined_z_simplified = simplify(X_combined_z)

print(f"δ[n] 的 z 变换: {X1_z_simplified}")
print(f"(0.5)^n u[n] 的 z 变换: {X2_z_simplified}")
print(f"线性组合 2δ[n] + 3(0.5)^n u[n] 的 z 变换: {X_combined_z_simplified}")

# 绘制序列及其线性组合
plt.figure(figsize=(15, 5))

plt.subplot(131)
plt.stem(n_vals, [x1(n) for n in n_vals])
plt.title('序列 x1[n] = δ[n]')
plt.xlabel('n')
plt.ylabel('x1[n]')
plt.grid(True)

plt.subplot(132)
plt.stem(n_vals, [x2(n) for n in n_vals])
plt.title('序列 x2[n] = (0.5)^n u[n]')
plt.xlabel('n')
plt.ylabel('x2[n]')
plt.grid(True)

plt.subplot(133)
plt.stem(n_vals, x_combined)
plt.title('线性组合 2x1[n] + 3x2[n]')
plt.xlabel('n')
plt.ylabel('组合序列')
plt.grid(True)

plt.tight_layout()
plt.savefig('linear_combination_z_transform.png', dpi=300)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, Sum, oo, simplify, Heaviside

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

# 创建符号变量
n, z = symbols('n z', integer=True)

# 原序列 x[n] = u[n]
def u(n):
    return Heaviside(n, 1)  # 使用Heaviside函数表示单位阶跃函数

# 右移k=2的序列 x[n-2]
k = 2
n_vals = np.arange(-5, 10)
x_n = [u(n) for n in n_vals]
x_shifted = [u(n-k) for n in n_vals]

# 计算z变换（使用符号求和）
X_z = Sum(u(n)*z**(-n), (n, 0, oo)).doit()
X_shifted_z = z**(-k) * X_z

# 简化表达式
X_z_simplified = simplify(X_z)
X_shifted_z_simplified = simplify(X_shifted_z)

print(f"原序列 u[n] 的 z 变换: {X_z_simplified}")
print(f"右移 {k} 位后的 z 变换: {X_shifted_z_simplified}")

# 绘制序列
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.stem(n_vals, x_n)
plt.title('原序列 u[n]')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.grid(True)

plt.subplot(1, 2, 2)
plt.stem(n_vals, x_shifted)
plt.title(f'右移 {k} 位序列 u[n-{k}]')
plt.xlabel('n')
plt.ylabel('x[n-{k}]')
plt.grid(True)

plt.tight_layout()
plt.savefig('shifted_sequence_z_transform.png', dpi=300)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, series, oo, simplify

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

# 创建符号变量
z = symbols('z')

# 定义 X(z)
X_z_expr = 1 / (1 - 0.5 * z**(-1))

# 在 z→∞ 处进行幂级数展开（使用 dir='+' 表示从正方向趋近于无穷）
X_series = series(X_z_expr, z, oo, n=10, dir='+')

# 提取系数（逆 z 变换结果）
x_n = []
for i in range(10):
    # 注意：在 z→∞ 展开时，系数对应 z^i 而非 z^(-i)
    coeff = X_series.coeff(z, i)
    x_n.append(coeff)

print("逆 z 变换结果 x[n]:", x_n)

# 绘制序列
n_vals = np.arange(0, 10)
plt.figure(figsize=(8, 6))
plt.stem(n_vals, x_n)
plt.title('逆 z 变换结果 x[n]')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.grid(True)
plt.tight_layout()
plt.savefig('inverse_z_transform_power_series.png', dpi=300)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, apart, simplify, Sum, Heaviside, lambdify

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

# 创建符号变量
z, n = symbols('z n', integer=True)

# 定义 X(z)/z
X_z_over_z = (z + 2) / ((z - 1) * (z - 2))

# 部分分式展开
X_partial_over_z = apart(X_z_over_z, z)
print(f"X(z)/z 的部分分式展开: {X_partial_over_z}")

# 乘以 z 得到 X(z) 的部分分式
X_partial = simplify(X_partial_over_z * z)
print(f"X(z) 的部分分式展开: {X_partial}")


# 手动计算逆 z 变换（查表法）
# X(z) = 3*z/(z-2) - 2*z/(z-1)
# 查表得: Z^-1{z/(z-a)} = a^n u[n]
def inverse_z_transform(X_partial_expr):
    # 分解为各个分式
    terms = []
    if '+' in str(X_partial_expr):
        terms = X_partial_expr.as_ordered_terms()
    else:
        terms = [X_partial_expr]

    x_n = 0
    for term in terms:
        # 提取系数和分式部分
        coeff = term.as_coeff_mul()[0]
        fraction = term / coeff

        # 识别分式形式 z/(z-a)
        if fraction == z / (z - 2):
            x_n += coeff * 2 ** n * Heaviside(n, 1)
        elif fraction == z / (z - 1):
            x_n += coeff * 1 ** n * Heaviside(n, 1)  # 1^n = 1

    return x_n


# 计算逆 z 变换
x_n_expr = inverse_z_transform(X_partial)
x_n_simplified = simplify(x_n_expr)
print(f"逆 z 变换结果: x[n] = {x_n_simplified}")

# 计算前10项数值
x_n_vals = [x_n_simplified.subs(n, k) for k in range(10)]
print("\n部分分式展开法逆 z 变换结果:")
for k, val in enumerate(x_n_vals):
    print(f"x[{k}] = {val}")

# 绘制序列
n_vals = np.arange(0, 10)
plt.figure(figsize=(8, 6))
plt.stem(n_vals, x_n_vals)
plt.title('逆 z 变换结果 x[n]')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.grid(True)
plt.tight_layout()
plt.savefig('inverse_z_transform_partial_fraction.png', dpi=300)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, solve, simplify, Heaviside

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

# 定义符号
n, z, Y_z = symbols('n z Y_z')

# 对差分方程两边取 z 变换
# y[n] - 0.5 y[n-1] = x[n]
# 左边 z 变换：Y(z) - 0.5 z^{-1} Y(z) + 0.5 y[-1] = Y(z)(1 - 0.5 z^{-1})
# 右边 z 变换：X(z) = z/(z-1) (u[n]的 z 变换)
X_z = z / (z - 1)
eq = Y_z * (1 - 0.5 * z**(-1)) - X_z
Y_z_sol = solve(eq, Y_z)[0]

# 简化 Y(z) 的表达式
Y_z_simplified = simplify(Y_z_sol)
print(f"Y(z) 的简化表达式: {Y_z_simplified}")

# 手动计算逆 z 变换
# Y(z) = z^2/((z-1)*(z-0.5))
# 使用部分分式展开和查表法
# Y(z)/z = z/((z-1)*(z-0.5)) = 2/(z-1) - 1/(z-0.5)
# 所以 Y(z) = 2*z/(z-1) - z/(z-0.5)
# 查表得: Z^-1{z/(z-a)} = a^n u[n]
def manual_inverse_z_transform(Y_z_expr):
    # 手动分解 Y(z) 为部分分式
    # 这里我们已经知道分解结果
    # Y(z) = 2*z/(z-1) - z/(z-0.5)
    y_n = 2 * (1**n) * Heaviside(n, 1) - (0.5**n) * Heaviside(n, 1)
    return y_n

# 计算逆 z 变换得到 y[n]
y_n_expr = manual_inverse_z_transform(Y_z_sol)
y_n_simplified = simplify(y_n_expr)
print(f"差分方程的解 y[n] = {y_n_simplified}")

# 计算前10项
y_n_vals = [y_n_simplified.subs(n, k) for k in range(10)]
print("\n差分方程的数值解:")
for k, val in enumerate(y_n_vals):
    print(f"y[{k}] = {val}")

# 绘制解序列
n_vals = np.arange(0, 10)
plt.figure(figsize=(8, 6))
plt.stem(n_vals, y_n_vals)
plt.title('差分方程的解 y[n]')
plt.xlabel('n')
plt.ylabel('y[n]')
plt.grid(True)
plt.tight_layout()
plt.savefig('difference_equation_solution.png', dpi=300)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

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

# 系统函数 H(z) = 1/(1 - 0.5z^-1)
omega = np.linspace(-np.pi, np.pi, 1000)  # 频率范围 [-π, π]
H = 1 / (1 - 0.5 * np.exp(-1j * omega))  # H(e^jω)

# 计算幅度和相位响应
magnitude = np.abs(H)
phase = np.angle(H)

# 找出特定频率点的响应
w_points = [-np.pi, -np.pi/2, 0, np.pi/2, np.pi]
H_points = 1 / (1 - 0.5 * np.exp(-1j * np.array(w_points)))
mag_points = np.abs(H_points)
phase_points = np.angle(H_points)

# 绘制幅度响应和相位响应
plt.figure(figsize=(16, 6))

# 幅度响应图
plt.subplot(1, 2, 1)
plt.plot(omega, magnitude, 'b-', linewidth=2)
plt.title('系统幅度响应 |H(e^jω)|')
plt.xlabel('角频率 ω (rad/sample)')
plt.ylabel('幅度')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi],
           [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$\pi/2$', r'$\pi$'])

# 标记特定频率点
plt.scatter(w_points, mag_points, color='red', s=50, zorder=5)
for w, m in zip(w_points, mag_points):
    plt.annotate(f'{m:.3f}', (w, m), textcoords="offset points",
                 xytext=(0,10), ha='center')

# 相位响应图
plt.subplot(1, 2, 2)
plt.plot(omega, phase, 'g-', linewidth=2)
plt.title('系统相位响应 ∠H(e^jω)')
plt.xlabel('角频率 ω (rad/sample)')
plt.ylabel('相位 (rad)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi],
           [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$\pi/2$', r'$\pi$'])
plt.yticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi],
           [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$\pi/2$', r'$\pi$'])

# 标记特定频率点
plt.scatter(w_points, phase_points, color='red', s=50, zorder=5)
for w, p in zip(w_points, phase_points):
    plt.annotate(f'{p:.3f}', (w, p), textcoords="offset points",
                 xytext=(0,10), ha='center')

plt.tight_layout()
plt.savefig('frequency_response.png', dpi=300, bbox_inches='tight')
plt.show()

# 打印特定频率点的响应值
print("特定频率点的频率响应:")
for w, mag, ph in zip(w_points, mag_points, phase_points):
    print(f"ω = {w/ np.pi:.2f}π rad: 幅度 = {mag:.4f}, 相位 = {ph:.4f} rad")

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.font_manager import FontProperties


def create_z_domain_mindmap():
    # 设置支持中文的字体
    plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS', 'Microsoft YaHei', 'WenQuanYi Micro Hei']
    plt.rcParams['axes.unicode_minus'] = False

    fig = plt.figure(figsize=(18, 12))
    ax = fig.add_subplot(111)
    ax.set_xlim(0, 12)
    ax.set_ylim(0, 12)
    ax.axis('off')

    # 使用标准字符替代特殊Unicode字符
    # 中心主题
    plt.text(6, 10, '离散系统的z域分析',
             fontsize=22, ha='center', va='center', fontweight='bold',
             bbox=dict(boxstyle="round,pad=0.5", fc="lightblue", ec="steelblue", lw=3))

    # 主分支
    main_branches = [
        (3, 8, '6.1 z变换'),
        (6, 8, '6.2 性质'),
        (9, 8, '6.3 逆z变换'),
        (12, 8, '6.4 系统分析')
    ]

    # 子分支 - 使用标准字符
    sub_branches = {
        '6.1 z变换': [
            (1.5, 6, '定义\n- 双边: X(z)=sum(x[n]z^{-n})\n- 单边: X(z)=sum_{n=0}^{∞}x[n]z^{-n}'),
            (2.0, 4, '收敛域\n- 右边序列: |z|>R\n- 左边序列: |z|<R\n- 双边序列: R1<|z|<R2'),
            (3.5, 2, '常用序列\n- δ[n] -> 1\n- u[n] -> z/(z-1)\n- a^n u[n] -> z/(z-a)')
        ],
        '6.2 性质': [
            (5.0, 6, '线性\na*x[n] + b*y[n] -> a*X(z) + b*Y(z)'),
            (6.0, 4, '移位特性\n- x[n-k] -> z^{-k}X(z)\n- x[n+k] -> z^k[X(z)-sum_{m=0}^{k-1}x[m]z^{-m}]'),
            (7.0, 2, '其他性质\n- 尺度: a^n x[n] -> X(z/a)\n- 卷积: x[n]*h[n] -> X(z)H(z)')
        ],
        '6.3 逆z变换': [
            (8.5, 6, '幂级数法\n- 右边序列: z^{-1}降幂\n- 左边序列: z升幂'),
            (9.0, 4, '部分分式\n- 一阶极点: X(z)/z = sum A_i/(z-p_i)\n- 高阶极点: 含重根项'),
            (10.0, 2, '留数法\nx[n] = sum Res[X(z)z^{n-1}]')
        ],
        '6.4 系统分析': [
            (11.0, 6, '差分方程\nY(z) = H(z)X(z)\n零状态/零输入响应'),
            (11.5, 4, '系统函数\nH(z) = Y(z)/X(z)\nH(z) = Z{h[n]}'),
            (12.0, 2, '稳定性\n- 因果系统: 极点全在单位圆内\n- 充要条件: sum |h[n]| < ∞'),
            (12.5, 0, '其他\n- 频率响应\n- 实现结构\n- s-z关系')
        ]
    }

    # 绘制主分支
    for x, y, text in main_branches:
        plt.plot([6, x], [10, y], 'k-', lw=2)
        circle = plt.Circle((x, y), 0.3, fc='lightgreen', ec='darkgreen', lw=2)
        ax.add_patch(circle)
        plt.text(x, y, text, fontsize=16, ha='center', va='center', fontweight='bold')

    # 绘制子分支
    for branch, points in sub_branches.items():
        branch_x, branch_y = next((x, y) for x, y, t in main_branches if t == branch)
        for x, y, text in points:
            plt.plot([branch_x, x], [branch_y, y], 'k--', lw=1.5)
            plt.text(x, y, text, fontsize=12, ha='center', va='center',
                     bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="gray", alpha=0.8))

    # 添加标题
    plt.title('离散系统的z域分析思维导图', fontsize=24, pad=20, y=0.95)

    # 添加单位圆示意图
    circle_ax = fig.add_axes([0.1, 0.1, 0.2, 0.2])
    circle_ax.set_aspect('equal')
    theta = np.linspace(0, 2 * np.pi, 100)
    circle_ax.plot(np.cos(theta), np.sin(theta), 'b-', label='单位圆')
    circle_ax.scatter([0.5, -0.3], [0.3, -0.4], marker='x', c='r', s=100, label='极点')
    circle_ax.scatter([0.7], [-0.2], marker='o', c='g', s=80, label='零点')
    circle_ax.set_title('零极点分布', fontsize=10)
    circle_ax.legend(fontsize=8)
    circle_ax.grid(True)
    circle_ax.set_xlim(-1.2, 1.2)
    circle_ax.set_ylim(-1.2, 1.2)

    # 添加系统框图
    box_ax = fig.add_axes([0.75, 0.1, 0.2, 0.2])
    box_ax.set_xlim(0, 5)
    box_ax.set_ylim(0, 5)
    box_ax.axis('off')

    # 绘制简单系统框图
    box_ax.text(1, 3, '输入\nx[n]', ha='center', va='center',
                bbox=dict(boxstyle="round", fc="lightyellow"))
    box_ax.arrow(1.5, 3, 1, 0, head_width=0.2, head_length=0.1, fc='k', ec='k')
    box_ax.text(2.5, 3, '系统\nH(z)', ha='center', va='center',
                bbox=dict(boxstyle="round", fc="lightblue"))
    box_ax.arrow(3.5, 3, 1, 0, head_width=0.2, head_length=0.1, fc='k', ec='k')
    box_ax.text(4.5, 3, '输出\ny[n]', ha='center', va='center',
                bbox=dict(boxstyle="round", fc="lightgreen"))
    box_ax.set_title('系统框图', fontsize=10)

    # 添加图例
    plt.figtext(0.5, 0.02, '图1：离散系统z域分析思维导图', ha='center', fontsize=14)

    plt.savefig('z_domain_mindmap.png', dpi=300, bbox_inches='tight')
    plt.show()


# 创建思维导图
create_z_domain_mindmap()