# 从理论到代码：手把手实现Sanathanan–Koerner矢量拟合算法 如果你正在处理电路仿真、信号系统建模或者任何需要从频域数据中提取有理函数模型的工作，那么“矢量拟合”这个词对你来说一定不陌生。面对一堆离散的频率响应数据点，如何找到一个既简洁又准确的传递函数来描述系统行为，是许多工程师和研究人员日常要解决的痛点。传统的直接拟合方法往往受困于数值不稳定，或者拟合出的模型物理意义不明确。这时，Sanathanan–Koerner算法提供了一条优雅的路径，它将一个棘手的非线性拟合问题，转化为一系列可迭代求解的线性最小二乘问题。 这篇文章就是为你准备的——无论你是刚接触系统辨识的新手，还是想将算法理论付诸实践的开发者。我们将彻底抛开复杂的公式推导，聚焦于如何用Python一步步构建一个健壮、可用的Sanathanan–Koerner算法实现。我会分享在编码过程中遇到的典型“坑”，比如数值条件数恶化如何处理、迭代何时终止，并提供完整的、可直接运行的代码示例。我们的目标很明确：让你不仅能理解算法，更能亲手把它用起来。 ## 1. 理解核心：Sanathanan–Koerner算法解决了什么？ 在深入代码之前，我们得先搞清楚要对付的“敌人”是什么。想象一下，你通过仿真或测量，得到了一组系统的频率响应数据：在N个不同的角频率点 ω₁, ω₂, ..., ω_N 上，你知道了对应的复数响应值 H(ω₁), H(ω₂), ..., H(ω_N)。这些数据可能来自一个复杂的RLC网络、一段传输线，或者任何线性时不变系统。 我们的目标是找到一个有理传递函数 H̃(s) 来逼近这些数据： ``` H̃(s) = N(s) / D(s) = (a₀ + a₁s + ... + aₙsⁿ) / (b₀ + b₁s + ... + bₙsⁿ) ``` 其中 s = jω。这本质上是一个模型参数（系数 a_i, b_i）估计问题。最直观的想法是最小化所有频率点上的误差平方和： ``` 目标函数 = Σ |H(ω_k) - N(jω_k)/D(jω_k)|² ``` 麻烦在于，目标函数关于分母系数 b_i 是非线性的。如果你尝试直接用非线性优化算法（如Levenberg-Marquardt）去求解，很容易陷入局部最优，或者对初始值极其敏感，计算量也很大。 Sanathanan和Koerner在1960年代提出的方法巧妙地绕开了这个非线性难题。它的核心思想是**迭代重加权**。算法从一个简单的线性化问题（即Levy法）开始，得到模型参数的初始估计。然后，在后续的迭代中，它利用上一轮迭代得到的分母多项式，对误差函数进行“重新加权”，从而将原问题转化为一个新的、关于当前迭代参数的线性最小二乘问题。 具体来说，在第 i 次迭代中，我们求解的是： ``` 最小化 Σ | [H(ω_k)*D⁽ⁱ⁾(jω_k) - N⁽ⁱ⁾(jω_k)] / D⁽ⁱ⁻¹⁾(jω_k) |² ``` 注意，分母上的 D⁽ⁱ⁻¹⁾(jω_k) 来自上一次迭代的结果，是已知的常数。因此，对于当前待求的系数 a⁽ⁱ⁾_i 和 b⁽ⁱ⁾_i 来说，这**仍然是一个线性最小二乘问题**。随着迭代进行，权重项 D⁽ⁱ⁻¹⁾(jω_k) 会不断更新，引导解向真实的最小二乘解收敛。 > **提示**：你可以把每次迭代看作是对模型“分母”的一次修正。初始的Levy拟合相当于假设分母为1（权重均匀），而后续迭代则根据当前对分母的最佳估计，重新调整各个频率点误差的“重要性”，从而逐步逼近真实的最优分母。 ## 2. 动手准备：Python环境与数据构造 理论清晰后，我们开始搭建实战环境。我将使用最通用的科学计算栈，确保代码的可复现性。 ### 2.1 环境配置与核心库 首先，确保你的Python环境安装了以下库。我推荐使用Anaconda来管理环境，避免依赖冲突。 ```bash # 创建并激活一个专用于信号处理的新环境（可选） conda create -n vector_fitting python=3.9 conda activate vector_fitting # 安装核心库 pip install numpy scipy matplotlib ``` 接下来，在Python脚本或Jupyter Notebook的开头，导入我们将要用到的模块： ```python import numpy as np import scipy.linalg as la import matplotlib.pyplot as plt from typing import Tuple, Optional # 设置绘图风格，让图表更美观 plt.style.use('seaborn-v0_8-darkgrid') np.set_printoptions(precision=4, suppress=True) ``` ### 2.2 生成用于测试的合成数据 为了验证我们的算法，我们需要一组“标准答案”已知的数据。这里，我们构造一个已知极点、零点的系统，计算出其精确的频率响应，并添加一点噪声来模拟真实测量。 ```python def generate_test_data( num_points: int = 200, freq_start: float = 1e0, freq_end: float = 1e4, noise_level: float = 0.01 ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ 生成用于测试矢量拟合算法的合成频率响应数据。 参数: num_points: 频率点数 freq_start: 起始频率 (Hz) freq_end: 终止频率 (Hz) noise_level: 添加到响应数据中的相对噪声水平 返回: freq: 角频率数组 (rad/s) s: 复频率点数组 (s = jω) H_exact: 无噪声的精确频率响应 H_noisy: 添加了噪声的频率响应（模拟测量数据） """ # 在对数尺度上均匀分布频率点 freq = np.logspace(np.log10(freq_start), np.log10(freq_end), num_points) omega = 2 * np.pi * freq # 转换为角频率 s = 1j * omega # 复频率 # 定义一个已知的传递函数：H(s) = (s + 1000) / [(s + 100)(s + 500)] # 零点: z = -1000 # 极点: p1 = -100, p2 = -500 # 增益: k = 1 # 计算精确响应 H_exact = (s + 1000) / ((s + 100) * (s + 500)) # 添加复高斯噪声，模拟测量误差 noise_real = np.random.randn(num_points) * noise_level * np.abs(H_exact) / np.sqrt(2) noise_imag = np.random.randn(num_points) * noise_level * np.abs(H_exact) / np.sqrt(2) H_noisy = H_exact + noise_real + 1j * noise_imag return omega, s, H_exact, H_noisy # 生成数据 omega, s, H_exact, H_meas = generate_test_data(num_points=100, noise_level=0.02) ``` 我们可以快速可视化一下生成的数据，看看它的幅频和相频特性： ```python fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8)) # 幅频响应 ax1.semilogx(omega, 20*np.log10(np.abs(H_exact)), 'b-', linewidth=2, label='Exact') ax1.semilogx(omega, 20*np.log10(np.abs(H_meas)), 'ro', markersize=4, alpha=0.6, label='Measured (Noisy)') ax1.set_ylabel('Magnitude (dB)') ax1.set_title('Synthetic Frequency Response Data for Testing') ax1.legend() ax1.grid(True, which="both", ls="-", alpha=0.2) # 相频响应 ax2.semilogx(omega, np.angle(H_exact, deg=True), 'b-', linewidth=2) ax2.semilogx(omega, np.angle(H_meas, deg=True), 'ro', markersize=4, alpha=0.6) ax2.set_xlabel('Angular Frequency (rad/s)') ax2.set_ylabel('Phase (degrees)') ax2.grid(True, which="both", ls="-", alpha=0.2) plt.tight_layout() plt.show() ``` ## 3. 算法实现：构建Sanathanan–Koerner拟合器 现在进入最核心的部分：实现算法本身。我们将把算法封装成一个类 `SanathananKoernerFitter`，这样便于管理状态、重复使用和扩展。 ### 3.1 初始化与问题构建 类的初始化需要指定我们想要拟合的有理函数阶数。这里我们假设分子和分母同阶，这是最常见的情况。 ```python class SanathananKoernerFitter: """ 使用Sanathanan-Koerner迭代方法实现矢量拟合（有理函数拟合）。 """ def __init__(self, order_n: int): """ 初始化拟合器。 参数: order_n: 传递函数分子和分母多项式的阶数。 """ if order_n < 1: raise ValueError("Polynomial order must be at least 1.") self.order_n = order_n # 初始化系数为None，在拟合后赋值 self.a_coeffs = None # 分子系数 [a0, a1, ..., an] self.b_coeffs = None # 分母系数 [b0, b1, ..., bn]，通常bn=1 def _build_matrix(self, s: np.ndarray, weights: Optional[np.ndarray] = None) -> np.ndarray: """ 为线性最小二乘问题构建系数矩阵A。 对于方程：H(s)*D(s) - N(s) ≈ 0，其中D(s)=Σ b_i s^i, N(s)=Σ a_i s^i。 将其重写为：Σ a_i s^i - H(s) * Σ b_i s^i ≈ 0。 对于每个频率点s_k，我们得到一个方程，未知数是[a0,...,an, b0,..., b_{n-1}]（假设bn=1）。 参数: s: 复频率点数组 (shape: (M,)) weights: 每个方程的权重数组 (shape: (M,))。如果为None，则所有权重为1。 返回: A: 系数矩阵 (shape: (M, 2*order_n + 1)) """ M = len(s) n = self.order_n A = np.zeros((M, 2 * n + 1), dtype=np.complex128) # 构建分子部分对应的列 (a_i * s^i) for i in range(n + 1): A[:, i] = s ** i # 构建分母部分对应的列 (-H(s) * b_i * s^i)，注意i从0到n-1，b_n固定为1 # 注意：H(s)在调用此函数时尚未传入，因此这部分需要在外部结合H(s)处理。 # 这里先构建 s^i 的部分，符号和H(s)的乘法在外部完成。 for i in range(n): # b0 到 b_{n-1} A[:, n + 1 + i] = - (s ** i) # 先占位，外部会乘以H(s) # 如果提供了权重，对矩阵的每一行进行加权 if weights is not None: if len(weights) != M: raise ValueError("Weights length must match number of frequency points.") # 将权重转换为列向量并广播到所有列 W = weights.reshape(-1, 1) A = A * W return A ``` 这里有一个关键设计：我们将矩阵构建中依赖于测量值 H(s) 的部分分离出来。这是因为在Sanathanan-Koerner迭代中，权重项（即上一次迭代的分母）是变化的，并且需要与 H(s) 相乘。我们将在拟合主循环中处理这部分逻辑。 ### 3.2 核心迭代拟合过程 这是算法的引擎。我们遵循标准的迭代流程：Levy初始化 -> 加权迭代 -> 收敛判断。 ```python def fit(self, s: np.ndarray, H: np.ndarray, max_iter: int = 20, tol: float = 1e-6) -> dict: """ 执行Sanathanan-Koerner迭代拟合。 参数: s: 复频率点数组 (1D complex array) H: 在频率点s处测量的频率响应数组 (1D complex array) max_iter: 最大迭代次数 tol: 迭代收敛容差（相邻迭代间系数变化的范数） 返回: result: 包含拟合结果的字典，包括系数、迭代历史等。 """ M = len(s) if len(H) != M: raise ValueError("Length of s and H must match.") n = self.order_n num_coeffs = 2 * n + 1 # a0...an, b0...b_{n-1}， b_n固定为1 # 初始化：第一次迭代使用Levy方法（权重全为1） print("Iteration 1 (Levy initialization)...") A = self._build_matrix(s, weights=None) # 无权重 # 构建右端项：对于方程 A * x = b，我们的方程是 A*x ≈ 0，但矩阵A中分母部分列需要乘以H(s) # 实际上，我们的方程是： Σ a_i s^i - H(s) * Σ b_i s^i ≈ 0 # 所以我们需要构造一个矩阵A_levy，其中分母部分的列已经乘以了-H(s) A_levy = A.copy() for i in range(n): # 对应b_i的列 A_levy[:, n + 1 + i] = -H * (s ** i) # 目标是最小化 ||A_levy * x||^2，这是一个齐次方程。 # 为了避免零解，我们需要施加一个约束。常见做法是固定分母的最高次项系数为1 (b_n=1)。 # 这等价于从待求解向量中移除b_n，并将其贡献移到方程右侧。 # 我们的矩阵A_levy的最后一列对应的是 -H(s) * s^n * b_n。设b_n=1，则这一列变为已知项，移到右侧。 # 因此，我们求解的方程变为： A_reduced * x_reduced = rhs # 其中，x_reduced = [a0,..., an, b0,..., b_{n-1}]^T # A_reduced 是A_levy的前 (2n+1) 列 # rhs = - (A_levy的最后一列) ，因为原来 A_levy * [x_reduced; 1] = 0 A_reduced = A_levy[:, :-1] # 去掉对应b_n的最后一列 rhs = -A_levy[:, -1].reshape(-1, 1) # 将b_n=1的贡献移到右边，变为负号 # 求解线性最小二乘问题 x, residuals, rank, sing_vals = la.lstsq(A_reduced, rhs, lapack_driver='gelsy') x = x.flatten() # 提取系数 a_current = x[:n+1] # a0, a1, ..., an b_current = np.append(x[n+1:], 1.0) # b0, b1, ..., b_{n-1}, b_n=1 # 存储迭代历史 history = { 'a_coeffs': [a_current.copy()], 'b_coeffs': [b_current.copy()], 'residual_norms': [np.linalg.norm(residuals) if len(residuals) > 0 else np.nan] } # Sanathanan-Koerner 迭代 for iter_num in range(2, max_iter + 1): # 计算当前迭代的权重：1 / D_{prev}(s)，其中D_{prev}(s)是上一次迭代的分母多项式值 D_prev = np.polyval(b_current[::-1], s) # np.polyval期望系数从高次到低次 weights = 1.0 / D_prev # 构建加权矩阵 A_weighted = self._build_matrix(s, weights=weights) # 同样，处理分母部分与H(s)的乘法 for i in range(n): A_weighted[:, n + 1 + i] = -H * (s ** i) * weights # 同样固定b_n=1，构建约化系统 A_w_reduced = A_weighted[:, :-1] rhs_w = -A_weighted[:, -1].reshape(-1, 1) # 对应b_n=1的项 # 求解加权最小二乘问题 x_new, residuals, rank, sing_vals = la.lstsq(A_w_reduced, rhs_w, lapack_driver='gelsy') x_new = x_new.flatten() a_new = x_new[:n+1] b_new = np.append(x_new[n+1:], 1.0) # 计算系数变化以检查收敛 delta_a = np.linalg.norm(a_new - a_current) / np.linalg.norm(a_current) delta_b = np.linalg.norm(b_new - b_current) / np.linalg.norm(b_current) delta_max = max(delta_a, delta_b) # 更新当前系数 a_current, b_current = a_new, b_new history['a_coeffs'].append(a_current.copy()) history['b_coeffs'].append(b_current.copy()) history['residual_norms'].append(np.linalg.norm(residuals) if len(residuals) > 0 else np.nan) print(f"Iteration {iter_num}: max coefficient change = {delta_max:.2e}") if delta_max < tol: print(f"Converged after {iter_num} iterations.") break else: print(f"Stopped after reaching maximum iterations ({max_iter}).") # 保存最终系数 self.a_coeffs = a_current self.b_coeffs = b_current result = { 'a_coeffs': self.a_coeffs, 'b_coeffs': self.b_coeffs, 'history': history, 'final_iter': len(history['a_coeffs']) } return result ``` ### 3.3 模型评估与后处理工具 拟合完成后，我们需要工具来评估模型质量，并可能进行进一步处理，如提取极点/零点。 ```python def predict(self, s: np.ndarray) -> np.ndarray: """ 使用拟合的系数计算有理函数在给定复频率点上的预测值。 参数: s: 复频率点数组 返回: H_pred: 预测的频率响应数组 """ if self.a_coeffs is None or self.b_coeffs is None: raise RuntimeError("Model has not been fitted yet. Call 'fit' first.") # 计算分子和分母多项式的值。注意np.polyval期望系数从高次到低次排列。 N = np.polyval(self.a_coeffs[::-1], s) D = np.polyval(self.b_coeffs[::-1], s) return N / D def calculate_error(self, s: np.ndarray, H: np.ndarray, metric: str = 'relative') -> np.ndarray: """ 计算拟合误差。 参数: s: 频率点 H: 测量响应 metric: 误差度量方式。'absolute'为绝对误差，'relative'为相对误差。 返回: error: 在每个频率点上的误差值。 """ H_pred = self.predict(s) if metric == 'absolute': error = np.abs(H - H_pred) elif metric == 'relative': # 避免除以零，使用测量值的幅度作为分母 denom = np.abs(H) denom[denom == 0] = 1e-12 # 一个很小的数 error = np.abs(H - H_pred) / denom else: raise ValueError("metric must be 'absolute' or 'relative'") return error def get_poles_zeros(self) -> Tuple[np.ndarray, np.ndarray]: """ 从拟合的有理函数系数中计算极点和零点。 返回: poles: 极点数组（分母多项式的根） zeros: 零点数组（分子多项式的根） """ if self.a_coeffs is None or self.b_coeffs is None: raise RuntimeError("Model has not been fitted yet.") # 求根。注意np.roots期望系数从高次到低次排列。 poles = np.roots(self.b_coeffs[::-1]) zeros = np.roots(self.a_coeffs[::-1]) return poles, zeros ``` ## 4. 实战演练：从拟合到分析 现在，让我们把所有的部件组装起来，进行一次完整的拟合流程演示。 ### 4.1 执行拟合并可视化迭代过程 ```python # 实例化拟合器，假设我们已知系统是2阶（或我们想用2阶模型去拟合） fitter = SanathananKoernerFitter(order_n=2) # 执行拟合 result = fitter.fit(s, H_meas, max_iter=15, tol=1e-9) # 获取拟合的系数 a_fit = result['a_coeffs'] b_fit = result['b_coeffs'] print("\nFitted numerator coefficients (a):", a_fit) print("Fitted denominator coefficients (b):", b_fit) # 计算拟合模型的预测值 H_fit = fitter.predict(s) # 绘制拟合结果对比 fig, axes = plt.subplots(3, 1, figsize=(12, 12), sharex=True) # 1. 幅频响应对比 ax = axes[0] ax.semilogx(omega, 20*np.log10(np.abs(H_exact)), 'k-', linewidth=3, label='Exact System', alpha=0.7) ax.semilogx(omega, 20*np.log10(np.abs(H_meas)), 'bo', markersize=5, alpha=0.5, label='Measured Data') ax.semilogx(omega, 20*np.log10(np.abs(H_fit)), 'r--', linewidth=2.5, label='SK Fit (n=2)') ax.set_ylabel('Magnitude (dB)') ax.legend(loc='best') ax.grid(True, which="both", ls="-", alpha=0.2) ax.set_title('Sanathanan-Koerner Fitting Result') # 2. 相频响应对比 ax = axes[1] ax.semilogx(omega, np.angle(H_exact, deg=True), 'k-', linewidth=3, alpha=0.7) ax.semilogx(omega, np.angle(H_meas, deg=True), 'bo', markersize=5, alpha=0.5) ax.semilogx(omega, np.angle(H_fit, deg=True), 'r--', linewidth=2.5) ax.set_ylabel('Phase (deg)') ax.grid(True, which="both", ls="-", alpha=0.2) # 3. 相对误差 ax = axes[2] rel_error = fitter.calculate_error(s, H_meas, metric='relative') ax.semilogx(omega, 100 * rel_error, 'g-', linewidth=2) ax.set_xlabel('Angular Frequency (rad/s)') ax.set_ylabel('Relative Error (%)') ax.grid(True, which="both", ls="-", alpha=0.2) ax.axhline(y=1.0, color='gray', linestyle=':', label='1% error line') ax.legend() plt.tight_layout() plt.show() ``` ### 4.2 分析迭代收敛性与模型极点/零点 查看迭代历史可以帮助我们理解算法的收敛行为，而检查极点/零点则可以验证模型的稳定性与物理合理性。 ```python # 分析迭代历史 history = result['history'] iter_nums = range(1, result['final_iter'] + 1) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # 绘制残差范数变化 ax1.plot(iter_nums, history['residual_norms'], 's-', linewidth=2, markersize=8) ax1.set_xlabel('Iteration Number') ax1.set_ylabel('Residual Norm (||Ax-b||)') ax1.set_title('Convergence of Residual') ax1.grid(True) ax1.set_xlim(0.5, result['final_iter']+0.5) ax1.set_xticks(iter_nums) # 绘制系数变化（以第一次迭代后的系数为基准） coeff_changes = [] for i in range(1, len(history['a_coeffs'])): delta_a = np.linalg.norm(history['a_coeffs'][i] - history['a_coeffs'][i-1]) delta_b = np.linalg.norm(history['b_coeffs'][i] - history['b_coeffs'][i-1]) coeff_changes.append(max(delta_a, delta_b)) ax2.semilogy(range(2, result['final_iter']+1), coeff_changes, 'o-', linewidth=2, markersize=8) ax2.axhline(y=1e-9, color='r', linestyle='--', label=f'Tolerance ({1e-9})') ax2.set_xlabel('Iteration Number') ax2.set_ylabel('Max Coefficient Change (log scale)') ax2.set_title('Coefficient Change Between Iterations') ax2.legend() ax2.grid(True) ax2.set_xlim(1.5, result['final_iter']+0.5) plt.tight_layout() plt.show() # 提取并打印极点/零点 poles, zeros = fitter.get_poles_zeros() print("\n=== Model Poles and Zeros ===") print("Poles (should be stable, i.e., have negative real parts):") for p in poles: print(f" {p:.4e}") print("\nZeros:") for z in zeros: print(f" {z:.4e}") # 检查稳定性：所有极点实部应为负（连续时间系统） if np.all(np.real(poles) < 0): print("\n✓ The fitted model is stable (all poles in left half-plane).") else: print("\n⚠ Warning: The fitted model has poles in the right half-plane or on the imaginary axis, indicating potential instability.") ``` ### 4.3 处理数值稳定性：缩放与正则化 在实际应用中，尤其是频率跨度大（如从Hz到GHz）或模型阶数高时，直接构造的Vandermonde类矩阵条件数会非常差，导致最小二乘求解不稳定。一个简单而有效的技巧是对频率数据进行**缩放**。 ```python def fit_with_scaling(self, s: np.ndarray, H: np.ndarray, max_iter: int = 20, tol: float = 1e-6) -> dict: """ 带频率缩放的拟合版本，以改善数值条件。 缩放因子通常取频率范围的几何平均值。 """ # 计算缩放因子：频率幅值的几何平均 scale = np.exp(np.mean(np.log(np.abs(s) + 1e-12))) # 避免log(0) if np.isclose(scale, 0): scale = 1.0 print(f"Using frequency scaling factor: {scale:.4e}") # 缩放频率数据 s_scaled = s / scale # 注意：缩放会影响传递函数。为了保持H(s)不变，我们需要同时缩放系数。 # 更简单的方法是：用缩放后的s进行拟合，最后再对系数进行反缩放。 result_scaled = self.fit(s_scaled, H, max_iter, tol) # 反缩放系数以得到原始s域的系数 n = self.order_n a_orig = self.a_coeffs.copy() b_orig = self.b_coeffs.copy() for i in range(n+1): a_orig[i] /= (scale ** i) b_orig[i] /= (scale ** i) self.a_coeffs = a_orig self.b_coeffs = b_orig result_scaled['a_coeffs'] = a_orig result_scaled['b_coeffs'] = b_orig return result_scaled # 将这个方法添加到我们的类中 SanathananKoernerFitter.fit_with_scaling = fit_with_scaling ``` 让我们用一个更具挑战性的例子来测试缩放的效果：一个更高阶的系统，频率跨度更大。 ```python # 生成一个更高阶、更宽频带的测试系统 def higher_order_system(s): # 一个4阶系统，具有两个实数极点和一对共轭复极点 # H(s) = (s^2 + 500s + 1e6) / [(s+100)(s+1000)(s^2 + 300s + 1e5)] num = np.polyval([1, 500, 1e6], s) # 分子: s^2 + 500s + 1e6 den = np.polyval([1, 1400, 441000, 4.1e7, 1e10], s) # 分母展开后系数 return num / den # 生成数据，频率跨越4个数量级 omega_hard = np.logspace(1, 5, 300) # 10^1 到 10^5 rad/s s_hard = 1j * omega_hard H_exact_hard = higher_order_system(s_hard) # 添加少量噪声 noise = 0.005 * (np.random.randn(len(s_hard)) + 1j*np.random.randn(len(s_hard))) * np.abs(H_exact_hard) H_meas_hard = H_exact_hard + noise # 尝试用4阶模型拟合 print("=== Fitting High-Order System with Wide Frequency Range ===") fitter_hard = SanathananKoernerFitter(order_n=4) # 1. 不使用缩放（可能数值不稳定） print("\n--- Attempt without scaling ---") try: result_no_scale = fitter_hard.fit(s_hard, H_meas_hard, max_iter=15, tol=1e-8) H_fit_no_scale = fitter_hard.predict(s_hard) error_no_scale = np.mean(np.abs(H_fit_no_scale - H_exact_hard) / np.abs(H_exact_hard)) print(f"Mean relative error (no scaling): {error_no_scale:.2%}") except np.linalg.LinAlgError as e: print(f"Fitting failed without scaling: {e}") # 2. 使用缩放 print("\n--- Attempt with frequency scaling ---") fitter_hard_scaled = SanathananKoernerFitter(order_n=4) result_scaled = fitter_hard_scaled.fit_with_scaling(s_hard, H_meas_hard, max_iter=15, tol=1e-8) H_fit_scaled = fitter_hard_scaled.predict(s_hard) error_scaled = np.mean(np.abs(H_fit_scaled - H_exact_hard) / np.abs(H_exact_hard)) print(f"Mean relative error (with scaling): {error_scaled:.2%}") # 可视化对比 fig, ax = plt.subplots(1, 2, figsize=(15, 5)) # 幅频响应 ax[0].semilogx(omega_hard, 20*np.log10(np.abs(H_exact_hard)), 'k-', label='Exact', linewidth=3, alpha=0.7) ax[0].semilogx(omega_hard, 20*np.log10(np.abs(H_meas_hard)), 'bo', alpha=0.4, markersize=3, label='Measured') ax[0].semilogx(omega_hard, 20*np.log10(np.abs(H_fit_scaled)), 'r--', label='SK Fit (with scaling)', linewidth=2) ax[0].set_xlabel('Angular Frequency (rad/s)') ax[0].set_ylabel('Magnitude (dB)') ax[0].legend() ax[0].grid(True, which="both", ls="-", alpha=0.2) ax[0].set_title('Fitting Result with Scaling') # 误差对比 rel_error_scaled = np.abs(H_fit_scaled - H_exact_hard) / np.abs(H_exact_hard) ax[1].semilogx(omega_hard, 100*rel_error_scaled, 'g-', linewidth=2) ax[1].axhline(y=1.0, color='gray', linestyle=':', label='1% error') ax[1].set_xlabel('Angular Frequency (rad/s)') ax[1].set_ylabel('Relative Error (%)') ax[1].legend() ax[1].grid(True, which="both", ls="-", alpha=0.2) ax[1].set_title('Fitting Error (with Scaling)') plt.tight_layout() plt.show() ``` 从我的经验来看，对于宽频带拟合，**频率缩放几乎是必需的**。它能将矩阵的条件数降低好几个数量级，使得求解过程稳定，最终拟合误差也显著减小。如果不做处理，高阶拟合很容易因为数值问题而完全失败，或者得到没有物理意义的极点（例如实部为正的不稳定极点）。 ## 5. 进阶话题与工程实践建议 掌握了基本实现后，我们来看看如何让这个工具在实际项目中更可靠、更强大。 ### 5.1 模型阶数选择：一个实际的难题 应该选择几阶模型？这是一个偏差-方差权衡的问题。阶数太低，模型无法捕捉系统动态（欠拟合）；阶数太高，会拟合噪声，导致模型在训练数据外表现差（过拟合）。这里介绍两种实用方法： **方法一：基于误差下降的启发式方法** 随着阶数增加，拟合误差会下降，但到某个点后下降会变得平缓。我们可以寻找这个“拐点”。 ```python def find_optimal_order(s, H, max_order=8, plot=True): """ 通过扫描不同模型阶数，寻找一个合理的阶数。 """ orders = range(1, max_order+1) train_errors = [] # 简单起见，我们用所有数据拟合并评估。更严谨的做法应使用交叉验证。 for n in orders: fitter = SanathananKoernerFitter(order_n=n) # 使用缩放以保持数值稳定 _ = fitter.fit_with_scaling(s, H, max_iter=15, tol=1e-9) H_pred = fitter.predict(s) # 计算相对误差的均方根 rel_error = np.abs(H - H_pred) / np.abs(H) rms_error = np.sqrt(np.mean(rel_error**2)) train_errors.append(rms_error) print(f"Order {n}: RMS relative error = {rms_error:.3%}") if plot: plt.figure(figsize=(10, 6)) plt.plot(orders, train_errors, 'bo-', linewidth=2, markersize=8) plt.xlabel('Model Order (n)') plt.ylabel('RMS Relative Error (on training data)') plt.title('Model Order Selection: Error vs. Complexity') plt.grid(True) plt.xticks(orders) plt.show() # 一个简单的启发式：选择误差下降开始变缓的阶数 errors = np.array(train_errors) # 计算误差的相对下降率 rel_reduction = -np.diff(errors) / errors[:-1] # 找到下降率显著小于前一个点的位置（例如，下降率小于平均下降率的一半） if len(rel_reduction) > 0: avg_reduction = np.mean(rel_reduction) # 找到第一个下降率小于平均下降率 * 0.5 的阶数，其前一个阶数可能是较优选择 for i, red in enumerate(rel_reduction): if red < avg_reduction * 0.5: suggested_order = i + 1 # 因为i从0开始，对应从1阶到2阶的变化 print(f"\nSuggested order based on reduction slowdown: {suggested_order}") return suggested_order print(f"\nNo clear slowdown detected. Using maximum order tested: {max_order}") return max_order # 使用之前生成的数据进行测试 optimal_n = find_optimal_order(s, H_meas, max_order=6) ``` **方法二：基于极点稳定性的后验检查** 拟合完成后，检查极点的位置。一个物理可实现的稳定系统，其极点应位于复平面的左半部分（连续时间系统）。如果拟合出的极点有正实部，可能意味着： 1. 模型阶数过高，拟合了噪声。 2. 数据质量差或频率范围不足。 3. 数值问题。 这时，你可能需要降低模型阶数，或检查数据预处理（如缩放）是否得当。 ### 5.2 处理多端口系统与矩阵值传递函数 我们目前处理的是单输入单输出（SISO）系统。对于多输入多输出（MIMO）系统，其传递函数是一个矩阵。一种直接的方法是**对矩阵的每个元素独立进行矢量拟合**。虽然这会忽略元素间的相关性，但对于许多工程应用来说已经足够，并且实现简单。 ```python def fit_mimo_tf(s: np.ndarray, H_matrix: np.ndarray, order_n: int) -> list: """ 对MIMO系统的传递函数矩阵进行拟合。 H_matrix的形状应为 (num_outputs, num_inputs, num_freq_points) 返回一个系数矩阵的列表，每个元素对应一个SISO传递函数的系数对。 """ num_outputs, num_inputs, num_freq = H_matrix.shape if len(s) != num_freq: raise ValueError("Frequency points dimension mismatch.") all_results = [] for i in range(num_outputs): row_results = [] for j in range(num_inputs): print(f"Fitting SISO transfer function H[{i},{j}]...") fitter = SanathananKoernerFitter(order_n=order_n) H_ij = H_matrix[i, j, :] result = fitter.fit_with_scaling(s, H_ij) row_results.append((fitter.a_coeffs, fitter.b_coeffs)) all_results.append(row_results) print("MIMO fitting complete.") return all_results # 示例：假设我们有一个2x2的系统 # H_matrix_example = np.random.randn(2, 2, len(s)) + 1j*np.random.randn(2, 2, len(s)) # 示例数据 # mimo_coeffs = fit_mimo_tf(s, H_matrix_example, order_n=2) ``` 对于大型MIMO系统，独立拟合可能计算量较大。更高级的方法（如**矩阵分式描述**或**公共极点拟合**）可以强制所有传递函数元素共享同一组极点，这能减少参数数量并保持模型的一致性，但实现起来复杂得多，通常需要求解非线性优化问题。 ### 5.3 代码封装与生产级考虑 为了在日常工作中方便使用，我们可以将整个流程封装成一个更完整的模块，并加入一些工程化特性： ```python class VectorFittingSK: """ 一个生产就绪的Sanathanan-Koerner矢量拟合类，包含错误处理、日志和更多选项。 """ def __init__(self, order_n: int, use_scaling: bool = True, fix_dc: bool = False): self.order_n = order_n self.use_scaling = use_scaling self.fix_dc = fix_dc # 是否强制拟合曲线在DC (s=0)处与数据匹配 self.fitter = SanathananKoernerFitter(order_n) def fit(self, frequencies: np.ndarray, response_data: np.ndarray, **kwargs): """ 主拟合接口。 frequencies: 频率数组 (Hz) response_data: 频率响应数据（复数） """ omega = 2 * np.pi * frequencies s = 1j * omega if self.fix_dc: # 找到DC点（频率为0或最接近0的点） dc_idx = np.argmin(np.abs(frequencies)) # 可以在此处添加约束，确保H_fit(s=0) = H_data(s=0) # 这通常通过向最小二乘问题添加一个额外的方程来实现 pass # 此处省略实现细节 if self.use_scaling: result = self.fitter.fit_with_scaling(s, response_data, **kwargs) else: result = self.fitter.fit(s, response_data, **kwargs) return result def export_to_spice(self, filename: str): """ 将拟合的有理函数导出为SPICE子电路或其他仿真器兼容的格式。 这是一个高级功能，需要根据目标仿真器的语法生成网表。 """ poles, zeros = self.fitter.get_poles_zeros() # 根据极点和零点生成等效电路（例如，使用RLC元件或受控源） # 此处仅为示意 with open(filename, 'w') as f: f.write("* Exported rational function from VectorFittingSK\n") f.write(f"* Poles: {poles}\n") f.write(f"* Zeros: {zeros}\n") # ... 具体的SPICE网表生成代码 print(f"Model exported to {filename}") # 使用示例 # vf = VectorFittingSK(order_n=3, use_scaling=True) # result = vf.fit(freq_array, measured_response_array, max_iter=20, tol=1e-10) # vf.export_to_spice('model.cir') ``` 在实际部署时，还需要考虑**异常处理**（如矩阵奇异、迭代不收敛）、**输入验证**（数据维度、频率单调性检查）以及**性能优化**（对于大量频率点，矩阵构建可以使用向量化操作加速）。对于超大规模问题，可能需要使用迭代求解器（如LSQR）而不是直接求解法方程。 最后，别忘了**可视化是调试和理解结果的最佳工具**。除了幅频/相频图，绘制**极点-零点图**能直观显示系统稳定性，绘制**时域阶跃响应**（通过逆拉普拉斯变换或数值积分）可以验证模型的因果性和动态行为是否合理。这些步骤能帮你建立起对拟合模型的信心，确保它不仅仅在频域匹配数据，在时域也能表现出符合物理直觉的行为。