diff --git a/core/basisQR.py b/core/basisQR.py index 2f16635..c32d05c 100644 --- a/core/basisQR.py +++ b/core/basisQR.py @@ -5,7 +5,7 @@ import skrf as rf from skrf import VectorFitting from core.freqency import auto_select -class BasicBasis: +class BasicBasisQR: def __init__(self,H,freqs,poles,weights=None): self.least_squares_rms_error = None self.least_squares_condition = None @@ -21,7 +21,8 @@ class BasicBasis: self.A = self.matrix_A(self.poles) self.B = self.vector_B(self.poles) self.D = 1.0 - self.Cr,self.Cw = self.fit_denominator(self.H, d0=1.0, weights=weights) + self.Cw = self.fit_denominator(self.H, d0=1.0, weights=weights) + self.Cr = None z = np.linalg.eigvals(self.A - self.B @ self.Cw) p_next = -z @@ -128,29 +129,40 @@ class BasicBasis: b_im = np.imag(d0 * H) A = np.vstack([A_re, A_im]).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]).astype(float) - x = np.linalg.inv(A.T @ A) @ A.T @ b + x = np.linalg.inv(R22) @ (Q2.T @ b) - self.least_squares_rms_error = np.sqrt(np.mean((A @ x - b)**2)) - self.least_squares_condition = np.linalg.cond(A) + self.least_squares_rms_error = np.sqrt(np.mean((Q2 @ R22 @ x - b)**2)) + self.least_squares_condition = np.linalg.cond(Q2 @ R22) - Cn,Cd = self.vector_C(x) - return Cn,Cd + Cw = self.vector_Cw(x) + return Cw - def vector_C(self,x): - Cn = np.asarray([x[:len(x)//2]], float).reshape(1,-1) - Cd = np.asarray([x[len(x)//2:]], float).reshape(1,-1) - return Cn, Cd + def vector_Cw(self,x): + return x.T - def evaluate(self,freqs,poles,Cn,Cd,d0=1.0): + def non_bias_Cr(self,d0 = 1.0): + A = np.asarray(self.Phi) + den = np.diag((d0 + self.Phi @ self.Cw.T).ravel()) + b = np.asarray(den) @ self.H.reshape(-1,1) + Cr, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None) + return Cr + + def evaluate(self,freqs, d0=1.0): s = 1j * 2*np.pi * np.asarray(freqs, float).ravel() - phi = self.generate_basis(s, poles) - num = phi @ Cn.T - den = d0 + phi @ Cd.T + phi = self.generate_basis(s, self.poles) + den = d0 + phi @ self.Cw.T + if self.Cr is None: + self.Cr = self.non_bias_Cr(d0=d0) + num = phi @ self.Cr H = num / den return H.ravel() @@ -162,7 +174,7 @@ if __name__ == "__main__": Dt_1 = np.ones((len(freqs),1),np.complex128) # Levi step (no weighting): - basis = BasicBasis(H11,freqs,poles=poles) + basis = BasicBasisQR(H11,freqs,poles=poles) Dt = basis.Dt poles = basis.next_poles @@ -180,7 +192,7 @@ if __name__ == "__main__": eigenval_condition = [] eigenval_rms_error = [] for i in range(K): - basis = BasicBasis(H11,freqs,poles=poles,weights=Dt) + basis = BasicBasisQR(H11,freqs,poles=poles,weights=Dt) Dt_1 = Dt Dt = basis.Dt poles = basis.next_poles @@ -198,7 +210,7 @@ if __name__ == "__main__": eigenval_rms_error.append(basis.eigenval_rms_error) # H11_evaluated = basis.evaluate_pole_residue(network.f[1:],poles,basis.C[0]) - H11_evaluated = basis.evaluate(network.f[2:], poles, basis.Cr[0],basis.Cw[0], d0=1.0) + H11_evaluated = basis.evaluate(network.f[2:], d0=1.0) import matplotlib.pyplot as plt fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) @@ -234,7 +246,7 @@ if __name__ == "__main__": ax4.set_ylabel("Magnitude") ax4.legend(loc="best") fig.tight_layout() - plt.savefig(f"basic_basis.png") + plt.savefig(f"basic_basis_QR.png") diff --git a/img.png b/img.png deleted file mode 100644 index fd59237..0000000 Binary files a/img.png and /dev/null differ