diff --git a/core/relaxed_basisQR.py b/core/relaxed_basisQR.py index b2591b8..75dd73c 100644 --- a/core/relaxed_basisQR.py +++ b/core/relaxed_basisQR.py @@ -6,17 +6,21 @@ from skrf import VectorFitting from core.freqency import auto_select class RelaxedBasicBasisQR: - def __init__(self,H,freqs,poles,weights=None,passivity=True,enforce_dc=True,fit_constant=True): + def __init__(self,H,freqs,poles,weights=None,passivity=True,dc_enforce=True,fit_constant=True): self.least_squares_rms_error = None self.least_squares_condition = None self.eigenval_condition = None self.eigenval_rms_error = None - self.dc_enforce = enforce_dc + self.dc_tol = 1e-18 + + self.dc_enforce = dc_enforce self.fit_constant = fit_constant - self.H = H + # self.H = H + # self.freqs = freqs self.freqs = freqs + self.H = H self.s = self.freqs * 2j * np.pi self.P = len(poles) self.poles = poles @@ -122,66 +126,78 @@ class RelaxedBasicBasisQR: - gamma (complex) Optional 'weights' (K,) apply row scaling: SK weighting if 1/|D_prev|. """ - if weights is None: - weights = np.diag(np.ones(len(H), np.complex128)) - else: - weights = np.diag([1/res for res in weights]) H = np.asarray(H, np.complex128).reshape(-1,1) K, N = self.Phi.shape one = np.ones((K, 1), np.complex128) - psi = weights @ self.Phi - psi_w = np.hstack([weights@one, psi]) - Phi_w = np.hstack([one, self.Phi]) + # SK weighting (applied only to the (73) rows we keep in LS) + if weights is None: + weights = np.diag(np.ones(len(H), np.complex128)) + else: + weights = np.diag([1/res for res in weights]) - beta = np.linalg.norm((weights@H)) * len(H) - mean_row = (beta / K) * np.sum(Phi_w, axis=0) - - Hpsi_w = H * psi_w - - + Phi = self.Phi + Phi_w = np.hstack([one, Phi]) + M = np.hstack([Phi, -(H * Phi_w)]) # (K, 2N+1), complex dc_tol = 1e-18 has_dc = self.dc_enforce and self.freqs[0] < dc_tol if has_dc: # Enforce DC response exactly: - k0 = int(np.argmin(np.abs(self.freqs))) + k0 = int(np.argmin(np.abs(self.freqs))) + keep = np.ones(K, dtype=bool); keep[k0] = False + M_w = weights @ M + A_re = np.real(M_w[keep, :]) + A_im = np.imag(M_w[keep, :]) mask = np.ones(K, dtype=bool); mask[k0] = False - # Weighted rows for k≠0 - A_re = np.hstack([np.real(psi)[mask:], np.real(-Hpsi_w)]) - A_im = np.hstack([np.imag(psi)[mask:], np.imag(-Hpsi_w)]) + # exact (unweighted) DC rows: + A_dc_re = np.real(M[k0, :]).reshape(1, -1) + A_dc_im = np.imag(M[k0, :]).reshape(1, -1) else: - A_re = np.hstack([np.real(psi), np.real(-Hpsi_w)]) - A_im = np.hstack([np.imag(psi), np.imag(-Hpsi_w)]) - A_w_re = np.concatenate([np.zeros(N, float), np.real(mean_row).astype(float)]).reshape(1, -1) + M_w = weights @ M + A_re = np.real(M_w) + A_im = np.imag(M_w) + A_dc_re = A_dc_im = None - b_re = np.zeros_like(H) - b_im = np.zeros_like(H) - b_w_re = beta + beta = float(np.sqrt(np.sum(np.abs(H)**2))) + mean_row = (beta / K) * np.sum(Phi_w, axis=0) + A_w0 = np.concatenate([np.zeros(N, float), + np.real(mean_row).astype(float)] + ).reshape(1, -1) + b_w0 = np.array([beta], float) - A = np.vstack([A_re, A_im, A_w_re]).astype(float) - Q,R = np.linalg.qr(A) - R22 = R[A_re.shape[1]//2:, A_re.shape[1]//2:] - Q2 = Q[:, A_re.shape[1]//2:] - # rown = np.linalg.norm(A, axis=1) - # rown = np.sqrt(rown) - # A = rown[:,None] * A - b = np.concatenate([b_re, b_im, [[b_w_re]]]).astype(float) + # ---- build final stacked-real system ---- + A_blocks = [A_re, A_im] - x = np.linalg.inv(R22) @ (Q2.T @ b) + # if A_dc_re is not None: + # A_blocks += [A_dc_re, A_dc_im] - self.least_squares_rms_error = np.sqrt(np.mean((Q2 @ R22 @ x - b)**2)) - self.least_squares_condition = np.linalg.cond(Q2 @ R22) + A_blocks += [A_w0] + A = np.vstack(A_blocks).astype(float) - Cw,w0 = self.vector_Cw(x,psi_w0=psi_w[:,0]) - return Cw,w0 + m = A_re.shape[0] + A_im.shape[0] + b = np.zeros(m, float) + # if A_dc_re is not None: + # b = np.concatenate([b, np.zeros(2, float)]) # DC rows → 0 + b = np.concatenate([b, b_w0]) + + # ---- QR solve for x = [c_H (N); c_w (N+1)] ---- + Q, R = np.linalg.qr(A, mode="reduced") + x = np.linalg.solve(R, Q.T @ b) + + # diagnostics + resid = A @ x - b + self.least_squares_rms_error = float(np.sqrt(np.mean(resid**2))) + self.least_squares_condition = float(np.linalg.cond(R)) + + # split cw and return + cw = x[N:] # last (N+1) entries = [w0, w_1..w_N] + w0 = float(cw[0]) + Cw = cw[1:].reshape(1, N) # row vector (1, N) + return Cw, w0 - def vector_Cw(self,x,psi_w0): - w0 = x[0] - C = x[1:] - return C.T,w0 def non_bias_Cr(self,w0): A = np.asarray(self.Phi) @@ -201,11 +217,11 @@ class RelaxedBasicBasisQR: return H.ravel() if __name__ == "__main__": - start_point = 2 + start_point = 0 network = rf.Network("/tmp/paramer/simulation/3000/3000.s2p") K = 5 H11,freqs = auto_select([network.y[i][0][0] for i in range(start_point,len(network.y))],network.f[start_point:],max_points=20) - poles = generate_starting_poles(2,beta_min=freqs[0]/1.1,beta_max=freqs[-1]*1.1) + poles = generate_starting_poles(2,beta_min=1e4,beta_max=freqs[-1]*1.1) Dt_1 = np.ones((len(freqs),1),np.complex128) # Levi step (no weighting): @@ -250,9 +266,14 @@ if __name__ == "__main__": fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) ax00 = axes[0][0] - ax00.plot(network.f[start_point:], np.abs([network.y[i][0][0] for i in range(start_point,len(network.y))]), 'o', ms=4, color='red', label='Samples') - ax00.plot(network.f[start_point:], np.abs(H11_evaluated), '-', lw=2, color='k', label='Fit') - ax00.plot(freqs, np.abs(H11), 'x', ms=4, color='blue', label='Input Samples') + full_freqences = network.f[start_point:] + sliced_freqences = freqs + sampled_points = [network.y[i][0][0] for i in range(start_point,len(network.y))] + fitted_points = H11_evaluated + input_points = H11 + ax00.plot(full_freqences, np.abs(sampled_points), 'o', ms=4, color='red', label='Samples') + ax00.plot(full_freqences, np.abs(fitted_points), '-', lw=2, color='k', label='Fit') + ax00.plot(sliced_freqences, np.abs(input_points), 'x', ms=4, color='blue', label='Input Samples') ax00.set_title("Response i=0, j=0") ax00.set_ylabel("Magnitude") ax00.legend(loc="best")