Add psd math

src/backscattering_analyzer/analyzer.py

@@ -10,7 +10,7 @@ from numpy import loadtxt, array
 from scipy.io.matlab import loadmat
 
 # maths
-from numpy import mean, zeros, pi, sin, cos, arange
+from numpy import mean, zeros, pi, sin, cos, arange, sqrt
 from scipy.signal import welch as psd
 from scipy.interpolate import CubicSpline
 
@@ -121,7 +121,9 @@ class Analyzer:
         self.movement = array(
             [
                 self.bench_movement[0],
-                self.bench_movement[1] - self.mirror_movement[1],
+                (self.settings.calib_bench * self.bench_movement[1])
+                    - (self.settings.calib_mirror
+                        * self.mirror_movement[1]),
             ]
         )
         self.movement[1] -= mean(self.movement[1])
@@ -161,20 +163,29 @@ class Analyzer:
             ]
         )
 
+    def interpolate_abciss(self, abcisses):
+        """
+        Return the axis that would be used by the interpolate method
+        """
+        rates = 1 / array(
+            [abciss[1] - abciss[0] for abciss in abcisses]
+        )
+        start = max([abciss[0] for abciss in abcisses])
+        end = min([abciss[-1] for abciss in abcisses])
+
+        return arange(start, end, 1 / max(rates))
+
     def interpolate(self, signals):
         """
         Interpolate multiple signals with a single abciss list, which
         has the smallest interval and the biggest rate
         """
-        rates = 1 / (signals[:, 0, 1] - signals[:, 0, 0])
-        start = max(signals[:, 0, 0])
-        end = min(signals[:, 0, -1])
 
         splines = [
             CubicSpline(signal[0], signal[1]) for signal in signals
         ]
 
-        x = arange(start, end, 1 / max(rates))
+        x = self.interpolate_abciss([signal[0] for signal in signals])
 
         signals = [array([x, spline(x)]) for spline in splines]
 
@@ -187,44 +198,87 @@ class Analyzer:
         """
         coupling = self.modelisation[self.settings.coupling_name()]
 
-        nperseg = int(len(coupling[0]) * 2) - 1
+        freq_sample = 1 / (self.movement[0, 1] - self.movement[0, 0])
+        nperseg = (
+            int(2 * len(coupling[0])) - 1 + len(coupling[0])
+        )  # minimum to keep same length than coupling, but not for
+        # the same frequency range, so more
 
-        rate = 1 / (self.movement[0, 1] - self.movement[0, 0])
+        result = zeros(
+            len(
+                self.interpolate_abciss(
+                    [
+                        psd(
+                            self.movement[1],
+                            fs=freq_sample,
+                            nperseg=nperseg,
+                        )[0],
+                        self.modelisation["freq"][0],
+                    ]
+                )
+            )
+        )
 
-        result = zeros(coupling.shape[1])
+        frequencies = 0
-        scattering_factor = [0, 0]
 
         for index in range(len(coupling)):
-            argument = (
+            phase = (index + 1) * 4 * pi / self.settings.wavelength
-                self.movement[1]
+
-                * 4
+            factor_n = array(
-                * pi
+                [
-                * (index + 1)
+                    *psd(
-                / self.settings.wavelength
+                        sin(phase * self.movement[1]),
+                        fs=freq_sample,
+                        nperseg=nperseg,
+                    )
+                ]
+            )
+            coupling_n = array(
+                [
+                    self.modelisation["freq"][0],
+                    abs(coupling[0]),
+                ]
+            )
+            factor_d = array(
+                [
+                    *psd(
+                        cos(phase * self.movement[1]),
+                        fs=freq_sample,
+                        nperseg=nperseg,
+                    )
+                ]
+            )
+            coupling_d = array(
+                [
+                    self.modelisation["freq"][0],
+                    abs(coupling[1]),
+                ]
             )
 
+            unpack = self.interpolate(
+                [factor_n, coupling_n, factor_d, coupling_d]
+            )
+            factor_n, coupling_n, factor_d, coupling_d = (
+                unpack[0],
+                unpack[1],
+                unpack[2],
+                unpack[3],
+            )
+
+            frequencies = factor_n[0]
+
             result += (
-                scattering_factor[index]
+                sqrt(self.settings.scattering_factor[index])
-                * (self.settings.wavelength / 4 * pi) ** 2
+                * self.settings.power_in
+                / self.settings.power_out
                 * (
-                    (
+                    coupling_n[1] * factor_n[1]
-                        psd(
+                    + coupling_d[1] * factor_d[1]
-                            sin(argument),
-                            fs=rate,
-                            nperseg=nperseg,
-                        )[1]
-                        * abs(coupling[0])
-                    )
-                    ** 2
-                    + (
-                        psd(
-                            cos(argument),
-                            fs=rate,
-                            nperseg=nperseg,
-                        )[1]
-                        * abs(coupling[1])
-                    )
-                    ** 2
                 )
             )
-        return result
+        return array(
+            [
+                frequencies,
+                result,
+            ]
+        )

src/backscattering_analyzer/settings.py

@@ -17,7 +17,12 @@ class Settings:
         self.date = "2023_03_24"
         self.folder = Path("/home/demagny/data")
         self.modelisation = "scatterCouplingO4.mat"
-        self.wavelength = 1.064e-6
+        self.calib_bench = 1.15
+        self.calib_mirror = 1.15
+        self.wavelength = 1.064e-6  # m
+        self.power_in = 23  # W
+        self.power_out = 8e-3  # W
+        self.scattering_factor = [1e-17, 0]  # parameter to change
 
         index = 0
         while index < len(options):