From 9070bb10e8e9f29e14c361c98f7d9eac6f11033e Mon Sep 17 00:00:00 2001
From: linarphy <linarphy@linarphy.net>
Date: Tue, 14 May 2024 15:19:44 +0200
Subject: [PATCH] Add psd math
---
src/backscattering_analyzer/analyzer.py | 138 ++++++++++++++++--------
src/backscattering_analyzer/settings.py | 7 +-
2 files changed, 102 insertions(+), 43 deletions(-)
diff --git a/src/backscattering_analyzer/analyzer.py b/src/backscattering_analyzer/analyzer.py
index bcced61..c576f63 100644
--- a/src/backscattering_analyzer/analyzer.py
+++ b/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(coupling.shape[1])
- scattering_factor = [0, 0]
-
- for index in range(len(coupling)):
- argument = (
- self.movement[1]
- * 4
- * pi
- * (index + 1)
- / self.settings.wavelength
- )
-
- result += (
- scattering_factor[index]
- * (self.settings.wavelength / 4 * pi) ** 2
- * (
- (
+ result = zeros(
+ len(
+ self.interpolate_abciss(
+ [
psd(
- sin(argument),
- fs=rate,
+ self.movement[1],
+ fs=freq_sample,
nperseg=nperseg,
- )[1]
- * abs(coupling[0])
- )
- ** 2
- + (
- psd(
- cos(argument),
- fs=rate,
- nperseg=nperseg,
- )[1]
- * abs(coupling[1])
- )
- ** 2
+ )[0],
+ self.modelisation["freq"][0],
+ ]
)
)
- return result
+ )
+
+ frequencies = 0
+
+ for index in range(len(coupling)):
+ phase = (index + 1) * 4 * pi / self.settings.wavelength
+
+ factor_n = array(
+ [
+ *psd(
+ 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 += (
+ sqrt(self.settings.scattering_factor[index])
+ * self.settings.power_in
+ / self.settings.power_out
+ * (
+ coupling_n[1] * factor_n[1]
+ + coupling_d[1] * factor_d[1]
+ )
+ )
+ return array(
+ [
+ frequencies,
+ result,
+ ]
+ )
diff --git a/src/backscattering_analyzer/settings.py b/src/backscattering_analyzer/settings.py
index a5735ad..13c4a2c 100644
--- a/src/backscattering_analyzer/settings.py
+++ b/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):