447 lines
18 KiB
Python
447 lines
18 KiB
Python
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
|
||
from plotly.subplots import make_subplots
|
||
import plotly.graph_objs
|
||
from core.freqency import auto_select
|
||
from scipy.linalg import block_diag
|
||
|
||
|
||
network = rf.Network("/tmp/paramer/simulation/3500/3500.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")
|
||
|
||
def formula_67(s,y):
|
||
diag_values = [y[i][0][0] for i in range(len(y))]
|
||
H = np.diag(diag_values)
|
||
P = 2
|
||
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("Dt", D_t)
|
||
# print("D_t/D_t_pre",D_t/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.abs(N_t / D_t))
|
||
|
||
class formula_70_psi:
|
||
"""
|
||
VF-(70) with final pole-residue model.
|
||
After fit():
|
||
self.poles : (P,) complex
|
||
self.res : (P, M) complex # residues per response (columns)
|
||
self.h : (M,) complex # optional constants
|
||
self.g : (M,) complex # optional proportional terms
|
||
"""
|
||
|
||
# -------- internals --------
|
||
|
||
def __init__(self, s, P, H, beta_min, beta_max ,alpha_scale=0.01, n_iter=20, d0=1.0, include_const=True, include_linear=False, verbose=True):
|
||
self.s = s
|
||
self.P = P
|
||
self.H = H
|
||
self.beta_min = beta_min
|
||
self.beta_max = beta_max
|
||
self.alpha_scale = alpha_scale
|
||
self.z0 = self._generate_starting_poles(P, beta_min=beta_min, beta_max=beta_max,alpha_scale=alpha_scale)
|
||
self.z0 = np.array(self.z0, dtype=np.complex128)
|
||
self.n_iter = n_iter
|
||
self.d0 = d0
|
||
self.include_const = include_const
|
||
self.include_linear = include_linear
|
||
self.verbose = verbose
|
||
|
||
self.rel_iteration = []
|
||
self.cond_iteration = []
|
||
self.conf_poles_recognize = []
|
||
self.rms_error_iteration = []
|
||
self.rms_error_poles_recongnize = []
|
||
|
||
@staticmethod
|
||
def _generate_starting_poles(P, beta_min:float, beta_max:float, alpha_scale=0.01):
|
||
"""
|
||
仅生成复共轭对: 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] (正虚部先, 后跟共轭)
|
||
"""
|
||
betas = 2*np.pi*np.linspace(beta_min, beta_max, P)
|
||
poles = []
|
||
for b in betas:
|
||
alpha = alpha_scale * b
|
||
p = -alpha + 1j * b
|
||
poles += [p, np.conj(p)]
|
||
return poles
|
||
|
||
def _generate_laguerre_basis(self, s: np.ndarray, z: np.ndarray):
|
||
poles = z
|
||
poles = sorted(poles, key=lambda p: np.real(p))
|
||
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("复极点缺少共轭")
|
||
pn = poles[i]
|
||
pc = poles[i + 1]
|
||
if not np.isclose(pc, np.conj(pn)):
|
||
pc, pn = pn,pc
|
||
if not np.isclose(pc, np.conj(pn)):
|
||
raise ValueError("复极点未按 (p, p*) 顺序排列 (正虚部在前)")
|
||
poles[i], poles[i+1] = pc, pn # swap
|
||
sigma = -np.real(pn) # >0
|
||
scale = np.sqrt(2 * sigma)
|
||
r = np.abs(pn)
|
||
denom = (s - pn) * (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 - pn)
|
||
product = product * (s + pn) / (s - pc)
|
||
|
||
basis[i + 1] = phi_p
|
||
basis[i + 2] = phi_pc
|
||
i += 2
|
||
continue
|
||
|
||
# 复对次 (负虚部) —— 应该被首元素处理，出现表示顺序错误
|
||
if np.iscomplex(pn) and np.imag(pn) < 0:
|
||
raise ValueError("检测到负虚部复极点但其共轭尚未处理，请将正虚部成员放在前面。")
|
||
|
||
# 实极点
|
||
sigma = -np.real(pn)
|
||
if sigma <= 0:
|
||
raise ValueError("实极点实部应为负 (稳定)。")
|
||
scale = np.sqrt(2 * sigma)
|
||
phi = scale / (s - pn) * product
|
||
# 更新乘积
|
||
product = product * (s + pn) / (s - pn)
|
||
i += 1
|
||
basis[i + 1] = phi
|
||
return basis
|
||
|
||
|
||
def _orthonormal_psi(self,s,z):
|
||
s = np.asarray(s, np.complex128).reshape(-1)
|
||
z = np.asarray(z, np.complex128).reshape(-1)
|
||
return self._generate_laguerre_basis(s,z).T
|
||
|
||
@staticmethod
|
||
def _psi(s, z):
|
||
s = np.asarray(s, np.complex128).reshape(-1)
|
||
z = np.asarray(z, np.complex128).reshape(-1)
|
||
return 1.0 / (s[:, None] + z[None, :])
|
||
|
||
# @staticmethod
|
||
# def _lhp(z):
|
||
# z = np.asarray(z, np.complex128).reshape(-1).copy()
|
||
# z[z.real > 0] = -np.conj(z[z.real > 0])
|
||
# return z
|
||
|
||
def _build_70(self, s, H_list, z_ref, d0):
|
||
Hs = [np.asarray(h, np.complex128).reshape(-1) for h in H_list]
|
||
K, P, M = len(s), len(z_ref), len(Hs)
|
||
Psi = self._psi(s, z_ref)
|
||
Z = np.zeros((K, P))
|
||
rows, rhs = [], []
|
||
for m, H in enumerate(Hs):
|
||
Hp = H[:, None]
|
||
L_re = np.real(Hp * Psi); L_im = np.imag(Hp * Psi) # acts on c
|
||
R_re = -np.real(Psi); R_im = -np.imag(Psi) # acts on r^(m)
|
||
rows.append(np.hstack([L_re] + [R_re if j == m else Z for j in range(M)]))
|
||
rhs.append(-np.real(d0 * H))
|
||
rows.append(np.hstack([L_im] + [R_im if j == m else Z for j in range(M)]))
|
||
rhs.append(-np.imag(d0 * H))
|
||
A = np.vstack(rows)
|
||
b = np.concatenate(rhs)
|
||
cond_poles_recognize = np.linalg.cond(A)
|
||
rms_error_poles_recongnize = np.sqrt(np.mean(np.abs(A @ np.zeros(A.shape[1]) - b)**2))
|
||
return A, b, M, P, cond_poles_recognize, rms_error_poles_recongnize
|
||
|
||
def _step_70(self, s, H_list, z_ref, d0=1.0, scale=True):
|
||
A, b, M, P, cond_poles_recognize, rms_error_poles_recongnize = self._build_70(s, H_list, z_ref, d0)
|
||
if scale:
|
||
coln = np.maximum(np.linalg.norm(A, axis=0), 1e-12)
|
||
x, *_ = np.linalg.lstsq(A / coln, b, rcond=None)
|
||
x = x / coln
|
||
cond_iteration = np.linalg.cond(A)
|
||
rms_error_iteration = np.sqrt(np.mean(np.abs(A @ x - b)**2))
|
||
else:
|
||
x, *_ = np.linalg.lstsq(A, b, rcond=None)
|
||
cond_iteration = np.linalg.cond(A)
|
||
rms_error_iteration = np.sqrt(np.mean(np.abs(A @ x - b)**2))
|
||
|
||
c = x[:P]
|
||
res_ratio = np.empty((P, M), np.complex128)
|
||
off = P
|
||
for m in range(M):
|
||
res_ratio[:, m] = x[off:off+P]; off += P
|
||
# relocate poles with test matrix
|
||
Sigma = np.diag(-np.asarray(z_ref, np.complex128))
|
||
T = Sigma - (np.ones((P, 1), np.complex128) @ (c.reshape(1, -1) / d0))
|
||
z_new = -np.linalg.eigvals(T)
|
||
# z_new = self._lhp(z_new)
|
||
# cond = np.linalg.cond(T)
|
||
return z_new, c, cond_iteration, rms_error_iteration, cond_poles_recognize, rms_error_poles_recongnize, res_ratio
|
||
|
||
# -------- public API --------
|
||
def fit(self):
|
||
"""
|
||
s : (K,) complex samples (j*2π*f_sel)
|
||
H : (K,) complex or (K,M) or list of M vectors
|
||
z0: initial poles (P,) complex (LHP + conjugate pairs recommended)
|
||
"""
|
||
# normalize responses -> list
|
||
if isinstance(self.H, (list, tuple)):
|
||
H_list = [np.asarray(h, np.complex128).reshape(-1) for h in self.H]
|
||
else:
|
||
H_arr = np.asarray(self.H, np.complex128)
|
||
if H_arr.ndim == 1:
|
||
H_list = [H_arr]
|
||
elif H_arr.ndim == 2 and H_arr.shape[0] == len(self.s):
|
||
H_list = [H_arr[:, i].copy() for i in range(H_arr.shape[1])]
|
||
else:
|
||
raise ValueError("H must be (K,), list of (K,), or (K,M) with M responses.")
|
||
M = len(H_list)
|
||
|
||
# z = self._lhp(np.asarray(self.z0, np.complex128))
|
||
z = self.z0
|
||
# SK/VF relocations
|
||
for it in range(self.n_iter):
|
||
z_next, c_last, cond_iteration, rms_error_iteration, cond_poles_recognize, rms_error_poles_recongnize, _ = self._step_70(self.s, H_list, z, d0=self.d0, scale=True)
|
||
rel = np.linalg.norm(z_next) / max(1.0, np.linalg.norm(z))
|
||
if self.verbose:
|
||
print(f"[VF-70] iter {it+1:02d}/{self.n_iter:02d} Δz_rel={rel}")
|
||
self.cond_iteration.append(cond_iteration)
|
||
self.rms_error_iteration.append(rms_error_iteration)
|
||
self.conf_poles_recognize.append(cond_poles_recognize)
|
||
self.rms_error_poles_recongnize.append(rms_error_poles_recongnize)
|
||
self.rel_iteration.append(rel)
|
||
z = z_next
|
||
|
||
# ---- Finalize: refit residues with fixed poles z (pole–residue model) ----
|
||
Phi = self._psi(self.s, z) # K x P
|
||
# Build design matrix for extras
|
||
extras = []
|
||
if self.include_const: extras.append(np.ones(len(self.s), np.complex128))
|
||
if self.include_linear: extras.append(self.s.astype(np.complex128))
|
||
if extras:
|
||
X_base = np.column_stack([Phi] + extras) # K x (P + E)
|
||
else:
|
||
X_base = Phi
|
||
|
||
res = np.empty((len(z), M), np.complex128)
|
||
h = np.zeros(M, np.complex128)
|
||
g = np.zeros(M, np.complex128)
|
||
|
||
for m, Hm in enumerate(H_list):
|
||
theta, *_ = np.linalg.lstsq(X_base, Hm, rcond=None) # complex LS
|
||
res[:, m] = theta[:len(z)]
|
||
e = len(theta) - len(z)
|
||
if e >= 1: h[m] = theta[len(z)]
|
||
if e >= 2: g[m] = theta[len(z)+1]
|
||
|
||
# store the rational function
|
||
self.poles = z # (P,)
|
||
self.res = res # (P, M)
|
||
self.h = h # (M,)
|
||
self.g = g # (M,)
|
||
return self
|
||
|
||
# Evaluate the stored **rational** model on any grid
|
||
def evaluate(self, s_eval, m=None):
|
||
if not hasattr(self, "poles"):
|
||
raise RuntimeError("Model not fitted. Call fit(...) first.")
|
||
s_eval = np.asarray(s_eval, np.complex128).reshape(-1)
|
||
Phi = self._psi(s_eval, self.poles) # K_eval x P
|
||
if m is None:
|
||
H = Phi @ self.res
|
||
if np.any(self.h): H += self.h
|
||
if np.any(self.g): H += s_eval[:, None] * self.g
|
||
return H # (K_eval, M) or (K_eval,1)
|
||
m = int(m)
|
||
H = Phi @ self.res[:, m]
|
||
H += self.h[m]
|
||
H += s_eval * self.g[m]
|
||
return H # (K_eval,)
|
||
|
||
# def plot_rel_and_cond(self):
|
||
# fig, (ax1, ax2) = plt.subplots(2,1,figsize=(15,20))
|
||
# ax1.plot(np.log(self.rel), 'g-', label='rel')
|
||
# ax2.plot(np.log(self.cond), 'b-', label='cond')
|
||
# ax1.set_xlabel('Iteration')
|
||
# ax1.set_ylabel('Relative Change', color='g')
|
||
# ax2.set_ylabel('Condition Number', color='b')
|
||
# ax1.tick_params(axis='y', labelcolor='g')
|
||
# ax2.tick_params(axis='y', labelcolor='b')
|
||
# fig.tight_layout()
|
||
# plt.title('Relative Change and Condition Number per Iteration')
|
||
# # plt.show()
|
||
# plt.savefig(f"img_rel_cond.png")
|
||
def plot_rel_and_cond(self):
|
||
fig = make_subplots(rows=3, cols=2)
|
||
fig.add_trace(plotly.graph_objs.Scatter(y=self.rel_iteration, mode='lines+markers', name='rel'), row=1, col=1)
|
||
fig.add_trace(plotly.graph_objs.Scatter(y=self.cond_iteration, mode='lines+markers', name='cond_iteration'), row=2, col=1)
|
||
fig.add_trace(plotly.graph_objs.Scatter(y=self.rms_error_iteration, mode='lines+markers', name='rms_error_iteration'), row=3, col=1)
|
||
fig.add_trace(plotly.graph_objs.Scatter(y=self.conf_poles_recognize, mode='lines+markers', name='cond_poles_recognize'), row=2, col=2)
|
||
fig.add_trace(plotly.graph_objs.Scatter(y=self.rms_error_poles_recongnize, mode='lines+markers', name='rms_error_poles_recognize'), row=3, col=2)
|
||
fig.update_xaxes(title_text='Iteration', row=1, col=1)
|
||
fig.update_yaxes(title_text='Relative Change', row=1, col=1)
|
||
fig.update_yaxes(title_text='Condition Number (Iteration)', row=2, col=1)
|
||
fig.update_yaxes(title_text='RMS Error (Iteration)', row=3, col=1)
|
||
fig.update_yaxes(title_text='Condition Number (Pole Recognition)', row=2, col=2)
|
||
fig.update_yaxes(title_text='RMS Error (Pole Recognition)', row=3, col=2)
|
||
fig.update_layout(height=800, width=800, title_text='Relative Change and Condition Number per Iteration')
|
||
fig.write_image("img_rel_cond.png")
|
||
# fig.show()
|
||
|
||
def evaluate_on_freq(self, freq_eval, m=None):
|
||
return self.evaluate(1j * 2*np.pi * np.asarray(freq_eval, float), m=m)
|
||
|
||
|
||
|
||
if __name__ == "__main__":
|
||
# formula_67(s,y)
|
||
H11 = np.array([y[i,0,0] for i in range(len(y))])
|
||
H11_slice,freqs_slice = auto_select(H11,freqs,max_points=20)
|
||
s_slice = freqs_slice * 2j * np.pi
|
||
P_pairs = 2
|
||
|
||
K = 10
|
||
f70 = formula_70_psi(s_slice,P_pairs, H11_slice, beta_min=1e8, beta_max=freqs_slice[-1],alpha_scale=0.01, n_iter=K, d0=1.0, verbose=True)
|
||
model = f70.fit()
|
||
model.plot_rel_and_cond()
|
||
|
||
Hfit_dense = model.evaluate_on_freq(freqs)
|
||
|
||
fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=False)
|
||
|
||
# (1) Magnitude plot like your screenshot
|
||
ax0 = axes[0]
|
||
ax0.plot(freqs, np.abs(H11), 'o', ms=4, color='red', label='Samples')
|
||
ax0.plot(freqs, np.abs(Hfit_dense), '-', lw=2, color='k', label='Fit')
|
||
ax0.plot(freqs_slice, np.abs(H11_slice), 'x', ms=4, color='blue', label='Input Samples')
|
||
ax0.set_title("Response i=0, j=0")
|
||
ax0.set_ylabel("Magnitude")
|
||
ax0.legend(loc="best")
|
||
|
||
# (2) RMS error vs iteration
|
||
# ax1 = axes[1]
|
||
# its = np.arange(1, K+1)
|
||
# ax1.plot(its, hist["rms_rel"], '-o', lw=2)
|
||
# ax1.set_xlabel("Iteration")
|
||
# ax1.set_ylabel("RMS error (relative)")
|
||
# ax1.grid(True, alpha=0.3)
|
||
# ax1.set_title(f"RMS(final) = {hist['rms_rel'][-1]:.3e}")
|
||
|
||
fig.tight_layout()
|
||
plt.savefig(f"img_formula_70.png")
|
||
|
||
|