feat: OVF formula 67
This commit is contained in:
mayge
174
core/orthonormal_basis.py
Normal file
174
core/orthonormal_basis.py
Normal file
@@ -0,0 +1,174 @@
|
||||
import numpy as np
|
||||
|
||||
# ------------------------------------------------------------
|
||||
# 按论文公式 (9)(10)(11) 生成 Muntz–Laguerre 正交有理基 (解析形式)：
|
||||
#
|
||||
# 给定稳定极点集合 {p_k} (Re(p_k)<0)。论文记法中使用 -a_k，其中 Re(a_k)>0。
|
||||
# 对应关系：p_k = -a_k ⇒ a_k = -p_k, Re(a_k)= -Re(p_k) >0
|
||||
#
|
||||
# 连续内积意义下（沿 jω 轴积分）这些 φ_k 解析正交。离散频率采样后数值上
|
||||
# 可能偏离，可再用加权 QR 做数值再正交（可选）。
|
||||
#
|
||||
# 公式在“稳定极点 p 表达”下的改写：
|
||||
# (原) 实极点: φ_p(s) = sqrt(2 Re(a_p)) / (s + a_p) * Π (s - a_i^*)/(s + a_i)
|
||||
# 变换 a_p = -p ⇒ Re(a_p)= -Re(p) = σ >0 且 (s + a_p) = (s - p)
|
||||
# 且乘积 (s - a_i^*)/(s + a_i) = (s + p_i^*)/(s - p_i)
|
||||
# ⇒ φ_p(s) = sqrt(-2 Re(p)) / (s - p) * Π_{i<p} (s + p_i^*)/(s - p_i)
|
||||
#
|
||||
# 复对 (p, p*)，取 imag(p)>0 的 p 作为首：
|
||||
# (原) φ_p = sqrt(2 Re(a_p)) (s - |a_p|)/[(s + a_p)(s + a_p^*)] * Π(...)
|
||||
# (原) φ_{p+1}= sqrt(2 Re(a_p)) (s + |a_p|)/[(s + a_p)(s + a_p^*)] * Π(...)
|
||||
# 代入 a_p=-p：
|
||||
# Re(a_p)= -Re(p)=σ>0, (s + a_p) = (s - p), (s + a_p^*)=(s - p^*)
|
||||
# |a_p| = |p|
|
||||
# 乘积同上 ⇒ Π_{i<p} (s + p_i^*)/(s - p_i)
|
||||
#
|
||||
# ⇒ φ_p(s) = sqrt(-2 Re(p)) (s - |p|)/[(s - p)(s - p^*)] * Π_{i<p} (s + p_i^*)/(s - p_i)
|
||||
# φ_{p^*}(s)= sqrt(-2 Re(p)) (s + |p|)/[(s - p)(s - p^*)] * Π_{i<p} (s + p_i^*)/(s - p_i)
|
||||
#
|
||||
# product 递推维护：prod_k = Π_{i≤k} (s + p_i^*)/(s - p_i)
|
||||
# 对复对要顺序乘两次（p 与 p*）。
|
||||
# ------------------------------------------------------------
|
||||
|
||||
|
||||
def generate_laguerre_basis(poles:list|np.ndarray, s: np.ndarray):
|
||||
basis = np.zeros((len(poles)+1,len(s)),dtype=complex)
|
||||
|
||||
product = np.ones(len(s),dtype=complex)
|
||||
|
||||
basis[0] = np.ones(len(s),dtype=complex) # φ_0 = 1
|
||||
i = 0
|
||||
while i < len(poles):
|
||||
if np.real(poles[i]) >= 0:
|
||||
raise ValueError(f"极点必须在左半平面: {poles[i]}")
|
||||
|
||||
# 复对首 (正虚部)
|
||||
if np.iscomplex(poles[i]) and np.imag(poles[i]) > 0:
|
||||
if i + 1 >= len(poles):
|
||||
raise ValueError("复极点缺少共轭")
|
||||
pc = poles[i + 1]
|
||||
if not np.isclose(pc, np.conj(poles[i])):
|
||||
raise ValueError("复极点未按 (p, p*) 顺序排列 (正虚部在前)")
|
||||
sigma = -np.real(poles[i]) # >0
|
||||
scale = np.sqrt(2 * sigma)
|
||||
r = np.abs(poles[i])
|
||||
denom = (s - poles[i]) * (s - pc)
|
||||
|
||||
# 两个基函数
|
||||
phi_p = scale * (s - r) / denom * product
|
||||
phi_pc = scale * (s + r) / denom * product
|
||||
|
||||
# product 先乘 (s + p^*)/(s - p)，再乘 (s + p)/(s - p^*)
|
||||
product = product * (s + pc) / (s - poles[i])
|
||||
product = product * (s + poles[i]) / (s - pc)
|
||||
|
||||
basis[i + 1] = phi_p
|
||||
basis[i + 2] = phi_pc
|
||||
i += 2
|
||||
continue
|
||||
|
||||
# 复对次 (负虚部) —— 应该被首元素处理，出现表示顺序错误
|
||||
if np.iscomplex(poles[i]) and np.imag(poles[i]) < 0:
|
||||
raise ValueError("检测到负虚部复极点但其共轭尚未处理，请将正虚部成员放在前面。")
|
||||
|
||||
# 实极点
|
||||
sigma = -np.real(poles[i])
|
||||
if sigma <= 0:
|
||||
raise ValueError("实极点实部应为负 (稳定)。")
|
||||
scale = np.sqrt(2 * sigma)
|
||||
phi = scale / (s - poles[i]) * product
|
||||
# 更新乘积
|
||||
product = product * (s + poles[i]) / (s - poles[i])
|
||||
i += 1
|
||||
basis[i + 1] = phi
|
||||
return basis
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# class MuntzLaguerreIterator:
|
||||
# def __init__(self, s: np.ndarray, stable_poles: list | np.ndarray):
|
||||
# """
|
||||
# s: 复频率数组 (Nf,), s = j 2π f
|
||||
# stable_poles: 稳定极点列表 (Re<0). 复共轭对要求正虚部在前 (p, p*).
|
||||
# """
|
||||
# self.s = np.asarray(s, dtype=complex)
|
||||
# self.poles = list(stable_poles)
|
||||
# self.N = len(self.poles)
|
||||
# self.k = 0
|
||||
# # 初始化乘积 Π_{i<p} (s + p_i^*)/(s - p_i)
|
||||
# self.product = np.ones_like(self.s, dtype=complex)
|
||||
|
||||
# def __iter__(self):
|
||||
# return self
|
||||
|
||||
# def __next__(self):
|
||||
# if self.k >= self.N:
|
||||
# raise StopIteration
|
||||
|
||||
# p = self.poles[self.k]
|
||||
# if np.real(p) >= 0:
|
||||
# raise ValueError(f"极点必须在左半平面: {p}")
|
||||
|
||||
# # 复对首 (正虚部)
|
||||
# if np.iscomplex(p) and np.imag(p) > 0:
|
||||
# if self.k + 1 >= self.N:
|
||||
# raise ValueError("复极点缺少共轭")
|
||||
# pc = self.poles[self.k + 1]
|
||||
# if not np.isclose(pc, np.conj(p)):
|
||||
# raise ValueError("复极点未按 (p, p*) 顺序排列 (正虚部在前)")
|
||||
# sigma = -np.real(p) # >0
|
||||
# scale = np.sqrt(2 * sigma)
|
||||
# r = np.abs(p)
|
||||
# denom = (self.s - p) * (self.s - pc)
|
||||
|
||||
# # 两个基函数
|
||||
# phi_p = scale * (self.s - r) / denom * self.product
|
||||
# phi_pc = scale * (self.s + r) / denom * self.product
|
||||
|
||||
# # product 先乘 (s + p^*)/(s - p)，再乘 (s + p)/(s - p^*)
|
||||
# self.product = self.product * (self.s + pc) / (self.s - p)
|
||||
# self.product = self.product * (self.s + p) / (self.s - pc)
|
||||
|
||||
# self.k += 2
|
||||
# return [phi_p, phi_pc]
|
||||
|
||||
# # 复对次 (负虚部) —— 应该被首元素处理，出现表示顺序错误
|
||||
# if np.iscomplex(p) and np.imag(p) < 0:
|
||||
# raise ValueError("检测到负虚部复极点但其共轭尚未处理，请将正虚部成员放在前面。")
|
||||
|
||||
# # 实极点
|
||||
# sigma = -np.real(p)
|
||||
# if sigma <= 0:
|
||||
# raise ValueError("实极点实部应为负 (稳定)。")
|
||||
# scale = np.sqrt(2 * sigma)
|
||||
# phi = scale / (self.s - p) * self.product
|
||||
# # 更新乘积
|
||||
# self.product = self.product * (self.s + p) / (self.s - p)
|
||||
|
||||
# self.k += 1
|
||||
# return [phi]
|
||||
|
||||
# def generate_muntz_laguerre_basis(s: np.ndarray, init_poles: list | np.ndarray):
|
||||
# """
|
||||
# 生成完整基函数列表: [φ_0=1, φ_1, φ_2, ...]
|
||||
# """
|
||||
# basis = [np.ones_like(s, dtype=complex)]
|
||||
# for block in MuntzLaguerreIterator(s, init_poles):
|
||||
# basis.extend(block)
|
||||
# return basis
|
||||
|
||||
# if __name__ == "__main__":
|
||||
# # 示例稳定极点 (复对正虚部在前)
|
||||
# stable_poles = [
|
||||
# -0.8e9,
|
||||
# -1.0e9 + 2.5e9j,
|
||||
# -1.0e9 - 2.5e9j,
|
||||
# -2.2e9
|
||||
# ]
|
||||
# freqs = np.linspace(1e8, 8e9, 400)
|
||||
# s = 1j * 2 * np.pi * freqs
|
||||
|
||||
# basis = generate_muntz_laguerre_basis(s, stable_poles)
|
||||
# print(f"生成 {len(basis)} 个基函数,{basis}")
|
||||
538
core/sk_iter.py
538
core/sk_iter.py
@@ -1,250 +1,326 @@
|
||||
import numpy as np
|
||||
from dataclasses import dataclass
|
||||
from typing import List, Sequence, Tuple, Optional
|
||||
from typing import List, Dict, Tuple
|
||||
|
||||
@dataclass
|
||||
class OPVFConfig:
|
||||
n_poles: int
|
||||
max_iter: int = 5
|
||||
tol: float = 1e-3
|
||||
add_constraint: bool = True
|
||||
real_split: bool = True
|
||||
lambda_reg: float = 0.0
|
||||
verbose: bool = True
|
||||
|
||||
# ---------- 参数正交基 ----------
|
||||
class ParamOrthoBasis:
|
||||
def __init__(self, params: np.ndarray, total_degree: int):
|
||||
# params: (Ns, d)
|
||||
self.mean = params.mean(0)
|
||||
self.std = params.std(0) + 1e-15
|
||||
self.exps = self._gen_exps(params.shape[1], total_degree)
|
||||
M = self._build_monomial_matrix(params) # (Ns,M)
|
||||
Q,R = np.linalg.qr(M)
|
||||
self.Q = Q
|
||||
self.R = R
|
||||
def _gen_exps(self, dim, D):
|
||||
exps=[]
|
||||
def rec(cur, i, rem):
|
||||
if i==dim:
|
||||
exps.append(tuple(cur)); return
|
||||
for k in range(rem+1):
|
||||
cur.append(k); rec(cur, i+1, rem-k); cur.pop()
|
||||
rec([],0,D)
|
||||
return exps
|
||||
def _build_monomial_matrix(self, params):
|
||||
X = (params - self.mean)/self.std
|
||||
Ns = X.shape[0]
|
||||
M = np.zeros((Ns, len(self.exps)))
|
||||
for idx,e in enumerate(self.exps):
|
||||
v = np.ones(Ns)
|
||||
for j,p in enumerate(e):
|
||||
if p: v *= X[:,j]**p
|
||||
M[:,idx]=v
|
||||
return M
|
||||
def eval(self, g: np.ndarray) -> np.ndarray:
|
||||
x = (g - self.mean)/self.std
|
||||
mono=[]
|
||||
for e in self.exps:
|
||||
val=1.0
|
||||
for j,p in enumerate(e):
|
||||
if p: val*= x[j]**p
|
||||
mono.append(val)
|
||||
mono=np.array(mono)
|
||||
# φ = mono @ R^{-1}
|
||||
return mono @ np.linalg.inv(self.R)
|
||||
|
||||
# ---------- 频率正交(有理)基 ----------
|
||||
class RationalFreqBasis:
|
||||
def __init__(self, freqs: np.ndarray, poles: np.ndarray):
|
||||
# 构造原始矩阵 G: 列0 常数 1, 后面 1/(s-a_k)
|
||||
s = 1j*2*np.pi*freqs
|
||||
cols=[np.ones_like(s)]
|
||||
for a in poles:
|
||||
cols.append(1.0/(s - a))
|
||||
G = np.vstack(cols).T # (Nf, K+1)
|
||||
# 加权可加 w_f，这里统一 1
|
||||
Q, R = np.linalg.qr(G)
|
||||
self.Q = Q # Φ(f) (Nf, K+1) 频率正交基取样
|
||||
self.R = R # 变换: G = Q R
|
||||
self.freqs = freqs
|
||||
def eval(self) -> np.ndarray:
|
||||
return self.Q # 直接返回已正交基取样 (Nf, K+1)
|
||||
|
||||
# ---------- 主模型 ----------
|
||||
class OrthonormalParametricVF:
|
||||
# ---------- 1. 生成初始极点（按论文 Re(-a_p) 小、Im 线性分布） ----------
|
||||
def generate_starting_poles(n_pairs: int, beta_min: float, beta_max: float, alpha_scale: float = 0.01):
|
||||
"""
|
||||
H(s, μ) = N(s,μ)/D(s,μ)
|
||||
N(s,μ)= Σ_k c_k(μ) φ_k(s); D(s,μ)= Σ_k ĉ_k(μ) φ_k(s)
|
||||
其中 φ_k(s) 是频率正交有理基; c_k(μ), ĉ_k(μ) 在参数正交基上展开:
|
||||
c_k(μ)= Σ_q C[k,q] ψ_q(μ); ĉ_k(μ)= Σ_q Ctilde[k,q] ψ_q(μ)
|
||||
仅生成复共轭对: p = -alpha + j beta, p*。
|
||||
n_pairs: 复对数量 (总极点数 = 2*n_pairs)
|
||||
beta_min,beta_max: 想要覆盖的虚部范围 (单位: rad/s)
|
||||
alpha_scale: alpha = alpha_scale * beta (文中 {α_p}=0.01{β_p})
|
||||
返回: list[complex] (正虚部先, 后跟共轭)
|
||||
"""
|
||||
def __init__(self, cfg: OPVFConfig, freq_basis: RationalFreqBasis, param_basis: ParamOrthoBasis):
|
||||
self.cfg = cfg
|
||||
self.fb = freq_basis
|
||||
self.pb = param_basis
|
||||
Kp1 = self.fb.Q.shape[1]
|
||||
Qp = self.pb.Q.shape[1]
|
||||
self.C = np.zeros((Kp1, Qp), dtype=complex) # 分子
|
||||
self.Ct = np.zeros((Kp1, Qp), dtype=complex) # 分母
|
||||
# 初始化分母：ĉ_0 ≈ 1，其余 0
|
||||
self.Ct[0,0] = 1.0
|
||||
betas = 2*np.pi*np.linspace(beta_min, beta_max, n_pairs)
|
||||
poles = []
|
||||
for b in betas:
|
||||
alpha = alpha_scale * b
|
||||
p = -alpha + 1j * b
|
||||
poles += [p, np.conj(p)]
|
||||
print(f"生成 {len(poles)} 个初始极点 (复对) {poles}]")
|
||||
return poles
|
||||
|
||||
def _assemble_phi_param(self, params: np.ndarray) -> np.ndarray:
|
||||
# 返回 (Ns, Qp)
|
||||
return self.pb.Q
|
||||
# ---------- 2. Muntz–Laguerre 频率基 (与你现有 orthonormal_basis 中一致的核心) ----------
|
||||
def muntz_laguerre_basis(freqs_hz: np.ndarray, stable_poles: List[complex]) -> np.ndarray:
|
||||
"""
|
||||
返回 Φ_ml (Nf, P) 列顺序: φ_0=1, 后续依次 (单极点 / 复对两列)
|
||||
"""
|
||||
s = 1j * 2 * np.pi * freqs_hz
|
||||
basis = [np.ones_like(s, dtype=complex)]
|
||||
product = np.ones_like(s, dtype=complex)
|
||||
k = 0
|
||||
N = len(stable_poles)
|
||||
while k < N:
|
||||
p = stable_poles[k]
|
||||
if np.real(p) >= 0:
|
||||
raise ValueError("极点需在左半平面")
|
||||
if np.imag(p) > 0: # 复对首
|
||||
if k + 1 >= N or not np.isclose(stable_poles[k+1], np.conj(p)):
|
||||
raise ValueError("复对未配对")
|
||||
pc = stable_poles[k+1]
|
||||
sigma = -np.real(p)
|
||||
scale = np.sqrt(2 * sigma)
|
||||
r = np.abs(p)
|
||||
denom = (s - p) * (s - pc)
|
||||
phi_p = scale * (s - r) / denom * product
|
||||
phi_pc = scale * (s + r) / denom * product
|
||||
basis.extend([phi_p, phi_pc])
|
||||
# update product
|
||||
product = product * (s + pc)/(s - p) * (s + p)/(s - pc)
|
||||
k += 2
|
||||
elif np.imag(p) < 0:
|
||||
raise ValueError("负虚部复极点必须跟在正虚部之后")
|
||||
else: # 实极点
|
||||
sigma = -np.real(p)
|
||||
scale = np.sqrt(2 * sigma)
|
||||
phi = scale / (s - p) * product
|
||||
basis.append(phi)
|
||||
product = product * (s + p)/(s - p)
|
||||
k += 1
|
||||
return np.column_stack(basis) # (Nf, P)
|
||||
|
||||
def fit(self, H_samples: List[np.ndarray], params: np.ndarray):
|
||||
"""
|
||||
H_samples: list 长度 Ns, 每个 (Nf,) 复
|
||||
params: (Ns,d)
|
||||
"""
|
||||
Φf = self.fb.eval() # (Nf, K+1)
|
||||
Ns = len(H_samples)
|
||||
Nf, Kp1 = Φf.shape
|
||||
Φμ = self._assemble_phi_param(params) # (Ns, Qp)
|
||||
Qp = Φμ.shape[1]
|
||||
# ---------- 3. 原始部分分式基 ----------
|
||||
def raw_partial_fraction_basis(freqs_hz: np.ndarray, stable_poles: List[complex]) -> np.ndarray:
|
||||
s = 1j * 2 * np.pi * freqs_hz
|
||||
cols = [np.ones_like(s, dtype=complex)]
|
||||
for p in stable_poles:
|
||||
cols.append(1.0 / (s - p))
|
||||
return np.column_stack(cols) # (Nf, P)
|
||||
|
||||
# t=0: 只解 C (式3, D=1)
|
||||
self._solve_t0(H_samples, Φf, Φμ)
|
||||
# ---------- 4. 变换矩阵 Φ_ml T = G_raw ----------
|
||||
def compute_transform_matrix(Phi_ml: np.ndarray, G_raw: np.ndarray, weights: np.ndarray | None = None) -> np.ndarray:
|
||||
# 加权最小二乘 T = (Φ^H W Φ)^{-1} Φ^H W G
|
||||
if weights is None:
|
||||
WPhi = Phi_ml
|
||||
WG = G_raw
|
||||
M = Phi_ml.conj().T @ WPhi
|
||||
RHS = Phi_ml.conj().T @ WG
|
||||
else:
|
||||
w = weights[:, None]
|
||||
M = Phi_ml.conj().T @ (w * Phi_ml)
|
||||
RHS = Phi_ml.conj().T @ (w * G_raw)
|
||||
T = np.linalg.solve(M, RHS)
|
||||
return T # (P,P), Φ T = G
|
||||
|
||||
D_prev = self._eval_D(Φf, Φμ) # (Nf, Ns)
|
||||
# ---------- 5. 参数多项式正交基 ----------
|
||||
def generate_param_basis(param_matrix: np.ndarray, total_degree: int):
|
||||
"""
|
||||
param_matrix: (Ns, d)
|
||||
返回:
|
||||
Qμ: (Ns, V) 正交列
|
||||
Rμ: (V, V) 上三角
|
||||
exps: 多指数列表
|
||||
"""
|
||||
X = param_matrix
|
||||
Ns, d = X.shape
|
||||
# 生成多指数
|
||||
exps=[]
|
||||
def rec(cur,i,rem):
|
||||
if i==d:
|
||||
exps.append(tuple(cur)); return
|
||||
for k in range(rem+1):
|
||||
cur.append(k); rec(cur,i+1,rem-k); cur.pop()
|
||||
rec([],0,total_degree)
|
||||
V = len(exps)
|
||||
M = np.zeros((Ns, V))
|
||||
for idx,e in enumerate(exps):
|
||||
v = np.ones(Ns)
|
||||
for j,p in enumerate(e):
|
||||
if p: v *= X[:,j]**p
|
||||
M[:,idx]=v
|
||||
Q,R = np.linalg.qr(M)
|
||||
return Q,R,exps
|
||||
|
||||
for it in range(1, self.cfg.max_iter+1):
|
||||
A, b = self._build_iter_system(H_samples, Φf, Φμ, D_prev)
|
||||
if self.cfg.lambda_reg>0:
|
||||
lam = self.cfg.lambda_reg
|
||||
A = np.vstack([A, lam*np.eye(A.shape[1])])
|
||||
b = np.concatenate([b, np.zeros(A.shape[1])])
|
||||
x, *_ = np.linalg.lstsq(A, b, rcond=None)
|
||||
self._unpack_iter_solution(x, Kp1, Qp)
|
||||
D_new = self._eval_D(Φf, Φμ)
|
||||
rel = np.max(np.abs(D_new-D_prev)/(np.abs(D_prev)+1e-12))
|
||||
if self.cfg.verbose:
|
||||
print(f"[OPVF-Iter {it}] rel_change={rel:.3e}")
|
||||
D_prev = D_new
|
||||
if rel < self.cfg.tol:
|
||||
if self.cfg.verbose:
|
||||
print(f"[OPVF] converged at {it}")
|
||||
break
|
||||
return self
|
||||
# ---------- 6. t=0 构建设计矩阵并解 C ----------
|
||||
def fit_t0(H_list: List[np.ndarray], Phi_f: np.ndarray, Qμ: np.ndarray) -> np.ndarray:
|
||||
"""
|
||||
H_list: 长度 Ns, 每项 (Nf,) 复数
|
||||
Phi_f: (Nf, P)
|
||||
Qμ: (Ns, V)
|
||||
返回: C (P,V)
|
||||
"""
|
||||
Ns = len(H_list); Nf,P = Phi_f.shape; V = Qμ.shape[1]
|
||||
cols = P*V
|
||||
A = np.zeros((Ns*Nf, cols), dtype=complex)
|
||||
b = np.zeros(Ns*Nf, dtype=complex)
|
||||
r=0
|
||||
for n,Hn in enumerate(H_list):
|
||||
blk = np.einsum('fp,v->fpv', Phi_f, Qμ[n]).reshape(Nf, cols)
|
||||
A[r:r+Nf,:]=blk
|
||||
b[r:r+Nf]=Hn
|
||||
r+=Nf
|
||||
x, *_ = np.linalg.lstsq(A, b, rcond=None)
|
||||
C = x.reshape(P, V)
|
||||
return C
|
||||
|
||||
def _solve_t0(self, H_samples, Φf, Φμ):
|
||||
Ns = len(H_samples); Nf,Kp1 = Φf.shape; Qp = Φμ.shape[1]
|
||||
# 设计矩阵行数 Ns*Nf；未知数 Kp1*Qp
|
||||
cols = Kp1*Qp
|
||||
A = np.zeros((Ns*Nf, cols), dtype=complex)
|
||||
b = np.zeros(Ns*Nf, dtype=complex)
|
||||
r = 0
|
||||
for n,Hn in enumerate(H_samples):
|
||||
phiμ = Φμ[n] # (Qp,)
|
||||
# Φf (Nf,Kp1) 外积 phiμ -> (Nf, Kp1*Qp)
|
||||
blk = np.einsum('fk,q->fkq', Φf, phiμ).reshape(Nf, cols)
|
||||
A[r:r+Nf,:] = blk
|
||||
b[r:r+Nf] = Hn
|
||||
r += Nf
|
||||
x, *_ = np.linalg.lstsq(A, b, rcond=None)
|
||||
self.C = x.reshape(Kp1, Qp)
|
||||
# 分母保持初始 (ĉ_0=1)
|
||||
|
||||
def _build_iter_system(self, H_samples, Φf, Φμ, D_prev):
|
||||
Ns = len(H_samples); Nf,Kp1 = Φf.shape; Qp=Φμ.shape[1]
|
||||
# 未知：C (Kp1*Qp) + Ct (Kp1*Qp) 但固定 Ct[0,0]=1，可去掉该变量
|
||||
mask_fix = np.zeros((Kp1,Qp), dtype=bool)
|
||||
mask_fix[0,0]=True
|
||||
idx_map={}
|
||||
# ---------- 7. t=1 (首轮 SK) 解 C, Ct ----------
|
||||
def fit_t1_SK(H_list: List[np.ndarray], Phi_f: np.ndarray, Qμ: np.ndarray,
|
||||
C_prev: np.ndarray, max_iter: int =1, tol: float =1e-3):
|
||||
"""
|
||||
只做一轮 (或少量) SK，初始 D=1。
|
||||
返回: C_new, Ct_new
|
||||
"""
|
||||
Ns = len(H_list); Nf,P = Phi_f.shape; V = Qμ.shape[1]
|
||||
# 初始化分母系数 Ct: Ct[0,0]=1
|
||||
Ct = np.zeros((P,V), dtype=complex)
|
||||
Ct[0,0]=1.0
|
||||
C = C_prev.copy()
|
||||
for it in range(max_iter):
|
||||
# 构建线性系统: (N - D*H)=0
|
||||
# 未知顺序: C(:), Ct(:) 去掉固定 Ct[0,0]
|
||||
mask_fix = np.zeros((P,V), dtype=bool); mask_fix[0,0]=True
|
||||
col_map={}
|
||||
col=0
|
||||
for i in range(Kp1):
|
||||
for j in range(Qp):
|
||||
idx_map[('C',i,j)] = col; col+=1
|
||||
for i in range(Kp1):
|
||||
for j in range(Qp):
|
||||
for i in range(P):
|
||||
for j in range(V):
|
||||
col_map[('C',i,j)]=col; col+=1
|
||||
for i in range(P):
|
||||
for j in range(V):
|
||||
if mask_fix[i,j]: continue
|
||||
idx_map[('Ct',i,j)] = col; col+=1
|
||||
cols = col
|
||||
rows = Ns*Nf
|
||||
A = np.zeros((rows, cols), dtype=complex)
|
||||
b = np.zeros(rows, dtype=complex)
|
||||
col_map[('Ct',i,j)]=col; col+=1
|
||||
A = np.zeros((Ns*Nf, col), dtype=complex)
|
||||
b = np.zeros(Ns*Nf, dtype=complex)
|
||||
r=0
|
||||
for n,Hn in enumerate(H_samples):
|
||||
phiμ = Φμ[n]
|
||||
Dp = D_prev[:,n] # (Nf,)
|
||||
invDp = 1.0/(Dp + 1e-12)
|
||||
# 分子块: (Φf φμ^T) / D_prev
|
||||
Num_blk = np.einsum('fk,q->fkq', Φf, phiμ).reshape(Nf, Kp1*Qp)
|
||||
Num_blk = (Num_blk.T * invDp).T
|
||||
# 评价当前 D_prev= Σ Ct φ ψ
|
||||
for n,Hn in enumerate(H_list):
|
||||
# 分子块 (Φ_f ⊗ ψ_n)
|
||||
blk = np.einsum('fp,v->fpv', Phi_f, Qμ[n]).reshape(Nf, P*V)
|
||||
# 填 C 部分
|
||||
for i in range(Kp1):
|
||||
for j in range(Qp):
|
||||
A[r:r+Nf, idx_map[('C',i,j)]] = Num_blk[:, i*Qp + j]
|
||||
# 分母块: -(H * Φf φμ^T)/D_prev
|
||||
Den_blk = (Num_blk.T * Hn).T * (-1.0)
|
||||
for i in range(Kp1):
|
||||
for j in range(Qp):
|
||||
for i in range(P):
|
||||
for j in range(V):
|
||||
A[r:r+Nf, col_map[('C',i,j)]] = blk[:, i*V + j]
|
||||
# 分母块: - H * blk
|
||||
for i in range(P):
|
||||
for j in range(V):
|
||||
if mask_fix[i,j]: continue
|
||||
A[r:r+Nf, idx_map[('Ct',i,j)]] = Den_blk[:, i*Qp + j]
|
||||
# 右端 0
|
||||
r += Nf
|
||||
|
||||
# 约束行 (式5)(8) 可选：Re( Σ_k D^{(t)}/D^{(t-1)} ) = K+1
|
||||
if self.cfg.add_constraint:
|
||||
row_c = np.zeros(cols, dtype=complex)
|
||||
# D^{(t)} ≈ Σ Ct_k φ_k(s) ; 用代表样本 n=0, 对所有频率平均
|
||||
n0=0
|
||||
phiμ0=Φμ[n0]
|
||||
mean_phiμ = phiμ0 # 可换平均
|
||||
# φ_k(s) 平均
|
||||
mean_phif = Φf.mean(0) # (Kp1,)
|
||||
for i in range(Kp1):
|
||||
for j in range(Qp):
|
||||
if mask_fix[i,j]: continue
|
||||
row_c[idx_map[('Ct',i,j)]] = mean_phif[i]*mean_phiμ[j]
|
||||
rhs = (Kp1) # K+1
|
||||
A = np.vstack([A, np.real(row_c) if self.cfg.real_split else row_c])
|
||||
b = np.concatenate([b, [rhs]])
|
||||
|
||||
# 拆实虚
|
||||
if self.cfg.real_split:
|
||||
A_real = np.vstack([np.real(A), np.imag(A)])
|
||||
b_real = np.concatenate([np.real(b), np.imag(b)])
|
||||
else:
|
||||
A_real, b_real = A, b
|
||||
return A_real, b_real
|
||||
|
||||
def _unpack_iter_solution(self, x, Kp1, Qp):
|
||||
# 重新填回 C, Ct (保持 Ct[0,0]=1)
|
||||
# 构建与 _build_iter_system 相同的 idx_map
|
||||
mask_fix = np.zeros((Kp1,Qp), dtype=bool); mask_fix[0,0]=True
|
||||
idx_C = Kp1*Qp
|
||||
# 注意：我们当时 C 索引从 0..(Kp1*Qp-1)
|
||||
self.C = x[:idx_C].reshape(Kp1, Qp)
|
||||
Ct_new = self.Ct.copy()
|
||||
pos = idx_C
|
||||
for i in range(Kp1):
|
||||
for j in range(Qp):
|
||||
A[r:r+Nf, col_map[('Ct',i,j)]] = - Hn * blk[:, i*V + j]
|
||||
r+=Nf
|
||||
# 求解
|
||||
x, *_ = np.linalg.lstsq(A, b, rcond=None)
|
||||
# 拆回
|
||||
idx=0
|
||||
for i in range(P):
|
||||
for j in range(V):
|
||||
C[i,j]=x[idx]; idx+=1
|
||||
Ct_new = Ct.copy()
|
||||
for i in range(P):
|
||||
for j in range(V):
|
||||
if mask_fix[i,j]: continue
|
||||
Ct_new[i,j]=x[pos]; pos+=1
|
||||
self.Ct = Ct_new
|
||||
Ct_new[i,j]=x[idx]; idx+=1
|
||||
# 收敛性简单检查：分母变化
|
||||
diff = np.max(np.abs(Ct_new - Ct))
|
||||
Ct = Ct_new
|
||||
if diff < tol:
|
||||
break
|
||||
return C, Ct
|
||||
|
||||
def _eval_C_mu(self, phiμ):
|
||||
# 返回 c_k(μ) (Kp1,)
|
||||
return self.C @ phiμ
|
||||
def _eval_Ct_mu(self, phiμ):
|
||||
return self.Ct @ phiμ
|
||||
def _eval_D(self, Φf, Φμ):
|
||||
# D(f, sample)= Σ_k ĉ_k(μ_n) φ_k(f)
|
||||
Kp1 = Φf.shape[1]
|
||||
Ns = Φμ.shape[0]
|
||||
D = np.zeros((Φf.shape[0], Ns), dtype=complex)
|
||||
for n in range(Ns):
|
||||
ct_mu = self._eval_Ct_mu(Φμ[n])
|
||||
D[:,n] = Φf @ ct_mu
|
||||
return D
|
||||
# ---------- 8. 假极点检测 (在 raw 基上) ----------
|
||||
def detect_spurious_poles(C_ml: np.ndarray, Ct_ml: np.ndarray, T: np.ndarray,
|
||||
Qμ: np.ndarray,
|
||||
tau_cancel=1e-2, tau_small=1e-3, eps=1e-14):
|
||||
"""
|
||||
C_ml, Ct_ml: (P,V) (Muntz-Laguerre 基)
|
||||
T: Φ_ml T = G_raw
|
||||
Qμ: (Ns,V)
|
||||
"""
|
||||
# raw 系数: c_raw = T^{-1} c_ml
|
||||
Tinv = np.linalg.inv(T)
|
||||
C_raw = Tinv @ C_ml
|
||||
Ct_raw = Tinv @ Ct_ml
|
||||
P,V = C_ml.shape
|
||||
Ns = Qμ.shape[0]
|
||||
r_vals = C_raw @ Qμ.T
|
||||
q_vals = Ct_raw @ Qμ.T
|
||||
metrics={}
|
||||
scales=[]
|
||||
for k in range(1,P): # 跳过常数列
|
||||
rv = r_vals[k]; qv = q_vals[k]
|
||||
diff = np.max(np.abs(rv - qv))
|
||||
scale = np.max([np.max(np.abs(rv)), np.max(np.abs(qv)), eps])
|
||||
eta = diff / scale
|
||||
metrics[k] = {"diff": diff, "scale": scale, "eta": eta}
|
||||
scales.append(scale)
|
||||
S_max = max(scales) if scales else 1.0
|
||||
cancel_spurious = [k for k in range(1,P)
|
||||
if metrics[k]["eta"] < tau_cancel and metrics[k]["scale"] > tau_small * S_max]
|
||||
small_spurious = [k for k in range(1,P)
|
||||
if metrics[k]["scale"] <= tau_small * S_max]
|
||||
return {
|
||||
"cancel_spurious": cancel_spurious,
|
||||
"small_spurious": small_spurious,
|
||||
"metrics": metrics,
|
||||
"S_max": S_max,
|
||||
"C_raw": C_raw,
|
||||
"Ct_raw": Ct_raw
|
||||
}
|
||||
|
||||
def evaluate(self, freqs: np.ndarray, g: np.ndarray):
|
||||
assert np.allclose(freqs, self.fb.freqs), "频率需在训练网格上 (示例简化)"
|
||||
Φf = self.fb.Q # (Nf,K+1)
|
||||
phiμ = self.pb.eval(g) # (Qp,)
|
||||
num = Φf @ (self.C @ phiμ)
|
||||
den = Φf @ (self.Ct @ phiμ)
|
||||
return num / (den + 1e-15)
|
||||
# ---------- 9. 统一入口 ----------
|
||||
@dataclass
|
||||
class OPVFResult:
|
||||
C_ml: np.ndarray
|
||||
Ct_ml: np.ndarray
|
||||
C_raw: np.ndarray
|
||||
Ct_raw: np.ndarray
|
||||
T: np.ndarray
|
||||
Qμ: np.ndarray
|
||||
Rμ: np.ndarray
|
||||
exps: List[Tuple[int]]
|
||||
poles: List[complex]
|
||||
spurious: Dict
|
||||
|
||||
def opvf_from_H(H_list: List[np.ndarray],
|
||||
freqs_hz: np.ndarray,
|
||||
param_matrix: np.ndarray,
|
||||
total_degree: int,
|
||||
poles: List[complex],
|
||||
do_t1: bool = True) -> OPVFResult:
|
||||
"""
|
||||
H_list: 长度 Ns; 每项 shape (Nf,)
|
||||
freqs_hz: (Nf,)
|
||||
param_matrix: (Ns,d) 归一化前参数值 (内部不自动标准化以保持可控)
|
||||
poles: 初始稳定极点列表
|
||||
"""
|
||||
# 频率基
|
||||
Phi_ml = muntz_laguerre_basis(freqs_hz, poles) # (Nf,P)
|
||||
G_raw = raw_partial_fraction_basis(freqs_hz, poles) # (Nf,P)
|
||||
# 简单权 (ω 权): w = Δf (因 (1/2π)∫ φφ* dω => ∑ Δf)
|
||||
w = _trap_weights(freqs_hz)
|
||||
T = compute_transform_matrix(Phi_ml, G_raw, w)
|
||||
|
||||
# 参数基
|
||||
Qμ, Rμ, exps = generate_param_basis(param_matrix, total_degree)
|
||||
|
||||
# t=0
|
||||
C_ml = fit_t0(H_list, Phi_ml, Qμ)
|
||||
|
||||
if do_t1:
|
||||
C_ml, Ct_ml = fit_t1_SK(H_list, Phi_ml, Qμ, C_ml, max_iter=1)
|
||||
else:
|
||||
# D=1
|
||||
P,V = C_ml.shape
|
||||
Ct_ml = np.zeros((P,V), dtype=complex)
|
||||
Ct_ml[0,0]=1.0
|
||||
|
||||
# 假极点检测
|
||||
spurious = detect_spurious_poles(C_ml, Ct_ml, T, Qμ)
|
||||
|
||||
return OPVFResult(
|
||||
C_ml=C_ml,
|
||||
Ct_ml=Ct_ml,
|
||||
C_raw=spurious["C_raw"],
|
||||
Ct_raw=spurious["Ct_raw"],
|
||||
T=T,
|
||||
Qμ=Qμ,
|
||||
Rμ=Rμ,
|
||||
exps=exps,
|
||||
poles=poles,
|
||||
spurious=spurious
|
||||
)
|
||||
|
||||
# ---------- 辅助: 梯形权 ----------
|
||||
def _trap_weights(f: np.ndarray):
|
||||
if len(f)==1: return np.ones(1)
|
||||
df = np.diff(f)
|
||||
w = np.zeros_like(f)
|
||||
w[0]=0.5*df[0]; w[-1]=0.5*df[-1]
|
||||
if len(f)>2:
|
||||
w[1:-1]=0.5*(df[:-1]+df[1:])
|
||||
return w
|
||||
|
||||
# ---------- 简单测试占位 ----------
|
||||
if __name__ == "__main__":
|
||||
# 虚构数据: 2 个参数样本 (Ns=2), 频率 200 点
|
||||
freqs = np.linspace(1e8, 5e9, 200)
|
||||
# 真正模型 (示例): H = Σ_k R_k/(s - p_k)
|
||||
true_poles = [-0.5e3 + 1.2e9j, -0.5e3 - 1.2e9j]
|
||||
s = 1j*2*np.pi*freqs
|
||||
def synth(Rs):
|
||||
return Rs[0]/(s - true_poles[0]) + Rs[1]/(s - true_poles[1])
|
||||
H_list = [synth([0.8+0.1j, 1.2-0.2j]), synth([0.9+0.05j, 1.1+0.1j])]
|
||||
params = np.array([[0.0],[1.0]]) # 1 维参数
|
||||
# 给一个冗余极点集合 (含真实 + 额外)
|
||||
start_poles = generate_starting_poles(n_pairs=10, beta_min=5e8, beta_max=2.0e9, alpha_scale=0.000001)
|
||||
res = opvf_from_H(H_list, freqs, params, total_degree=1, poles=start_poles)
|
||||
print("C_ml shape:", res.C_ml.shape)
|
||||
print("假极点(cancel):", res.spurious["cancel_spurious"])
|
||||
print("假极点(small):", res.spurious["small_spurious"])
|
||||
@@ -1,3 +1,5 @@
|
||||
## 留数定理
|
||||
|
||||
### 1. 正交归一化条件
|
||||
|
||||
首先给出了一组基函数 \(\phi_m(s)\)，要求它们满足正交归一化条件：
|
||||
|
||||
13
docs/有理函数正交化.md
Normal file
13
docs/有理函数正交化.md
Normal file
@@ -0,0 +1,13 @@
|
||||
$$
|
||||
\begin{align}
|
||||
&N(s) = \sum^N_{i=0}n_is^i=\xi_1\prod^N_{i=0}(s-p_{ni})\\
|
||||
&D(s) = \sum^M_{i=0}d_is^i=\xi_2\prod^M_{i=0}(s-p_{mi})\\
|
||||
&H(s) = \xi\frac{\prod^N_{i=0}(s-p_{ni})}{\prod^M_{i=0}(s-p_{mi})}\\
|
||||
&let \space \phi_x=\frac{1}{s-p_x}\\
|
||||
&H(s) = \frac{\sum_{i=0}^Nc_{ni}\phi_i(s)}{\sum_{j=0}^Mc_{mj}\phi_j(s)}\\
|
||||
&当M=N时，通分之后得到\\
|
||||
&H(s) = \frac{\frac{\sum_{i=0}^Nc_i\prod_{a=0,a\not=i}(s-p_a)}{\prod^N_{i=0}(s-p_i)}}{\frac{\sum_{j=0}^Nc_j\prod_{b=0,b\not=j}(s-p_b)}{\prod^N_{j=0}(s-p_j)}}\\
|
||||
&此时由于H(s)分母的分子连乘部分存在缺项，所以p_b为假极点
|
||||
\end{align}
|
||||
$$
|
||||
|
||||
2
env
Normal file
2
env
Normal file
@@ -0,0 +1,2 @@
|
||||
export PYTHONPATH=$(pwd)/core:$PYTHONPATH
|
||||
export PYTHONPATH=$(pwd)/models:$PYTHONPATH
|
||||
128
main.py
128
main.py
@@ -1,15 +1,117 @@
|
||||
from models.capa import Capa
|
||||
W = []
|
||||
L = []
|
||||
i = 15.52
|
||||
|
||||
while i <= 100:
|
||||
W.append(i)
|
||||
L.append(i)
|
||||
i = int(i*1.05*100 +0.5) / 100.0
|
||||
|
||||
capa = Capa()
|
||||
capa.sweep(W,L)
|
||||
print(capa.network(result_dir=capa.results[0].result_dir,id=capa.results[0].id))
|
||||
from core.orthonormal_basis import generate_laguerre_basis
|
||||
from core.sk_iter import generate_starting_poles
|
||||
import numpy as np
|
||||
import skrf as rf
|
||||
import matplotlib.pyplot as plt
|
||||
import sympy as sp
|
||||
from scipy.linalg import null_space
|
||||
|
||||
|
||||
network = rf.Network("/tmp/paramer/simulation/4000/4000.s2p")
|
||||
freqs = network.f
|
||||
s = freqs * 2j * np.pi
|
||||
vf = rf.VectorFitting(network)
|
||||
vf.vector_fit(n_poles_real=2, n_poles_cmplx=1,parameter_type="y")
|
||||
|
||||
poles = vf.poles
|
||||
residues = vf.residues
|
||||
y = vf.network.y
|
||||
|
||||
fig, ax = plt.subplots(2, 2)
|
||||
fig.set_size_inches(12, 8)
|
||||
vf.plot("mag",0,0,freqs,ax=ax[0][0],parameter="y")
|
||||
rms_error11 = vf.get_rms_error(0,0,"y")
|
||||
print("rms_error11",rms_error11)
|
||||
vf.plot("mag",1,0,freqs,ax=ax[1][0],parameter="y")
|
||||
rms_error21 = vf.get_rms_error(1,0,"y")
|
||||
print("rms_error21",rms_error21)
|
||||
vf.plot("mag",0,1,freqs,ax=ax[0][1],parameter="y")
|
||||
rms_error12 = vf.get_rms_error(0,1,"y")
|
||||
print("rms_error12",rms_error12)
|
||||
vf.plot("mag",1,1,freqs,ax=ax[1][1],parameter="y")
|
||||
rms_error22 = vf.get_rms_error(1,1,"y")
|
||||
print("rms_error22",rms_error22)
|
||||
fig.tight_layout()
|
||||
plt.show()
|
||||
plt.savefig(f"img.png")
|
||||
|
||||
diag_values = [y[i][0][0] for i in range(len(y))]
|
||||
H = np.diag(diag_values)
|
||||
P = 4
|
||||
start_poles = generate_starting_poles(P,beta_min=1e8,beta_max=freqs[-1])
|
||||
basis = generate_laguerre_basis(start_poles,s).T
|
||||
print("start_poles",start_poles)
|
||||
print("basis",basis)
|
||||
|
||||
# first step iteration
|
||||
# A*x = b
|
||||
|
||||
A11 = np.real(basis)
|
||||
print("A11 shape:",A11.shape)
|
||||
A12 = np.real(- H @ basis[:,1:])
|
||||
print("A12 shape:",A12.shape)
|
||||
# print("A11",A11)
|
||||
A21 = np.imag(basis)
|
||||
print("A21 shape:",A21.shape)
|
||||
A22 = np.imag(- H @ basis[:,1:])
|
||||
print("A22 shape:",A22.shape)
|
||||
A1 = np.hstack([A11,A12])
|
||||
A2 = np.hstack([A21,A22])
|
||||
A = np.vstack([A1,A2])
|
||||
print ("A shape:",A.shape)
|
||||
|
||||
b1 = np.real(H @ basis[:,0])
|
||||
b2 = np.imag(H @ basis[:,0])
|
||||
b = np.hstack([b1,b2])
|
||||
Q, R = np.linalg.qr(A, mode='reduced')
|
||||
print("Q :",Q)
|
||||
print("R :",R)
|
||||
print("b shape:",b.shape)
|
||||
# x = np.linalg.solve(R, Q.T @ b)
|
||||
x_star, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
|
||||
|
||||
# print("x_qr",x_star)
|
||||
print("residuals",residuals)
|
||||
# print("rank",rank)
|
||||
# print("s",s)
|
||||
|
||||
x_ne = np.linalg.inv(A.T @ A) @ A.T @ b
|
||||
print("x_ne",x_ne)
|
||||
|
||||
# sk iteration
|
||||
target_vec = np.array([1+0j] + list(x_ne[len(x_ne)//2+1:]))
|
||||
D_t = basis @ target_vec
|
||||
print("D_t",D_t)
|
||||
K = 25
|
||||
|
||||
for i in range(K):
|
||||
print(f"Iteration {i+1}/{K}")
|
||||
A11 = np.real(basis / D_t[:, np.newaxis])
|
||||
A12 = np.real(- H @ basis[:,1:] / D_t[:, np.newaxis])
|
||||
A21 = np.imag(basis / D_t[:, np.newaxis])
|
||||
A22 = np.imag(- H @ basis[:,1:] / D_t[:, np.newaxis])
|
||||
A1 = np.hstack([A11,A12])
|
||||
A2 = np.hstack([A21,A22])
|
||||
A = np.vstack([A1,A2])
|
||||
b1 = np.real(H @ basis[:,0])
|
||||
b2 = np.imag(H @ basis[:,0])
|
||||
b = np.hstack([b1,b2])
|
||||
x_star, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
|
||||
# print("x_lstsq",x_star)
|
||||
# print("residuals",residuals)
|
||||
# print("rank",rank)
|
||||
# print("s",s)
|
||||
x_ne = np.linalg.inv(A.T @ A) @ A.T @ b
|
||||
# print("x_ne",x_ne)
|
||||
target_vec = np.array([1+0j] + list(x_ne[len(x_ne)//2+1:]))
|
||||
D_t_pre = D_t
|
||||
D_t = basis @ target_vec
|
||||
print("D_t/D_t_pre",np.linalg.norm(D_t)/np.linalg.norm(D_t_pre))
|
||||
# print("D_t",D_t)
|
||||
n_target_vec = np.array(list(x_ne[:len(x_ne)//2+1]))
|
||||
N_t = basis @ n_target_vec
|
||||
# print("H = N_t / D_t",np.linalg.norm(H - N_t / D_t)/np.linalg.norm(H))
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -1,202 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
# ------------------------------------------------------------
|
||||
# 按论文公式 (9)(10)(11) 生成 Muntz–Laguerre 正交有理基 (解析形式)：
|
||||
#
|
||||
# 给定稳定极点集合 {p_k} (Re(p_k)<0)。论文记法中使用 -a_k，其中 Re(a_k)>0。
|
||||
# 对应关系：p_k = -a_k ⇒ a_k = -p_k, Re(a_k)= -Re(p_k) >0
|
||||
#
|
||||
# 连续内积意义下（沿 jω 轴积分）这些 φ_k 解析正交。离散频率采样后数值上
|
||||
# 可能偏离，可再用加权 QR 做数值再正交（可选）。
|
||||
#
|
||||
# 公式在“稳定极点 p 表达”下的改写：
|
||||
# (原) 实极点: φ_p(s) = sqrt(2 Re(a_p)) / (s + a_p) * Π (s - a_i^*)/(s + a_i)
|
||||
# 变换 a_p = -p ⇒ Re(a_p)= -Re(p) = σ >0 且 (s + a_p) = (s - p)
|
||||
# 且乘积 (s - a_i^*)/(s + a_i) = (s + p_i^*)/(s - p_i)
|
||||
# ⇒ φ_p(s) = sqrt(-2 Re(p)) / (s - p) * Π_{i<p} (s + p_i^*)/(s - p_i)
|
||||
#
|
||||
# 复对 (p, p*)，取 imag(p)>0 的 p 作为首：
|
||||
# (原) φ_p = sqrt(2 Re(a_p)) (s - |a_p|)/[(s + a_p)(s + a_p^*)] * Π(...)
|
||||
# (原) φ_{p+1}= sqrt(2 Re(a_p)) (s + |a_p|)/[(s + a_p)(s + a_p^*)] * Π(...)
|
||||
# 代入 a_p=-p：
|
||||
# Re(a_p)= -Re(p)=σ>0, (s + a_p) = (s - p), (s + a_p^*)=(s - p^*)
|
||||
# |a_p| = |p|
|
||||
# 乘积同上 ⇒ Π_{i<p} (s + p_i^*)/(s - p_i)
|
||||
#
|
||||
# ⇒ φ_p(s) = sqrt(-2 Re(p)) (s - |p|)/[(s - p)(s - p^*)] * Π_{i<p} (s + p_i^*)/(s - p_i)
|
||||
# φ_{p^*}(s)= sqrt(-2 Re(p)) (s + |p|)/[(s - p)(s - p^*)] * Π_{i<p} (s + p_i^*)/(s - p_i)
|
||||
#
|
||||
# product 递推维护：prod_k = Π_{i≤k} (s + p_i^*)/(s - p_i)
|
||||
# 对复对要顺序乘两次（p 与 p*）。
|
||||
# ------------------------------------------------------------
|
||||
|
||||
class MuntzLaguerreIterator:
|
||||
def __init__(self, s: np.ndarray, stable_poles: list | np.ndarray):
|
||||
"""
|
||||
s: 复频率数组 (Nf,), s = j 2π f
|
||||
stable_poles: 稳定极点列表 (Re<0). 复共轭对要求正虚部在前 (p, p*).
|
||||
"""
|
||||
self.s = np.asarray(s, dtype=complex)
|
||||
self.poles = list(stable_poles)
|
||||
self.N = len(self.poles)
|
||||
self.k = 0
|
||||
# 初始化乘积 Π_{i<p} (s + p_i^*)/(s - p_i)
|
||||
self.product = np.ones_like(self.s, dtype=complex)
|
||||
|
||||
def __iter__(self):
|
||||
return self
|
||||
|
||||
def __next__(self):
|
||||
if self.k >= self.N:
|
||||
raise StopIteration
|
||||
|
||||
p = self.poles[self.k]
|
||||
if np.real(p) >= 0:
|
||||
raise ValueError(f"极点必须在左半平面: {p}")
|
||||
|
||||
# 复对首 (正虚部)
|
||||
if np.iscomplex(p) and np.imag(p) > 0:
|
||||
if self.k + 1 >= self.N:
|
||||
raise ValueError("复极点缺少共轭")
|
||||
pc = self.poles[self.k + 1]
|
||||
if not np.isclose(pc, np.conj(p)):
|
||||
raise ValueError("复极点未按 (p, p*) 顺序排列 (正虚部在前)")
|
||||
sigma = -np.real(p) # >0
|
||||
scale = np.sqrt(2 * sigma)
|
||||
r = np.abs(p)
|
||||
denom = (self.s - p) * (self.s - pc)
|
||||
|
||||
# 两个基函数
|
||||
phi_p = scale * (self.s - r) / denom * self.product
|
||||
phi_pc = scale * (self.s + r) / denom * self.product
|
||||
|
||||
# product 先乘 (s + p^*)/(s - p)，再乘 (s + p)/(s - p^*)
|
||||
self.product = self.product * (self.s + pc) / (self.s - p)
|
||||
self.product = self.product * (self.s + p) / (self.s - pc)
|
||||
|
||||
self.k += 2
|
||||
return [phi_p, phi_pc]
|
||||
|
||||
# 复对次 (负虚部) —— 应该被首元素处理，出现表示顺序错误
|
||||
if np.iscomplex(p) and np.imag(p) < 0:
|
||||
raise ValueError("检测到负虚部复极点但其共轭尚未处理，请将正虚部成员放在前面。")
|
||||
|
||||
# 实极点
|
||||
sigma = -np.real(p)
|
||||
if sigma <= 0:
|
||||
raise ValueError("实极点实部应为负 (稳定)。")
|
||||
scale = np.sqrt(2 * sigma)
|
||||
phi = scale / (self.s - p) * self.product
|
||||
# 更新乘积
|
||||
self.product = self.product * (self.s + p) / (self.s - p)
|
||||
|
||||
self.k += 1
|
||||
return [phi]
|
||||
|
||||
def generate_muntz_laguerre_basis(s: np.ndarray, stable_poles: list | np.ndarray):
|
||||
"""
|
||||
生成完整基函数列表: [φ_0=1, φ_1, φ_2, ...]
|
||||
"""
|
||||
basis = [np.ones_like(s, dtype=complex)]
|
||||
for block in MuntzLaguerreIterator(s, stable_poles):
|
||||
basis.extend(block)
|
||||
return basis
|
||||
|
||||
# ---------- 可选：离散再正交 (加权 QR) ----------
|
||||
# def trapezoid_weights(freqs: np.ndarray):
|
||||
# if len(freqs) == 1:
|
||||
# return np.ones(1)
|
||||
# df = np.diff(freqs)
|
||||
# w = np.zeros_like(freqs, dtype=float)
|
||||
# w[0] = 0.5 * df[0]
|
||||
# w[-1] = 0.5 * df[-1]
|
||||
# if len(freqs) > 2:
|
||||
# w[1:-1] = 0.5 * (df[:-1] + df[1:])
|
||||
# return w
|
||||
|
||||
def weighted_qr_from_basis(basis_cols: list[np.ndarray], weights: np.ndarray | None = None):
|
||||
A = np.column_stack(basis_cols) # (Nf, M)
|
||||
if weights is None:
|
||||
sw = np.ones(A.shape[0])
|
||||
else:
|
||||
sw = np.sqrt(weights.real)
|
||||
Aw = sw[:, None] * A
|
||||
Qw, R = np.linalg.qr(Aw)
|
||||
Phi = Qw / (sw[:, None] + 1e-30)
|
||||
return Phi, R # Raw = Phi R
|
||||
|
||||
def check_discrete_orthogonality(Phi: np.ndarray, w: np.ndarray):
|
||||
G = Phi.conj().T @ (w[:, None] * Phi)
|
||||
off = np.max(np.abs(G - np.eye(G.shape[0])))
|
||||
return G, off
|
||||
|
||||
def verify_orthonormal(Phi: np.ndarray, w: np.ndarray, atol=1e-10, rtol=1e-8):
|
||||
"""
|
||||
返回:
|
||||
G : Gram 矩阵 (Φ^H W Φ)
|
||||
max_off : 最大非对角幅值
|
||||
diag_err : max |diag(G)-1|
|
||||
passed : 是否满足阈值
|
||||
"""
|
||||
G = Phi.conj().T @ (w[:, None] * Phi)
|
||||
I = np.eye(G.shape[0])
|
||||
diag_err = np.max(np.abs(np.diag(G) - 1.0))
|
||||
max_off = np.max(np.abs(G - I + np.diag(np.diag(G) - 1.0)))
|
||||
passed = (diag_err <= atol) and (max_off <= rtol)
|
||||
return G, max_off, diag_err, passed
|
||||
|
||||
def omega_weights(freqs_hz: np.ndarray):
|
||||
"""
|
||||
基于 ω=2πf 的梯形法得到 w_ω = Δω/(2π)，使得
|
||||
(1/2π) ∫_{-∞}^{∞} → Σ w_ω,k
|
||||
"""
|
||||
f = freqs_hz
|
||||
if len(f) == 1:
|
||||
return np.ones(1)
|
||||
df = np.diff(f)
|
||||
w_f = np.zeros_like(f)
|
||||
w_f[0] = 0.5 * df[0]
|
||||
w_f[-1] = 0.5 * df[-1]
|
||||
if len(f) > 2:
|
||||
w_f[1:-1] = 0.5 * (df[:-1] + df[1:])
|
||||
# dω = 2π df, (1/2π) * dω = df ⇒ 直接 w_f 就是 w_ω
|
||||
return w_f # 已等价于 Δω/(2π)
|
||||
|
||||
# ------------------ 示例 ------------------
|
||||
if __name__ == "__main__":
|
||||
# 示例稳定极点 (复对正虚部在前)
|
||||
stable_poles = [
|
||||
-0.8e9,
|
||||
-1.0e9 + 2.5e9j,
|
||||
-1.0e9 - 2.5e9j,
|
||||
-2.2e9
|
||||
]
|
||||
freqs = np.linspace(1e8, 8e9, 400)
|
||||
s = 1j * 2 * np.pi * freqs
|
||||
|
||||
basis = generate_muntz_laguerre_basis(s, stable_poles)
|
||||
print("解析基函数数量 =", len(basis))
|
||||
print("前5个基函数示例 (每个前10个频点):")
|
||||
for i in range(5):
|
||||
print(f"φ_{i}:", basis[i][:10])
|
||||
|
||||
w = omega_weights(freqs)
|
||||
Phi_num, R = weighted_qr_from_basis(basis, w)
|
||||
Gram, off = check_discrete_orthogonality(Phi_num, w)
|
||||
print("离散 Gram 最大非对角元素 =", off)
|
||||
print("R 形状:", R.shape)
|
||||
# 验证 Raw ≈ Phi R
|
||||
raw = np.column_stack(basis)
|
||||
err = np.max(np.abs(raw - Phi_num @ R))
|
||||
print("重构误差 ||Raw - Phi R||_∞ =", err)
|
||||
|
||||
# 验证正交性
|
||||
print("离散 Gram 矩阵 (前5x5):")
|
||||
print(Gram[:5, :5])
|
||||
|
||||
Gcheck, max_off, diag_err, ok = verify_orthonormal(Phi_num, w)
|
||||
print(f"Diag 误差={diag_err:.3e}, Max off={max_off:.3e}, Orthonormal={ok}")
|
||||
# 额外: 验证 R
|
||||
# raw = Φ R => R ≈ Φ^H W raw (因为 Φ^H W Φ = I)
|
||||
R_alt = Phi_num.conj().T @ (w[:,None] * raw)
|
||||
print("R 差异 ||R - R_alt||_max =", np.max(np.abs(R - R_alt)))
|
||||
9
test/test_numpy.py
Normal file
9
test/test_numpy.py
Normal file
@@ -0,0 +1,9 @@
|
||||
import numpy as np
|
||||
|
||||
# 创建一个复数矩阵
|
||||
A = np.array([[1+2j, 2-1j], [3+4j, 4+0j]])
|
||||
|
||||
Q, R = np.linalg.qr(A)
|
||||
|
||||
print("Q =\n", Q)
|
||||
print("R =\n", R)
|
||||
102
test/test_orthonormal_basis.py
Normal file
102
test/test_orthonormal_basis.py
Normal file
@@ -0,0 +1,102 @@
|
||||
import numpy as np
|
||||
from core.orthonormal_basis import generate_muntz_laguerre_basis
|
||||
|
||||
# ---------- 可选：离散再正交 (加权 QR) ----------
|
||||
# def trapezoid_weights(freqs: np.ndarray):
|
||||
# if len(freqs) == 1:
|
||||
# return np.ones(1)
|
||||
# df = np.diff(freqs)
|
||||
# w = np.zeros_like(freqs, dtype=float)
|
||||
# w[0] = 0.5 * df[0]
|
||||
# w[-1] = 0.5 * df[-1]
|
||||
# if len(freqs) > 2:
|
||||
# w[1:-1] = 0.5 * (df[:-1] + df[1:])
|
||||
# return w
|
||||
|
||||
def weighted_qr_from_basis(basis_cols: list[np.ndarray], weights: np.ndarray | None = None):
|
||||
A = np.column_stack(basis_cols) # (Nf, M)
|
||||
if weights is None:
|
||||
sw = np.ones(A.shape[0])
|
||||
else:
|
||||
sw = np.sqrt(weights.real)
|
||||
Aw = sw[:, None] * A
|
||||
Qw, R = np.linalg.qr(Aw)
|
||||
Phi = Qw / (sw[:, None] + 1e-30)
|
||||
return Phi, R # Raw = Phi R
|
||||
|
||||
def check_discrete_orthogonality(Phi: np.ndarray, w: np.ndarray):
|
||||
G = Phi.conj().T @ (w[:, None] * Phi)
|
||||
off = np.max(np.abs(G - np.eye(G.shape[0])))
|
||||
return G, off
|
||||
|
||||
def verify_orthonormal(Phi: np.ndarray, w: np.ndarray, atol=1e-10, rtol=1e-8):
|
||||
"""
|
||||
返回:
|
||||
G : Gram 矩阵 (Φ^H W Φ)
|
||||
max_off : 最大非对角幅值
|
||||
diag_err : max |diag(G)-1|
|
||||
passed : 是否满足阈值
|
||||
"""
|
||||
G = Phi.conj().T @ (w[:, None] * Phi)
|
||||
I = np.eye(G.shape[0])
|
||||
diag_err = np.max(np.abs(np.diag(G) - 1.0))
|
||||
max_off = np.max(np.abs(G - I + np.diag(np.diag(G) - 1.0)))
|
||||
passed = (diag_err <= atol) and (max_off <= rtol)
|
||||
return G, max_off, diag_err, passed
|
||||
|
||||
def omega_weights(freqs_hz: np.ndarray):
|
||||
"""
|
||||
基于 ω=2πf 的梯形法得到 w_ω = Δω/(2π)，使得
|
||||
(1/2π) ∫_{-∞}^{∞} → Σ w_ω,k
|
||||
"""
|
||||
f = freqs_hz
|
||||
if len(f) == 1:
|
||||
return np.ones(1)
|
||||
df = np.diff(f)
|
||||
w_f = np.zeros_like(f)
|
||||
w_f[0] = 0.5 * df[0]
|
||||
w_f[-1] = 0.5 * df[-1]
|
||||
if len(f) > 2:
|
||||
w_f[1:-1] = 0.5 * (df[:-1] + df[1:])
|
||||
# dω = 2π df, (1/2π) * dω = df ⇒ 直接 w_f 就是 w_ω
|
||||
return w_f # 已等价于 Δω/(2π)
|
||||
|
||||
def evaluate(basis,freqs):
|
||||
w = omega_weights(freqs)
|
||||
Phi_num, R = weighted_qr_from_basis(basis, w)
|
||||
Gram, off = check_discrete_orthogonality(Phi_num, w)
|
||||
print("离散 Gram 最大非对角元素 =", off)
|
||||
print("R 形状:", R.shape)
|
||||
# 验证 Raw ≈ Phi R
|
||||
raw = np.column_stack(basis)
|
||||
err = np.max(np.abs(raw - Phi_num @ R))
|
||||
print("重构误差 ||Raw - Phi R||_∞ =", err)
|
||||
|
||||
# 验证正交性
|
||||
print("离散 Gram 矩阵 (前5x5):")
|
||||
print(Gram[:5, :5])
|
||||
|
||||
Gcheck, max_off, diag_err, ok = verify_orthonormal(Phi_num, w)
|
||||
print(f"Diag 误差={diag_err:.3e}, Max off={max_off:.3e}, Orthonormal={ok}")
|
||||
# 额外: 验证 R
|
||||
# raw = Φ R => R ≈ Φ^H W raw (因为 Φ^H W Φ = I)
|
||||
R_alt = Phi_num.conj().T @ (w[:,None] * raw)
|
||||
print("R 差异 ||R - R_alt||_max =", np.max(np.abs(R - R_alt)))
|
||||
|
||||
# ------------------ 示例 ------------------
|
||||
if __name__ == "__main__":
|
||||
# 示例稳定极点 (复对正虚部在前)
|
||||
init_poles = [
|
||||
-1.0e3 + 2.5e9j,
|
||||
-1.0e3 - 2.5e9j,
|
||||
]
|
||||
freqs = np.linspace(1e8, 8e9, 40)
|
||||
s = 1j * 2 * np.pi * freqs
|
||||
|
||||
basis = generate_muntz_laguerre_basis(s, init_poles)
|
||||
print("解析基函数数量 =", len(basis))
|
||||
print("基函数:")
|
||||
for i in range(len(basis)):
|
||||
print(f"φ_{i}:", basis[i][:])
|
||||
|
||||
evaluate(basis, freqs)
|
||||
Reference in New Issue
Block a user