diff --git a/app/cocal_methods.py b/app/cocal_methods.py index 755c671b50702199cd404fc0b1532245f15fbed2..6ece6159b00d833c51daf1d3566d7baa1a8d7ee7 100644 --- a/app/cocal_methods.py +++ b/app/cocal_methods.py @@ -25,6 +25,7 @@ import xmlschema from pydub import AudioSegment from PyDynamic.misc.tools import complex_2_real_imag, real_imag_2_complex from PyDynamic.uncertainty.propagate_MonteCarlo import UMC_generic +from PyDynamic.uncertainty.propagate_DFT import DFT_multiply, DFT_deconv from scipy.interpolate import interp1d from sounddevice import InputStream, play, query_devices from soundfile import SoundFile @@ -378,7 +379,6 @@ class CocalMethods: return tree def generate_spectral_domain_visualization(self): - # provide shortcuts for plotting b = self.transfer_behavior["IIR"]["b"] a = self.transfer_behavior["IIR"]["a"] @@ -387,18 +387,25 @@ class CocalMethods: # further relevant parameters time_delta = np.mean(np.diff(self.ref_times[0])) frame_rate = 1.0 / time_delta - + Nb = len(b) ba = np.concatenate((b, a[1:])) mask = np.logical_and(self.ref_frequency < 5000.0, self.ref_frequency >= 0.0) + mask = np.logical_and(self.ref_frequency < 20000.0, self.ref_frequency >= 0.0) # visualize coefficient uncertainties by plotting the spectrum uncertainty via MC def draw_samples(size): return np.random.multivariate_normal(ba, ba_cov, size=size) - w = np.linspace(0, np.pi, 50) - f = w / (2 * np.pi) * frame_rate - w_exp = np.exp(-1j * w) # ??? + ## reduced frequency resolution + # w = np.linspace(0, np.pi, 50) + # f = w / (2 * np.pi) * frame_rate + # w_exp = np.exp(-1j * w) # ??? + + # full frequency resolution: + f = self.ref_frequency + w = self.ref_frequency / frame_rate * (2 * np.pi) + w_exp = np.exp(-1j * w) def evaluate(sample): return complex_2_real_imag(self.freqz_core(sample, Nb, w_exp)) @@ -406,7 +413,7 @@ class CocalMethods: umc_kwargs = { "draw_samples": draw_samples, "evaluate": evaluate, - "runs": 20, + "runs": 200, "blocksize": 8, "n_cpu": 1, "compute_full_covariance": True, @@ -417,29 +424,69 @@ class CocalMethods: h = real_imag_2_complex(h_ri) h_unc = real_imag_2_complex(np.sqrt(np.diag(h_cov))) + h_empirical = self.dut_spectrum / self.ref_spectrum + + # compensate + h_comp_ri, h_comp_cov = DFT_deconv( + complex_2_real_imag(h_empirical), + h_ri, + np.zeros((2 * len(h_empirical), 2 * len(h_empirical))), + h_cov, + ) + h_comp = real_imag_2_complex(h_comp_ri) + h_comp_unc = real_imag_2_complex(np.sqrt(np.diag(h_comp_cov))) + + # regularize + #regularized_filter = sig.cheby2() + #h_regularized = self.freqz_core() + #h_compensated_and_regularized = DFT_multiply( + # h_compensated, h_regularized, h_regularized_co + #) + # visualize result - fig, ax = plt.subplots(2, 1, sharex=True) - ax[0].plot(f, np.abs(h), label="fitted TF") - ax[0].fill_between( - f, - np.abs(h) - np.abs(h_unc), - np.abs(h) + np.abs(h_unc), + fig, ax = plt.subplots(1, 1, sharex=True, squeeze=False) + ax[0,0].plot( + f[mask], + np.abs(h[mask]), + label="fitted TF", + color="blue", + ) + ax[0,0].fill_between( + f[mask], + np.abs(h[mask]) - np.abs(h_unc[mask]), + np.abs(h[mask]) + np.abs(h_unc[mask]), alpha=0.3, label="unc of fitted TF", + color="blue", ) - ax[0].scatter( + ax[0,0].scatter( self.ref_frequency[mask], - np.abs(self.dut_spectrum[mask] / self.ref_spectrum[mask]), + np.abs(h_empirical[mask]), label="empirical TF", + s=2, + color="black", + ) + ax[0,0].scatter( + self.ref_frequency[mask], + np.abs(h_comp[mask]), + label="compensated TF", s=1, + color="red", ) - ax[0].plot(f, np.ones_like(f), "--r", label="ideal") - ax[0].plot(f, np.abs(1.0/h), label="compensated empirical") - ax[0].legend() - ax[0].set_xscale("log") - ax[0].set_yscale("log") - ax[1].plot(f, h_unc) - #plt.savefig(self.result_image_path) + ax[0,0].fill_between( + f[mask], + np.abs(h_comp[mask]) - np.abs(h_comp_unc[mask]), + np.abs(h_comp[mask]) + np.abs(h_comp_unc[mask]), + alpha=0.3, + label="unc of compensated TF", + color="red", + ) + + ax[0,0].plot(f, np.ones_like(f), "--r", label="ideal") + ax[0,0].legend() + ax[0,0].set_xscale("log") + ax[0,0].set_yscale("log") + # plt.savefig(self.result_image_path) plt.show()