From ec5e67536e00cd9b95292d7d402b49c0b57028b1 Mon Sep 17 00:00:00 2001 From: linarphy Date: Thu, 23 May 2024 18:37:40 +0200 Subject: [PATCH] Add fit --- src/backscattering_analyzer/__init__.py | 70 ++++++++++++++- src/backscattering_analyzer/analyzer.py | 110 +++++++++++++++--------- src/backscattering_analyzer/signal.py | 5 +- 3 files changed, 139 insertions(+), 46 deletions(-) diff --git a/src/backscattering_analyzer/__init__.py b/src/backscattering_analyzer/__init__.py index ee1aa3e..adb422c 100644 --- a/src/backscattering_analyzer/__init__.py +++ b/src/backscattering_analyzer/__init__.py @@ -1,7 +1,6 @@ from __future__ import annotations -from typing import Any from numpy.typing import NDArray -from numpy import arange, float64 +from numpy import arange, float64, zeros, pi, sqrt def interpolate_abciss(signals: tuple[Signal, ...]) -> NDArray[float64]: @@ -32,4 +31,71 @@ def interpolate(signals: tuple[Signal, ...]) -> tuple[Signal, ...]: return tuple(new_signals) +def compute_light( + scatter_factor: list[float], + coupling: list[Signal], + movement: Signal, + wavelength: float, + power_in: float, + power_out: float, +) -> Signal: + """ + Compute the projection from a given coupling and movement + """ + frequencies = interpolate_abciss( + (movement.sin().psd().sqrt(), coupling[0].abs()) + ) + parts = zeros( + ( + len(coupling), + len(frequencies), + ) + ) + + for index in range(2): + phase = (index + 1) * 4 * pi / wavelength + + factor_n = (movement * phase).sin().psd().sqrt() + coupling_n = coupling[0].abs() + factor_d = (movement * phase).cos().psd().sqrt() + coupling_d = coupling[1].abs() + + factor_n, coupling_n, factor_d, coupling_d = interpolate( + (factor_n, coupling_n, factor_d, coupling_d) + ) + + parts[index] = ( + sqrt(scatter_factor[index]) + * power_in + / power_out + * (coupling_n * factor_n + coupling_d * factor_d).y + ) + + return Signal( + frequencies, + sum(parts), + movement.settings, + ) + + +def opt_compute_light( + scatter_factor: float, + factor_n: Signal, + coupling_n: Signal, + factor_d: Signal, + coupling_d: Signal, + power_in: float, + power_out: float, +) -> NDArray[float64]: + """ + Optimized computing of light with pre-computed factor + """ + return ( + sqrt(scatter_factor) + * power_in + / power_out + * (coupling_n * factor_n + coupling_d * factor_d).y + ) + + from backscattering_analyzer.signal import Signal # no circular import diff --git a/src/backscattering_analyzer/analyzer.py b/src/backscattering_analyzer/analyzer.py index be0070c..c674e35 100644 --- a/src/backscattering_analyzer/analyzer.py +++ b/src/backscattering_analyzer/analyzer.py @@ -2,20 +2,23 @@ from sys import argv from backscattering_analyzer.settings import Settings from backscattering_analyzer.signal import Signal -from backscattering_analyzer import interpolate, interpolate_abciss -from numpy import loadtxt +from backscattering_analyzer import ( + compute_light, + opt_compute_light, + interpolate, +) +from numpy import loadtxt, logspace, where, zeros, argmin, intp, pi from scipy.io.matlab import loadmat -# maths -from numpy import zeros, pi, sqrt - class Analyzer: """ Utility class to study backscattering light in VIRGO """ - def __init__(self, arguments: None | list[str] | str = None) -> None: + def __init__( + self, arguments: None | list[str] | str = None + ) -> None: if arguments is None: options = [] elif isinstance(arguments, str): @@ -135,47 +138,68 @@ class Analyzer: def compute_light(self) -> None: """ - Compute psd of the computed projection with current bench - excitation + Compute projection with current bench excitation """ - results = zeros( - ( - len(self.coupling), - len( - interpolate_abciss( - (self.movement.psd(), self.coupling[0].abs()) - ) - ), + self.projection = compute_light( + scatter_factor=self.settings.scattering_factor, + coupling=self.coupling, + movement=self.movement, + wavelength=self.settings.wavelength, + power_in=self.settings.power_in, + power_out=self.settings.power_out, + ) + + def fit_scatter_factor( + self, start: int, stop: int, number: int + ) -> tuple[intp, float]: + """ + Find the best scatter factor (first order only) in the given + range + """ + import matplotlib.pyplot as plt + factors = logspace(start, stop, number) + sums = zeros(number) + + phase = 4 * pi / self.settings.wavelength + factor_n = (self.movement * phase).sin().psd().sqrt() + coupling_n = self.coupling[0].abs() + factor_d = (self.movement * phase).cos().psd().sqrt() + coupling_d = self.coupling[1].abs() + + coupling_d.cut(5, 40) # cut signal between 5 and 40 Hz + + factor_n, coupling_n, factor_d, coupling_d, data, reference = ( + interpolate( + ( + factor_n, + coupling_n, + factor_d, + coupling_d, + self.data_signal.psd().sqrt(), + self.reference_signal.psd().sqrt(), + ) ) ) - # frequencies depends of psd result, which we do not have yet - frequencies = zeros(results.shape[1]) + reference = reference.y + data = data.y - for index in range(len(self.coupling)): - phase = (index + 1) * 4 * pi / self.settings.wavelength - - factor_n = (self.movement * phase).sin().psd().sqrt() - coupling_n = self.coupling[0].abs() - factor_d = (self.movement * phase).cos().psd().sqrt() - coupling_d = self.coupling[1].abs() - - factor_n, coupling_n, factor_d, coupling_d = interpolate( - (factor_n, coupling_n, factor_d, coupling_d) + for index in range(number): + self.settings.log("{}".format(index)) + projection = opt_compute_light( + scatter_factor=factors[index], + factor_n=factor_n, + coupling_n=coupling_n, + factor_d=factor_d, + coupling_d=coupling_d, + power_in=self.settings.power_in, + power_out=self.settings.power_out, ) + diff = abs(projection + reference - data) + _ = plt.loglog(projection + reference) + _ = plt.loglog(data) + _ = plt.show() + sums[index] = sum(diff) + min_index = argmin(sums) - # no need to redefine it each time but simpler here - frequencies = factor_n.x - - results[index] = ( - sqrt(self.settings.scattering_factor[index]) - * self.settings.power_in - / self.settings.power_out - * (coupling_n * factor_n + coupling_d * factor_d).y - ) - - self.projection = Signal( - frequencies, - sum(results), - self.settings, - ) + return min_index, factors[min_index] diff --git a/src/backscattering_analyzer/signal.py b/src/backscattering_analyzer/signal.py index 092bfc8..1829df8 100644 --- a/src/backscattering_analyzer/signal.py +++ b/src/backscattering_analyzer/signal.py @@ -15,7 +15,10 @@ class Signal: """ def __init__( - self, x: NDArray[float64], value: NDArray[float64], settings: Settings + self, + x: NDArray[float64], + value: NDArray[float64], + settings: Settings, ) -> None: if x.shape != value.shape: raise Exception("x and y does not have the same dimension")