wip: first draft of plot showing compensated transfer function

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()