求助：基于Bode图确定截止频率、误差常数及裕度的代码修正

问题背景 #

通过CSV文件中的幅值、相位、频率数据绘制Bode图后，需计算系统的静态误差常数、增益裕度、相位裕度及截止频率，但原有代码的计算逻辑存在错误，导致参数结果不准确。

现有代码的核心问题 #

• 截止频率计算错误：np.argmax(magnitude > -3)仅返回第一个满足幅值>-3dB的索引，未定位到幅值从高于-3dB跌落至低于-3dB的交叉点
• 增益/相位裕度逻辑混淆：
  • 错误将幅值过0dB的频率作为增益裕度对应频率，实际增益裕度对应相位为-180°时的幅值
  • 错误取幅值为0dB处的相位作为相位裕度，实际相位裕度是该相位与-180°的差值
• 带宽计算逻辑混乱：未定义变量w_deg，且带宽计算逻辑不符合定义
• 未处理静态误差常数：缺少从低频段幅值推导静态误差常数的逻辑

修正后的完整代码 #

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

# 读取数据
data = pd.read_csv('example_P3.csv')
magnitude = data['Mag'].values
phase = data['Phase'].values
w_rad = data['W'].values
w_hz = w_rad / (2 * np.pi)

# --------------------------
# 1. 计算静态误差常数（Kp/Kv/Ka）
# --------------------------
# 取低频段数据（前10%的频率点，可根据实际调整）
low_freq_idx = int(len(w_rad) * 0.1)
low_w = w_rad[:low_freq_idx]
low_mag = magnitude[:low_freq_idx]

# 计算低频段幅值斜率（dB/dec）
log_w = np.log10(low_w)
log_mag = low_mag / 20  # 转换为线性幅值的对数
slope, _ = np.polyfit(log_w, log_mag, 1)
slope_db = slope * 20  # 转换为dB/dec

static_constant = None
constant_type = ""
if np.isclose(slope_db, 0, atol=1):
    # 0型系统，Kp = 10^(Mag_dB/20)
    avg_low_mag = np.mean(low_mag)
    Kp = 10 ** (avg_low_mag / 20)
    static_constant = Kp
    constant_type = "位置误差常数Kp"
elif np.isclose(slope_db, -20, atol=1):
    # 1型系统，Kv = 10^((Mag_dB + 20logω)/20)
    mid_idx = low_freq_idx // 2
    Kv = 10 ** ((low_mag[mid_idx] + 20 * np.log10(low_w[mid_idx])) / 20)
    static_constant = Kv
    constant_type = "速度误差常数Kv"
elif np.isclose(slope_db, -40, atol=1):
    # 2型系统，Ka = 10^((Mag_dB + 40logω)/20)
    mid_idx = low_freq_idx // 2
    Ka = 10 ** ((low_mag[mid_idx] + 40 * np.log10(low_w[mid_idx])) / 20)
    static_constant = Ka
    constant_type = "加速度误差常数Ka"

# --------------------------
# 2. 计算截止频率（-3dB频率）
# --------------------------
# 找到幅值从高于-3dB降到低于-3dB的交叉点
cross_idx = np.where(np.diff(np.sign(magnitude + 3)))[0]
if len(cross_idx) > 0:
    idx = cross_idx[0]
    # 线性插值计算精准截止频率
    f_interp = interp1d(magnitude[idx:idx+2], w_hz[idx:idx+2], kind='linear')
    cutoff_freq_hz = f_interp(-3)
else:
    # 若未找到交叉点，取幅值最接近-3dB的频率
    closest_idx = np.argmin(np.abs(magnitude + 3))
    cutoff_freq_hz = w_hz[closest_idx]

# --------------------------
# 3. 计算相位裕度（PM）和截止频率处的相位
# --------------------------
# 找到幅值穿越0dB的频率（ωc）
cross_0_idx = np.where(np.diff(np.sign(magnitude)))[0]
if len(cross_0_idx) > 0:
    idx_0 = cross_0_idx[0]
    f_interp_0 = interp1d(magnitude[idx_0:idx_0+2], w_hz[idx_0:idx_0+2], kind='linear')
    wc_hz = f_interp_0(0)
    # 插值获取该频率下的相位
    phase_interp = interp1d(w_hz, phase, kind='linear')
    phase_at_wc = phase_interp(wc_hz)
    phase_margin = phase_at_wc + 180  # PM = φ(ωc) - (-180°)
else:
    closest_0_idx = np.argmin(np.abs(magnitude))
    wc_hz = w_hz[closest_0_idx]
    phase_at_wc = phase[closest_0_idx]
    phase_margin = phase_at_wc + 180

# --------------------------
# 4. 计算增益裕度（GM）和相位穿越-180°的频率
# --------------------------
# 找到相位穿越-180°的频率
cross_180_idx = np.where(np.diff(np.sign(phase + 180)))[0]
gain_margin_db = np.nan
w180_hz = np.nan
if len(cross_180_idx) > 0:
    idx_180 = cross_180_idx[0]
    f_interp_180 = interp1d(phase[idx_180:idx_180+2], w_hz[idx_180:idx_180+2], kind='linear')
    w180_hz = f_interp_180(-180)
    # 插值获取该频率下的幅值
    mag_interp = interp1d(w_hz, magnitude, kind='linear')
    mag_at_w180 = mag_interp(w180_hz)
    gain_margin_db = -mag_at_w180  # GM(dB) = -Mag(ω180)

# --------------------------
# 绘图标注
# --------------------------
plt.figure(figsize=(12, 8))

# 幅值图
plt.subplot(2, 1, 1)
plt.semilogx(w_hz, magnitude, 'b-', label='幅值曲线')
plt.axhline(0, color='gray', linestyle=':', label='0 dB')
plt.axhline(-3, color='orange', linestyle=':', label='-3 dB')
plt.axvline(cutoff_freq_hz, color='y', linestyle='--', label=f'截止频率: {cutoff_freq_hz:.2f} Hz')
plt.axvline(wc_hz, color='c', linestyle='--', label=f'幅值穿越频率: {wc_hz:.2f} Hz')
plt.title('Bode图 - 幅值')
plt.ylabel('幅值 (dB)')
plt.legend()
plt.grid(True, which="both", ls="-")

# 相位图
plt.subplot(2, 1, 2)
plt.semilogx(w_hz, phase, 'r-', label='相位曲线')
plt.axhline(-180, color='gray', linestyle=':', label='-180°')
plt.axvline(wc_hz, color='c', linestyle='--', label=f'幅值穿越频率: {wc_hz:.2f} Hz')
if not np.isnan(w180_hz):
    plt.axvline(w180_hz, color='g', linestyle='--', label=f'相位穿越频率: {w180_hz:.2f} Hz')
plt.title('Bode图 - 相位')
plt.xlabel('频率 (Hz)')
plt.ylabel('相位 (°)')
plt.legend()
plt.grid(True, which="both", ls="-")

plt.tight_layout()
plt.show()

# 打印计算结果
print(f"静态误差常数({constant_type}): {static_constant:.4f}")
print(f"截止频率(-3dB): {cutoff_freq_hz:.2f} Hz")
print(f"相位裕度: {phase_margin:.2f} °")
if not np.isnan(gain_margin_db):
    print(f"增益裕度: {gain_margin_db:.2f} dB")
else:
    print("未找到相位穿越-180°的点，无法计算增益裕度")

关键参数的正确计算逻辑 #

1. 静态误差常数：
   • 通过低频段幅值斜率判断系统型别（0型/1型/2型），再根据对应公式计算位置/速度/加速度误差常数
2. 截止频率：
   • 定位幅值穿越-3dB的交叉点，通过线性插值获取精准频率，避免直接取离散点的误差
3. 相位裕度：
   • 先找到幅值穿越0dB的频率（ωc），计算该频率下的相位与-180°的差值，即PM = φ(ωc) + 180°
4. 增益裕度：
   • 找到相位穿越-180°的频率（ω180），计算该频率下幅值的dB值的绝对值，即GM(dB) = -Mag(ω180)

内容的提问来源于stack exchange，提问作者Luis Manuel Martínez Gómez