diff --git a/core/MultiplePortQR.py b/core/MultiplePortQR.py index 0f67eaa..58e1dff 100644 --- a/core/MultiplePortQR.py +++ b/core/MultiplePortQR.py @@ -3,7 +3,8 @@ from core.sk_iter import generate_starting_poles from scipy.linalg import block_diag import skrf as rf from skrf import VectorFitting -from core.freqency import auto_select +from core.freqency import auto_select_multple_ports +import matplotlib.pyplot as plt import random as rnd class MultiplePortQR: @@ -23,6 +24,7 @@ class MultiplePortQR: # self.freqs = freqs self.freqs = freqs self.H = H + self.ports = H.shape[1] self.s = self.freqs * 2j * np.pi self.P = len(poles) self.poles = poles @@ -128,7 +130,6 @@ class MultiplePortQR: - gamma (complex) Optional 'weights' (K,) apply row scaling: SK weighting if 1/|D_prev|. """ - H = np.asarray(H, np.complex128).reshape(-1,1) K, N = self.Phi.shape one = np.ones((K, 1), np.complex128) Phi = self.Phi @@ -137,40 +138,74 @@ class MultiplePortQR: keep = np.ones(K, dtype=bool) # 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]) - - if self.fit_constant: - Phi_w = np.hstack([one, Phi]) - M = np.hstack([Phi, -(H * Phi_w)]) # (K, 2N+1), complex - else: - M = np.hstack([Phi, -(H * Phi)]) # (K, 2N), complex if has_dc: # Enforce DC response exactly: k0 = int(np.argmin(np.abs(self.freqs))) keep[k0] = False - M_w = weights @ M - A_re = np.real(M_w[keep, :]) - A_im = np.imag(M_w[keep, :]) + + if self.fit_constant: + Phi_w = np.hstack([one, Phi]) + index = 0 + M_kp = None + for i in range(self.ports): + for j in range(self.ports): + M0 = np.zeros((K,N*self.ports**2),dtype=complex) + M0[:,index*N:(index+1)*N] = Phi + M0 = np.hstack([M0, -(H[:,i,j].reshape(-1,1) * Phi_w)]).reshape((K, -1))[keep,:] # (K, 2N), complex + index+=1 + M_kp = M0 if M_kp is None else np.vstack([M_kp, M0]) + assert M_kp is not None + else: + index = 0 + M_kp = None + for i in range(self.ports): + for j in range(self.ports): + M0 = np.zeros((K,N*ports**2),dtype=complex) + M0[:,index*N:(index+1)*N] = Phi + M0 = np.hstack([M0, -(H[:,i,j].reshape(-1,1) * Phi)]).reshape((K, -1))[keep,:] # (K, 2N), complex + index+=1 + M_kp = M0 if M_kp is None else np.vstack([M_kp, M0]) + assert M_kp is not None + + if weights is None: + weights_kp = np.diag(np.ones(len(freqs[keep]) * self.ports**2, np.complex128)) + else: + weights_kp0 = weights[keep] + weights0 = [] + for i in range(self.ports **2 ): + for res in weights_kp0: + weights0.append(1/res) + weights_kp = np.diag(np.array(weights0)) + + + if has_dc: + M_w_kp = weights_kp @ M_kp + A_re = np.real(M_w_kp) + A_im = np.imag(M_w_kp) mask = np.ones(K, dtype=bool); mask[k0] = False # exact (unweighted) DC rows: - A_dc_re = np.real(M[k0, :]).reshape(1, -1) - A_dc_im = np.imag(M[k0, :]).reshape(1, -1) + # A_dc_re = np.real(M_kp).reshape(1, -1) + # A_dc_im = np.imag(M_kp).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 + M_w_kp = weights_kp @ M_kp + A_re = np.real(M_w_kp) + A_im = np.imag(M_w_kp) + # A_dc_re = A_dc_im = None A_blocks = [A_re, A_im] if self.fit_constant: - 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), + Hk_kp = None + for i in range(self.ports): + for j in range(self.ports): + Hk_kp0 = H[:,i,j][keep] + Hk_kp = Hk_kp0 if Hk_kp is None else np.hstack([Hk_kp, Hk_kp0]) + assert Hk_kp is not None + Hk_sum = np.sum(np.abs(Hk_kp)**2) + beta = float(np.sqrt(Hk_sum)) + mean_row = (beta / weights_kp.shape[0]) * np.sum(Phi_w, axis=0) + A_w0 = np.concatenate([np.zeros(N*self.ports**2, float), np.real(mean_row).astype(float)] ).reshape(1, -1) b_w0 = np.array([beta], float) @@ -180,7 +215,14 @@ class MultiplePortQR: b = np.zeros(m, float) b = np.concatenate([b, b_w0]) else: - H_kp = (weights @ H)[keep,:] + H_kp = None + for i in range(self.ports): + for j in range(self.ports): + H_kp0 = weights_kp @ (H[:,i,j]).reshape(1,-1)[keep,:] + H_kp = H_kp0 if H_kp is None else np.hstack([H_kp, H_kp0]) + assert H_kp is not None + H_kp = H_kp.reshape(-1,1) + 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) @@ -197,11 +239,11 @@ class MultiplePortQR: Q, R = np.linalg.qr(A, mode="reduced") if self.fit_constant: - Q2 = Q[:,A.shape[1]//2:] - R22 = R[A.shape[1]//2:,A.shape[1]//2:] + Q2 = Q[:,Phi.shape[1] * self.ports**2:] + R22 = R[Phi.shape[1] * self.ports**2:,Phi.shape[1] * self.ports**2:] else: - Q2 = Q[:,A.shape[1]//2:] - R22 = R[A.shape[1]//2:,A.shape[1]//2:] + Q2 = Q[:,Phi.shape[1] * self.ports**2:] + R22 = R[Phi.shape[1] * self.ports**2:,Phi.shape[1] * self.ports**2:] x = np.linalg.solve(R22, Q2.T @ b) @@ -236,19 +278,27 @@ class MultiplePortQR: def non_bias_Cr(self,w0): A = np.asarray(self.Phi) den = np.diag((w0 + 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) + Cr = [] + for i in range(self.ports): + Cr.append([]) + for j in range(self.ports): + b = np.asarray(den) @ self.H[:,i,j].reshape(-1,1) + Cr_ij, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None) + Cr[i].append(Cr_ij) return Cr def evaluate(self,freqs, w0): + H = np.zeros((len(freqs),self.ports,self.ports),dtype=complex) s = 1j * 2*np.pi * np.asarray(freqs, float).ravel() phi = self.generate_basis(s, self.poles) den = w0 + phi @ self.Cw.T if self.Cr is None: self.Cr = self.non_bias_Cr(w0=w0) - num = phi @ self.Cr - H = num / den - return H.ravel() + for i in range(self.ports): + for j in range(self.ports): + num = phi @ self.Cr[i][j] + H[:,i,j] = (num / den).reshape(1,-1) + return H def noise(n:complex,coeff:float=0.05): noise_r = rnd.gauss(-coeff * n.real, coeff * n.real) @@ -257,19 +307,23 @@ def noise(n:complex,coeff:float=0.05): if __name__ == "__main__": start_point = 0 - network = rf.Network("/tmp/paramer/simulation/3000/3000.s2p") + network = rf.Network("/tmp/paramer/simulation/3500/3500.s2p") + ports = network.nports K = 10 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))] + noised_sampled_points = network.y[start_point:,:,:] + sampled_points = network.y[start_point:,:,:] - H11,freqs = auto_select(noised_sampled_points,full_freqences,max_points=20) + # noised_sampled_points = - network.y[start_point:,0,1].reshape(-1,1,1) + # sampled_points = network.y[start_point:,1,1].reshape(-1,1,1) + + H,freqs = auto_select_multple_ports(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) # Levi step (no weighting): - basis = MultiplePortQR(H11,freqs,poles=poles) + basis = MultiplePortQR(H,freqs,poles=poles) Dt = basis.Dt poles = basis.next_poles @@ -287,7 +341,7 @@ if __name__ == "__main__": eigenval_condition = [] eigenval_rms_error = [] for i in range(K): - basis = MultiplePortQR(H11,freqs,poles=poles,weights=Dt) + basis = MultiplePortQR(H,freqs,poles=poles,weights=Dt) Dt_1 = Dt Dt = basis.Dt poles = basis.next_poles @@ -305,55 +359,72 @@ 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[start_point:], w0=basis.w0) - import matplotlib.pyplot as plt - fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) - - ax00 = axes[0][0] - fitted_points = H11_evaluated + H_evaluated = basis.evaluate(full_freqences, w0=basis.w0) + fitted_points = H_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') - 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") + input_points = H + for i in range(ports): + for j in range(ports): + fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) + ax00 = axes[0][0] + ax00.plot(full_freqences, np.abs(sampled_points[:,i,j]), 'o', ms=4, color='red', label='Samples') + ax00.plot(full_freqences, np.abs(fitted_points[:,i,j]), '-', lw=2, color='k', label='Fit') + ax00.plot(sliced_freqences, np.abs(input_points[:,i,j]), 'x', ms=4, color='blue', label='Input Samples') + ax00.set_title(f"Response i={i+1}, j={j+1}") + ax00.set_ylabel("Magnitude") + ax00.legend(loc="best") - ax01 = axes[0][1] - ax01.set_title("Response i=0, j=0") - ax01.set_ylabel("Phase (deg)") - ax01.plot(network.f[start_point:], np.angle([network.y[i][0][0] for i in range(start_point,len(network.y))],deg=True), 'o', ms=4, color='red', label='Samples') - ax01.plot(network.f[start_point:], np.angle(H11_evaluated,deg=True), '-', lw=2, color='k', label='Fit') - ax01.plot(freqs, np.angle(H11,deg=True), 'x', ms=4, color='blue', label='Input Samples') - ax01.legend(loc="best") + ax01 = axes[0][1] + ax01.set_title(f"Response i={i+1}, j={j+1}") + ax01.set_ylabel("Phase (deg)") + ax01.plot(full_freqences, np.angle(sampled_points[:,i,j],deg=True), 'o', ms=4, color='red', label='Samples') + ax01.plot(full_freqences, np.angle(fitted_points[:,i,j],deg=True), '-', lw=2, color='k', label='Fit') + ax01.plot(sliced_freqences, np.angle(input_points[:,i,j],deg=True), 'x', ms=4, color='blue', label='Input Samples') + ax01.legend(loc="best") - ax10 = axes[1][0] - ax10.plot(least_squares_condition, label='Least Squares Condition') - ax10.set_title("least_squares_condition") - ax10.set_ylabel("Magnitude") - ax10.legend(loc="best") + # ax00 = axes[0][0] + # ax00.plot(full_freqences, np.real(sampled_points[:,i,j]), 'o', ms=4, color='red', label='Samples') + # ax00.plot(full_freqences, np.real(fitted_points[:,i,j]), '-', lw=2, color='k', label='Fit') + # ax00.plot(sliced_freqences, np.real(input_points[:,i,j]), 'x', ms=4, color='blue', label='Input Samples') + # ax00.set_title(f"Response i={i+1}, j={j+1}") + # ax00.set_ylabel("Real Part") + # ax00.legend(loc="best") - ax11 = axes[1][1] - ax11.plot(least_squares_rms_error, label='Least Squares RMS Error') - ax11.set_title("least_squares_rms_error") - ax11.set_ylabel("Magnitude") - ax11.legend(loc="best") + # ax01 = axes[0][1] + # ax01.set_title(f"Response i={i+1}, j={j+1}") + # ax01.set_ylabel("Imag Part") + # ax01.plot(full_freqences, np.imag(sampled_points[:,i,j]), 'o', ms=4, color='red', label='Samples') + # ax01.plot(full_freqences, np.imag(fitted_points[:,i,j]), '-', lw=2, color='k', label='Fit') + # ax01.plot(sliced_freqences, np.imag(input_points[:,i,j]), 'x', ms=4, color='blue', label='Input Samples') + # ax01.legend(loc="best") - ax20 = axes[2][0] - ax20.plot(eigenval_condition, label='Eigenvalue Condition') - ax20.set_title("eigenval_condition") - ax20.set_ylabel("Magnitude") - ax20.legend(loc="best") + ax10 = axes[1][0] + ax10.plot(least_squares_condition, label='Least Squares Condition') + ax10.set_title("least_squares_condition") + ax10.set_ylabel("Magnitude") + ax10.legend(loc="best") - ax21 = axes[2][1] - ax21.plot(eigenval_rms_error, label='Eigenvalue RMS Error') - ax21.set_title("eigenval_rms_error") - ax21.set_ylabel("Magnitude") - ax21.legend(loc="best") - fig.tight_layout() - plt.savefig(f"relaxed_basic_basis_QR.png") + ax11 = axes[1][1] + ax11.plot(least_squares_rms_error, label='Least Squares RMS Error') + ax11.set_title("least_squares_rms_error") + ax11.set_ylabel("Magnitude") + ax11.legend(loc="best") + + ax20 = axes[2][0] + ax20.plot(eigenval_condition, label='Eigenvalue Condition') + ax20.set_title("eigenval_condition") + ax20.set_ylabel("Magnitude") + ax20.legend(loc="best") + + ax21 = axes[2][1] + ax21.plot(eigenval_rms_error, label='Eigenvalue RMS Error') + ax21.set_title("eigenval_rms_error") + ax21.set_ylabel("Magnitude") + ax21.legend(loc="best") + fig.tight_layout() + plt.savefig(f"MultiplePortQR_port_{i+1}{j+1}.png") + print(f"Saved MultiplePortQR_port_{i+1}{j+1}.png") diff --git a/core/freqency.py b/core/freqency.py index 617d296..14f23a8 100644 --- a/core/freqency.py +++ b/core/freqency.py @@ -1,4 +1,135 @@ import numpy as np +# 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] + 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 @@ -6,49 +137,13 @@ def auto_select(H, freq, 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: + if N < 4 or max_points is None or max_points >= N: + # 直接返回所有点 return H.copy(), f.copy() eps = 1e-16 @@ -56,8 +151,6 @@ def auto_select(H, freq, 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: @@ -71,21 +164,17 @@ def auto_select(H, freq, 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: @@ -93,17 +182,14 @@ def auto_select(H, freq, 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: + if sel.size > max_points: priority = np.zeros(sel.size, dtype=int) if ensure_ends: priority[(sel == 0) | (sel == N-1)] = 3 @@ -112,7 +198,6 @@ def auto_select(H, freq, keep = [] budget = max_points - # keep highest-priority first for lev in (3, 2, 1, 0): cand = sel[priority == lev] if cand.size == 0: @@ -128,4 +213,36 @@ def auto_select(H, freq, break sel = np.array(sorted(set(keep)), dtype=int) - return H[sel], f[sel] \ No newline at end of file + if sel.size < max_points: + all_idx = set(range(N)) + missing = list(sorted(all_idx - set(sel))) + n_missing = max_points - sel.size + if n_missing > 0 and missing: + extra = np.linspace(0, len(missing)-1, n_missing, dtype=int) + sel = np.concatenate([sel, np.array(missing)[extra]]) + sel = np.array(sorted(set(sel)), dtype=int) + if sel.size < max_points: + left = list(sorted(all_idx - set(sel))) + if left: + sel = np.concatenate([sel, np.random.choice(left, max_points-sel.size, replace=False)]) + sel = np.array(sorted(set(sel)), dtype=int) + sel = sel[:max_points] + + return H[sel], f[sel] + +def auto_select_multple_ports(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): + ports = H.shape[1] + H_selected = np.zeros((max_points,ports,ports),dtype=complex) + for i in range(ports): + for j in range(ports): + H_selected[:,i,j], freq_selected = auto_select(H[:,i,j], freq, + n_baseline=n_baseline, peak_prominence=peak_prominence, + peak_window=peak_window, topgrad_q=topgrad_q, + max_points=max_points, ensure_ends=ensure_ends) + return H_selected, freq_selected diff --git a/core/relaxed_basisQR.py b/core/relaxed_basisQR.py index bb7b1b0..681655d 100644 --- a/core/relaxed_basisQR.py +++ b/core/relaxed_basisQR.py @@ -258,7 +258,7 @@ def noise(n:complex,coeff:float=0.05): if __name__ == "__main__": start_point = 0 network = rf.Network("/tmp/paramer/simulation/3000/3000.s2p") - K = 10 + K = 50 full_freqences = network.f[start_point:] noised_sampled_points = [(network.y[i][0][0]) for i in range(start_point,len(network.y))] @@ -354,6 +354,7 @@ if __name__ == "__main__": ax21.legend(loc="best") fig.tight_layout() plt.savefig(f"relaxed_basic_basis_QR.png") + print("Saved relaxed_basic_basis_QR.png")