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") 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: # ---------- helpers ---------- def _psi(self,s, z_prev): """Psi[k,p] = 1/(s_k + z_prev[p]) (shape: Kf x P)""" s = np.asarray(s, dtype=np.complex128).reshape(-1) z_prev = np.asarray(z_prev, dtype=np.complex128).reshape(-1) return 1.0 / (s[:, None] + z_prev[None, :]) def _enforce_lhp_conjugate(self, poles): """Reflect any RHP poles across the imag axis and pair conjugates loosely.""" z = np.array(poles, dtype=np.complex128) for i in range(len(z)): if z[i].real > 0: z[i] = -np.conj(z[i]) return z # ---------- build Ax=b for (70) ----------1 def build_vf70_system(self,s, H_list, z_prev, d0=1.0): """ s: (Kf,) complex sample points (e.g. j*2πf) H_list: list of (Kf,) complex responses, one per MIMO element you want to fit z_prev: (P,) complex poles from previous iteration d0: fixed scalar, usually 1.0 Returns: A (2*Kf*M, (M+1)*P), b (2*Kf*M,), plus Psi for convenience Unknown vector x = [c (P), r^(1) (P), ..., r^(M) (P)] """ H_list = [np.asarray(h, dtype=np.complex128).reshape(-1) for h in H_list] Kf = len(s) P = len(z_prev) M = len(H_list) Psi = self._psi(s, z_prev) # Kf x P zerosP = np.zeros((Kf, P)) # convenience rows = [] rhs = [] for m, H in enumerate(H_list): Hp = H[:, None] # Kf x 1 # Left block (shared c): H * Psi L_re = np.real(Hp * Psi) # Kf x P L_im = np.imag(Hp * Psi) # Kf x P # Right block for this response: -Psi R_re = -np.real(Psi) # Kf x P R_im = -np.imag(Psi) # Kf x P # Assemble columns: [c | r^(1) | r^(2) | ... | r^(M)] # Real rows cols_re = [L_re] for j in range(M): cols_re.append(R_re if j == m else zerosP) rows.append(np.hstack(cols_re)) # Imag rows cols_im = [L_im] for j in range(M): cols_im.append(R_im if j == m else zerosP) rows.append(np.hstack(cols_im)) # RHS (move d0*H to the right) rhs.append(-np.real(d0 * H)) rhs.append(-np.imag(d0 * H)) A = np.vstack(rows) # (2*Kf*M) x ((M+1)*P) b = np.concatenate(rhs) # (2*Kf*M,) return A, b, Psi # ---------- one relocation step ---------- def vf70_step_once(self,s, H_list, z_prev, d0=1.0, column_scale=True): """ Returns: z_new (P,), c (P,), residues (list of M arrays length P) """ P = len(z_prev) M = len(H_list) A, b, Psi = self.build_vf70_system(s, H_list, z_prev, d0=d0) # optional simple column scaling for conditioning if column_scale: col_norm = np.maximum(np.linalg.norm(A, axis=0), 1e-12) x, *_ = np.linalg.lstsq(A / col_norm, b, rcond=None) x = x / col_norm else: x, *_ = np.linalg.lstsq(A, b, rcond=None) # unpack unknowns c = x[:P] residues = [] off = P for _ in range(M): residues.append(x[off:off+P]) off += P # pole relocation (zeros of d0 + Psi @ c) Sigma = np.diag(-np.asarray(z_prev, dtype=np.complex128)) T = Sigma - (np.ones((P, 1), dtype=np.complex128) @ (c.reshape(1, -1) / d0)) z_new = -np.linalg.eigvals(T) z_new = self._enforce_lhp_conjugate(z_new) return z_new, c, residues # ---------- iterate K times (K == n_iter) ---------- def vf70_iterate(self,s, H_list, z0, n_iter=15, d0=1.0, verbose=True): z = np.array(z0, dtype=np.complex128) c = None residues = None for it in range(n_iter): z_next, c, residues = self.vf70_step_once(s, H_list, z, d0=d0) if verbose: rel = np.linalg.norm(z_next) / max(1.0, np.linalg.norm(z)) print(f"[VF-70] iter {it+1:02d}/{n_iter:02d} Δz_rel = {rel}") z = z_next return z, c, residues def eval_Hfit_from_70(self,s, z_prev, c, residues_m, d0=1.0): """ Evaluate H_fit for one response m using (70): H_fit = (Psi @ residues_m) / (d0 + Psi @ c) """ Psi = 1.0 / (s[:, None] + np.asarray(z_prev)[None, :]) D_ratio = d0 + Psi @ c N_part = Psi @ residues_m return N_part / D_ratio def rms_abs(self,y): return np.sqrt(np.mean(np.abs(y)**2)) def rms_error(self,H, Hfit): """ Returns (abs_rms, relative_rms). Relative RMS is normalized by RMS of H. """ err = Hfit - H abs_rms = self.rms_abs(err) ref = max(self.rms_abs(H), 1e-16) rel_rms = abs_rms / ref return abs_rms, rel_rms # --- run K (== n_iter) VF-(70) steps and record history --- def run_vf70_with_history(self,s, H, z0, K=25, d0=1.0, verbose=True): """ H: 1D complex array (single response). If you want multi-response, pass a list [H11, H21, ...] and adapt the evaluation below. """ # Wrap single response as a list for the step function H_list = [np.asarray(H, dtype=np.complex128)] z = np.array(z0, dtype=np.complex128) hist = { "z": [z.copy()], "c": [], "r": [], "Hfit": [], "rms_abs": [], "rms_rel": [], } for it in range(K): z_new, c, residues = self.vf70_step_once(s, H_list, z, d0=d0, column_scale=True) # Evaluate fitted response after this iteration Hfit = self.eval_Hfit_from_70(s, z, c, residues[0], d0=d0) abs_r, rel_r = self.rms_error(H, Hfit) if verbose: drel = np.linalg.norm(z_new - z) / max(1.0, np.linalg.norm(z)) print(f"[VF-70] iter {it+1:02d}/{K:02d} Δz_rel={drel:.3e} RMS={abs_r:.3e} RMSrel={rel_r:.3e}") # save hist["c"].append(c.copy()) hist["r"].append(residues[0].copy()) hist["Hfit"].append(Hfit.copy()) hist["rms_abs"].append(abs_r) hist["rms_rel"].append(rel_r) z = z_new hist["z"].append(z.copy()) return hist def evaluate(self, s_eval, z_ref, c, residues, d0=1.0): """ Evaluate the (70) model on any frequency grid. Parameters ---------- s_eval : (K_eval,) complex Points where to evaluate (e.g., 1j*2π*freq_eval). z_ref : (P,) complex The *reference* poles used in the LS system that produced (c, residues). For the last iteration in run_vf70_with_history, this is hist['z'][-2]. c : (P,) complex Denominator-update coefficients from (70). residues : (P,) complex OR list of (P,) for multi-response Residues for the response(s) from (70). d0 : scalar, default 1.0 Returns ------- H_eval : (K_eval,) complex OR (K_eval, M) for M responses """ s_eval = np.asarray(s_eval, dtype=np.complex128).reshape(-1) z_ref = np.asarray(z_ref, dtype=np.complex128).reshape(-1) Psi_eval = 1.0 / (s_eval[:, None] + z_ref[None, :]) # K_eval x P Den = d0 + Psi_eval @ np.asarray(c, dtype=np.complex128) # Single response if not isinstance(residues, (list, tuple)): r = np.asarray(residues, dtype=np.complex128).reshape(-1) Num = Psi_eval @ r return Num / Den # Multi-response M = len(residues) H_eval = np.empty((s_eval.size, M), dtype=np.complex128) for m, r_m in enumerate(residues): r_m = np.asarray(r_m, dtype=np.complex128).reshape(-1) H_eval[:, m] = (Psi_eval @ r_m) / Den return H_eval def auto_select(H, freq, n_baseline=64, # log-spaced backbone points peak_prominence=0.05, # fraction of |H| dB dynamic range for peak detection peak_window=5, # take ±peak_window samples around each peak topgrad_q=0.98, # keep top 2% largest slope/phase-change points max_points=25, # final cap on selected samples (None = no cap) ensure_ends=True): """ Select several significant sample points for vector fitting. Strategy: 1) Always keep endpoints (optional). 2) Add a log-spaced baseline over the band. 3) Detect resonance peaks in |H| (on a log scale) and keep small windows around them. 4) Add points with the largest magnitude slope and phase-change (w.r.t log-f). 5) De-duplicate, sort, and optionally thin to 'max_points' with priority to endpoints and detected peaks. Parameters ---------- H : (N,) complex array Frequency response samples. freq : (N,) float array Frequency axis [Hz], strictly increasing. n_baseline : int Count of log-spaced baseline samples across the band. peak_prominence : float Peak prominence threshold as a fraction of the dynamic range in log|H|. 0.05 ≈ keep peaks ≥ 5% of the range. peak_window : int Number of neighbor indices to include on each side of every detected peak. topgrad_q : float in (0,1) Quantile for selecting strong slope/phase points. 0.98 ⇒ keep the top 2% largest derivatives. max_points : int or None If not None, cap the total number of selected indices to this value. ensure_ends : bool Always include the first and last samples. Returns ------- H_sel : (K,) complex array freq_sel : (K,) float array """ H = np.asarray(H).reshape(-1) f = np.asarray(freq).reshape(-1) if H.size != f.size: raise ValueError("H and freq must have the same length.") N = f.size if N < 4: return H.copy(), f.copy() eps = 1e-16 mag = np.abs(H) logmag = np.log10(mag + eps) phase = np.unwrap(np.angle(H)) # log-frequency axis (scale-invariant derivatives) # keep it linear if any non-positive freq sneaks in if np.all(f > 0): lf = np.log(f) else: lf = f.copy() dlf = np.gradient(lf) d_logmag = np.gradient(logmag) / (dlf + 1e-16) d_phase = np.gradient(phase) / (dlf + 1e-16) idx = set() if ensure_ends: idx.update([0, N-1]) # 1) log-spaced baseline if n_baseline > 0: # map a log grid to nearest indices grid = np.linspace(lf.min(), lf.max(), n_baseline) base_idx = np.clip(np.searchsorted(lf, grid), 0, N-1) idx.update(np.unique(base_idx).tolist()) # 2) peaks in |H| try: from scipy.signal import find_peaks dyn = logmag.max() - logmag.min() prom = peak_prominence * (dyn + 1e-12) peaks, _ = find_peaks(logmag, prominence=prom) except Exception: # simple fallback: strict local maxima peaks = np.where((mag[1:-1] > mag[:-2]) & (mag[1:-1] > mag[2:]))[0] + 1 for p in peaks: lo = max(0, p - peak_window) hi = min(N, p + peak_window + 1) idx.update(range(lo, hi)) # 3) strongest slope / phase-change points thr_slope = np.quantile(np.abs(d_logmag), topgrad_q) thr_phase = np.quantile(np.abs(d_phase), topgrad_q) idx.update(np.where(np.abs(d_logmag) >= thr_slope)[0].tolist()) idx.update(np.where(np.abs(d_phase) >= thr_phase)[0].tolist()) # 4) finalize set sel = np.array(sorted(idx), dtype=int) # 5) optional thinning with priority to endpoints and peaks if max_points is not None and sel.size > max_points: priority = np.zeros(sel.size, dtype=int) if ensure_ends: priority[(sel == 0) | (sel == N-1)] = 3 if peaks.size: priority[np.isin(sel, peaks)] = np.maximum(priority[np.isin(sel, peaks)], 2) keep = [] budget = max_points # keep highest-priority first for lev in (3, 2, 1, 0): cand = sel[priority == lev] if cand.size == 0: continue if cand.size <= budget: keep.extend(cand.tolist()) budget -= cand.size else: step = max(1, int(np.ceil(cand.size / budget))) keep.extend(cand[::step][:budget].tolist()) budget = 0 if budget == 0: break sel = np.array(sorted(set(keep)), dtype=int) return H[sel], f[sel] 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 z0 = generate_starting_poles(P_pairs, beta_min=1e8, beta_max=freqs_slice[-1]) z0 = np.array(z0, dtype=np.complex128) f70 = formula_70() K = 25 hist = f70.run_vf70_with_history(s_slice, H11_slice, z0, K=K, d0=1.0, verbose=True) it_show = K-1 Hfit_final = hist["Hfit"][it_show] s_dense = 1j * 2*np.pi * freqs c_final = hist["c"][it_show] r_final = hist["r"][it_show] # single response z_ref = hist["z"][it_show] Hfit_dense = f70.evaluate(s_dense, z_ref, c_final, r_final, d0=1.0) 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")