From 363971e2c18574d1230a9a98a52ba37c7909cb13 Mon Sep 17 00:00:00 2001 From: linarphy Date: Mon, 27 May 2024 18:05:47 +0200 Subject: [PATCH] Use science_signal --- src/backscattering_analyzer/__init__.py | 64 ++++--------- src/backscattering_analyzer/analyzer.py | 121 +++++++++--------------- 2 files changed, 65 insertions(+), 120 deletions(-) diff --git a/src/backscattering_analyzer/__init__.py b/src/backscattering_analyzer/__init__.py index 04c6029..c27c8ae 100644 --- a/src/backscattering_analyzer/__init__.py +++ b/src/backscattering_analyzer/__init__.py @@ -1,34 +1,7 @@ -from __future__ import annotations +from science_signal.signal import Signal +from science_signal import interpolate, interpolate_abciss from numpy.typing import NDArray -from numpy import arange, float64, zeros, pi, sqrt - - -def interpolate_abciss(signals: tuple[Signal, ...]) -> NDArray[float64]: - """ - return the axis that would be used by the interpolate function - """ - rates: list[float] = [signal.rate for signal in signals] - - start: float = max([signal.x[0] for signal in signals]) - end: float = min([signal.x[-1] for signal in signals]) - - return arange(start, end, 1 / max(rates)) - - -def interpolate(signals: tuple[Signal, ...]) -> tuple[Signal, ...]: - """ - Interpolate multiple signals with a single abciss list, which has - the smallest interval and the bigget rate - """ - splines = [signal.spline() for signal in signals] - - x = interpolate_abciss(signals) - - new_signals = [ - Signal(x, spline(x), signals[0].settings) for spline in splines - ] - - return tuple(new_signals) +from numpy import float64, zeros, pi, sqrt def compute_light( @@ -43,7 +16,7 @@ def compute_light( Compute the projection from a given coupling and movement """ frequencies = interpolate_abciss( - (movement.sin().psd().sqrt(), coupling[0].abs()) + movement.sin().psd().sqrt(), abs(coupling[0]) ) parts = zeros( ( @@ -56,25 +29,28 @@ def compute_light( phase = (index + 1) * 4 * pi / wavelength factor_n = (movement * phase).sin().psd().sqrt() - coupling_n = coupling[0].abs() + coupling_n = abs(coupling[0]) factor_d = (movement * phase).cos().psd().sqrt() - coupling_d = coupling[1].abs() + coupling_d = abs(coupling[1]) factor_n, coupling_n, factor_d, coupling_d = interpolate( - (factor_n, coupling_n, factor_d, coupling_d) + 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 + parts[index] = opt_compute_light( + scatter_factor[index], + factor_n, + coupling_n, + factor_d, + coupling_d, + power_in, + power_out, + phase, ) return Signal( frequencies, sum(parts), - movement.settings, ) @@ -86,12 +62,14 @@ def opt_compute_light( coupling_d: Signal, power_in: float, power_out: float, + phase: float, ) -> NDArray[float64]: """ Optimized computing of light with pre-computed factor """ return ( sqrt(scatter_factor) + / phase * power_in / power_out * (coupling_n * factor_n + coupling_d * factor_d).y @@ -108,6 +86,7 @@ def fit_compute_light( power_out: float, data: Signal, reference: Signal, + phase: float, ) -> float: """ Scalar function used to find the right scattering factor @@ -124,13 +103,10 @@ def fit_compute_light( coupling_d=coupling_d, power_in=power_in, power_out=power_out, + phase=phase, ), - factor_n.settings, ) + reference - data ).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 0ef5340..d3e94f8 100644 --- a/src/backscattering_analyzer/analyzer.py +++ b/src/backscattering_analyzer/analyzer.py @@ -1,15 +1,16 @@ from sys import argv + +from scipy.interpolate import CubicSpline from backscattering_analyzer.settings import Settings -from backscattering_analyzer.signal import Signal +from science_signal.signal import Signal +from science_signal import interpolate from backscattering_analyzer import ( compute_light, fit_compute_light, opt_compute_light, - interpolate, ) -from numpy import loadtxt, logspace, pi +from numpy import argmin, loadtxt, logspace, pi, array from scipy.io.matlab import loadmat -from scipy.optimize import Bounds, minimize class Analyzer: @@ -53,9 +54,7 @@ class Analyzer: self.settings.log("loading bench movement") try: data = loadtxt(file).T - self.bench_movement = ( - Signal(data[0], data[1], self.settings) * 1e-6 - ) # um + self.bench_movement = Signal(data[0], data[1]) * 1e-6 # um except OSError: raise Exception("{file} does not exist".format(file=file)) @@ -67,9 +66,7 @@ class Analyzer: self.settings.log("loading mirror movement") try: data = loadtxt(file).T - self.mirror_movement = ( - Signal(data[0], data[1], self.settings) * 1e-6 - ) + self.mirror_movement = Signal(data[0], data[1]) * 1e-6 except OSError: raise Exception("{file} does not exist".format(file=file)) @@ -81,7 +78,7 @@ class Analyzer: self.settings.log("loading excited h(t)") try: data = loadtxt(file).T - self.data_signal = Signal(data[0], data[1], self.settings) + self.data_signal = Signal(data[0], data[1]) except OSError: raise Exception("{file} does not exist".format(file=file)) @@ -93,9 +90,7 @@ class Analyzer: self.settings.log("loading reference h(t)") try: data = loadtxt(file).T - self.reference_signal = Signal( - data[0], data[1], self.settings - ) + self.reference_signal = Signal(data[0], data[1]) except OSError: raise Exception("{file} does not exist".format(file=file)) @@ -114,7 +109,6 @@ class Analyzer: Signal( self.modelisation["freq"][0], coupling, - self.settings, ) for coupling in coupling_values ] @@ -134,13 +128,14 @@ class Analyzer: - self.mirror_movement * self.settings.calib_mirror ) .detrend("linear") - .low_pass_filter(5 * self.settings.vibration_frequency) + .filter(end=5 * self.settings.vibration_frequency) ) def compute_light(self) -> None: """ Compute projection with current bench excitation """ + self.settings.log("computing backscatterd light projection") self.projection = compute_light( scatter_factor=self.settings.scattering_factor, coupling=self.coupling, @@ -159,69 +154,64 @@ class Analyzer: guess = self.settings.scattering_factor[0] phase = 4 * pi / self.settings.wavelength factor_n = (self.movement * phase).sin().psd().sqrt() - coupling_n = self.coupling[0].abs() + coupling_n = abs(self.coupling[0]) factor_d = (self.movement * phase).cos().psd().sqrt() - coupling_d = self.coupling[1].abs() + coupling_d = abs(self.coupling[1]) - coupling_d = coupling_d.cut_x( - 10, 200 - ) # cut signal between 10 and 200 Hz + coupling_d = coupling_d.cut( + 11, + 25, + ) # cut signal between 10 and 200 Hz (updated to 11-15 to test) 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(), - ) - ) - ) - - bounds = Bounds(1e-30, 1e-3) - min_result = minimize( - fit_compute_light, - guess, - ( factor_n, coupling_n, factor_d, coupling_d, - self.settings.power_in, - self.settings.power_out, - data, - reference, - ), - method="Powell", - bounds=bounds, - ) - - if not min_result.success: - raise Exception(min_result.message) - - self.settings.log( - "found the best scattering factor in {} iterations".format( - min_result.nit + self.data_signal.psd().sqrt(), + self.reference_signal.psd().sqrt(), ) ) + x = logspace(-15, -3, 1000) + y = array( + [ + fit_compute_light( + x[i], + factor_n, + coupling_n, + factor_d, + coupling_d, + self.settings.power_in, + self.settings.power_out, + data, + reference, + phase, + ) + for i in range(len(x)) + ] + ) + + factor: float = x[argmin(y)] + import matplotlib.pyplot as plt projection = Signal( factor_n.x, opt_compute_light( - scatter_factor=min_result.x, + scatter_factor=factor, 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, + phase=phase, ), - self.settings, ) + _ = plt.loglog(x, y) + _ = plt.show() _ = plt.loglog(projection.x, projection.y, label="projection") _ = plt.loglog(reference.x, reference.y, label="référence") @@ -232,25 +222,4 @@ class Analyzer: _ = plt.legend() _ = plt.show() - """ - as the minimization does not work (why ?), visual help - """ - x = logspace(-50, 0, 1000) - y = [ - fit_compute_light( - x[i], - factor_n, - coupling_n, - factor_d, - coupling_d, - self.settings.power_in, - self.settings.power_out, - data, - reference, - ) - for i in range(len(x)) - ] - _ = plt.loglog(x, y) - _ = plt.show() - - return min_result.x + return factor