algorithm/ovf
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 1 Wiki Activity

feat: Add dc enforce

Browse Source
This commit is contained in:
mayge
2025-09-22 04:42:51 -04:00
parent 63d4de7076
commit 966adfe140
1 changed files with 69 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files

117
core/relaxed_basisQR.py
View File

@@ -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")