feat: 整合了basicQR和relaxedQR

core/relaxed_basisQR.py

@@ -4,9 +4,10 @@ from scipy.linalg import block_diag
 import skrf as rf
 from skrf import VectorFitting
 from core.freqency import auto_select
+import random as rnd
 
 class RelaxedBasicBasisQR:
-    def __init__(self,H,freqs,poles,weights=None,passivity=True,dc_enforce=True,fit_constant=True):
+    def __init__(self,H,freqs,poles,weights=None,passivity=True,dc_enforce=True,fit_constant=True,fit_proportional=False):
         self.least_squares_rms_error  = None
         self.least_squares_condition = None
         self.eigenval_condition = None
@@ -16,6 +17,7 @@ class RelaxedBasicBasisQR:
 
         self.dc_enforce = dc_enforce
         self.fit_constant = fit_constant
+        self.fit_proportional = fit_proportional
         
         # self.H = H
         # self.freqs = freqs
@@ -27,7 +29,7 @@ class RelaxedBasicBasisQR:
         self.Phi = self.generate_basis(self.s, self.poles)
         self.A = self.matrix_A(self.poles)
         self.B = self.vector_B(self.poles)
-        self.Cw,self.w0 = self.fit_denominator(self.H, weights=weights)
+        self.Cw,self.w0,self.e = self.fit_denominator(self.H, weights=weights)
         self.D = self.w0
 
         self.Cr = None
@@ -119,7 +121,7 @@ class RelaxedBasicBasisQR:
             B = Bb if B is None else np.vstack([B, Bb])
         return B
     
-    def fit_denominator(self, H, weights=None):
+    def fit_denominator(self, H, weights=None, d0 = 1.0):
         """
         Solve formula (70) on the real basis Φ to obtain:
           - d (real)  → packs into C for this state's block structure
@@ -129,6 +131,10 @@ class RelaxedBasicBasisQR:
         H = np.asarray(H, np.complex128).reshape(-1,1)
         K, N = self.Phi.shape
         one  = np.ones((K, 1), np.complex128)
+        Phi   = self.Phi
+        dc_tol = 1e-18
+        has_dc = self.dc_enforce and self.freqs[0] < dc_tol
+        keep = np.ones(K, dtype=bool)
 
         # SK weighting (applied only to the (73) rows we keep in LS)
         if weights is None:
@@ -136,17 +142,16 @@ class RelaxedBasicBasisQR:
         else:
             weights = np.diag([1/res for res in weights])
 
-        Phi   = self.Phi
+        if self.fit_constant:  
-        Phi_w = np.hstack([one, Phi])   
+            Phi_w = np.hstack([one, Phi]) 
-        M     = np.hstack([Phi, -(H * Phi_w)])    # (K, 2N+1), complex
+            M     = np.hstack([Phi, -(H * Phi_w)])    # (K, 2N+1), complex
-
+        else:
-        dc_tol = 1e-18
+            M     = np.hstack([Phi, -(H * Phi)])       # (K, 2N),   complex
-        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)))
-            keep = np.ones(K, dtype=bool); keep[k0] = False
+            keep[k0] = False
             M_w     = weights @ M
             A_re  = np.real(M_w[keep, :])
             A_im  = np.imag(M_w[keep, :])
@@ -154,49 +159,78 @@ class RelaxedBasicBasisQR:
             # 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:
             M_w    = weights @ M
             A_re = np.real(M_w)
             A_im = np.imag(M_w)
             A_dc_re = A_dc_im = None
 
-        beta     = float(np.sqrt(np.sum(np.abs(H)**2)))
+        A_blocks = [A_re, A_im]
-        mean_row = (beta / K) * np.sum(Phi_w, axis=0)
+
-        A_w0      = np.concatenate([np.zeros(N, float),
+        if self.fit_constant:
-                               np.real(mean_row).astype(float)]
+            beta     = float(np.sqrt(np.sum(np.abs(H)**2)))
-                              ).reshape(1, -1)
+            mean_row = (beta / K) * np.sum(Phi_w, axis=0)
-        b_w0      = np.array([beta], float)
+            A_w0      = np.concatenate([np.zeros(N, float),
+                                np.real(mean_row).astype(float)]
+                                ).reshape(1, -1)
+            b_w0      = np.array([beta], float)
+
+            A_blocks += [A_w0]
+            m = A_re.shape[0] + A_im.shape[0]
+            b = np.zeros(m, float)
+            b = np.concatenate([b, b_w0])
+        else:
+            H_kp = (weights @ H)[keep,:]
+            b_re = np.real(d0 * H_kp)
+            b_im = np.imag(d0 * H_kp)
+            b = np.concatenate([b_re.ravel(), b_im.ravel()]).astype(float)
 
         # ---- build final stacked-real system ----
-        A_blocks = [A_re, A_im]
+        
 
         # if A_dc_re is not None:
         #     A_blocks += [A_dc_re, A_dc_im]
-
-        A_blocks += [A_w0]
-        A = np.vstack(A_blocks).astype(float)
-
-        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)] ----
+        A = np.vstack(A_blocks).astype(float)
         Q, R = np.linalg.qr(A, mode="reduced")
-        x = np.linalg.solve(R, Q.T @ b)
+        
+        if self.fit_constant:
+            Q2 = Q[:,A.shape[1]//2:]
+            R22 = R[A.shape[1]//2:,A.shape[1]//2:]
+        else:
+            Q2 = Q[:,A.shape[1]//2:]
+            R22 = R[A.shape[1]//2:,A.shape[1]//2:]
+
+        x = np.linalg.solve(R22, Q2.T @ b)
 
         # diagnostics
-        resid = A @ x - b
+        resid = Q2 @ R22 @ 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]
+        # cw = x[N:]                     # last (N+1) entries = [w0, w_1..w_N]
-        w0 = float(cw[0])
+        # w0 = float(cw[0])
-        Cw = cw[1:].reshape(1, N)      # row vector (1, N)
+        # Cw = cw[1:].reshape(1, N)      # row vector (1, N)
-        return Cw, w0
+        return self.extract_Cw_d_e(x,N,d0)
+    
+    def extract_Cw_d_e(self,C,N,d0=1.0):
+        if self.fit_proportional and self.fit_constant:
+            d = C[1]
+            e = C[0]
+            return C[2:].reshape(1, -1), d, e
+        elif self.fit_proportional and not self.fit_constant:
+            d = 0.0
+            e = C[0]
+            return C[1:].reshape(1, -1), d, e
+        elif not self.fit_proportional and self.fit_constant:
+            d = C[0]
+            e = 0.0
+            return C[1:].reshape(1, -1), d, e
+        else:
+            return C.reshape(1, -1), d0, 0.0
     
     
     def non_bias_Cr(self,w0):
@@ -215,12 +249,22 @@ class RelaxedBasicBasisQR:
         num = phi @ self.Cr
         H = num / den
         return H.ravel()
+    
+def noise(n:complex,coeff:float=0.05):
+    noise_r = rnd.gauss(-coeff * n.real, coeff * n.real)
+    noise_i = rnd.gauss(-coeff * n.imag, coeff * n.imag)
+    return complex(n.real + noise_r, n.imag + noise_i)
 
 if __name__ == "__main__":
     start_point = 0
     network = rf.Network("/tmp/paramer/simulation/3000/3000.s2p")
-    K = 5
+    K = 10
-    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)
+
+    full_freqences = network.f[start_point:]
+    noised_sampled_points = [(network.y[i][0][0]) for i in range(start_point,len(network.y))]
+    sampled_points = [network.y[i][0][0] for i in range(start_point,len(network.y))]
+
+    H11,freqs = auto_select(noised_sampled_points,full_freqences,max_points=20)
     poles = generate_starting_poles(2,beta_min=1e4,beta_max=freqs[-1]*1.1)
     
     Dt_1 = np.ones((len(freqs),1),np.complex128)
@@ -266,10 +310,9 @@ if __name__ == "__main__":
     fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False)
 
     ax00 = axes[0][0]
-    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
+    sliced_freqences = freqs
+
     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')