diff --git a/app/cocal_methods.py b/app/cocal_methods.py index 8f956534cbb03a3e60219f9cb3bb3803fa5e263f..9f73b6b63be31e69c9d117b5fd20c6bda47c7287 100644 --- a/app/cocal_methods.py +++ b/app/cocal_methods.py @@ -805,6 +805,9 @@ class CocalMethods: mask = np.logical_and(self.ref_frequency < 15000.0, self.ref_frequency >= 0.0) mask_ri = np.r_[mask, mask] + # custom weights to reduce influence of higher frequencies + weights = np.abs(np.diff(np.log10(np.abs(self.ref_frequency)+1), prepend=1)) + # function to draw samples for Monte Carlo def draw_samples(size): return real_imag_2_complex( @@ -827,6 +830,7 @@ class CocalMethods: theta0=theta0, Nb=Nb, w_empirical_exp=w_empirical_exp, + weights=weights[mask], ) umc_kwargs = { @@ -876,7 +880,7 @@ class CocalMethods: # main calculation extracted from scipy.signal.freqz return npp.polyval(w_exp, b, tensor=False) / npp.polyval(w_exp, a, tensor=False) - def fit_filter(self, sample, H_dut, theta0, Nb, w_empirical_exp): + def fit_filter(self, sample, H_dut, theta0, Nb, w_empirical_exp, weights=None): # function to estimate filter H_ref = sample @@ -886,7 +890,7 @@ class CocalMethods: res = opt.minimize( self.evaluate_filter, x0=theta0, - args=(Nb, H_empirical, w_empirical_exp), + args=(Nb, H_empirical, w_empirical_exp, weights), method="Powell", ) @@ -898,14 +902,18 @@ class CocalMethods: return res.x - def evaluate_filter(self, theta, Nb, H_empirical, w_empirical_exp): + def evaluate_filter(self, theta, Nb, H_empirical, w_empirical_exp, weights=None): # evaluation function used in the numerical optimization inside fit_filter # H_ba = sig.freqz(b, a, w_empirical)[1] # quite slow H_ba = self.freqz_core(theta, Nb, w_empirical_exp) # faster - - # return np.linalg.norm(H_ba - H_empirical) # does not use log-scaling - # return np.linalg.norm(np.log(H_ba) - np.log(H_empirical)) - return np.linalg.norm( - np.log(np.abs(H_ba)) - np.log(np.abs(H_empirical)) - ) # only amplitude matters, phase discarded + + # compute difference + # diff = np.linalg.norm(H_ba - H_empirical) # does not use log-scaling + # diff = np.linalg.norm(np.log(H_ba) - np.log(H_empirical)) + diff = np.log(np.abs(H_ba)) - np.log(np.abs(H_empirical)) # only amplitude matters, phase discarded + + if weights is None: + return np.linalg.norm(diff) + else: + return np.linalg.norm(diff*weights)