linarphy/backscattering-analyzer
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Use science_signal

Browse source
This commit is contained in:
linarphy 2024-05-27 18:05:47 +02:00
parent 002c7ecd9d
commit 363971e2c1
No known key found for this signature in database
GPG key ID: E61920135EFF2295
2 changed files with 65 additions and 120 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

64
src/backscattering_analyzer/__init__.py
View file

@ -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.typing import NDArray
from numpy import arange, float64, zeros, pi, sqrt from numpy import 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)
def compute_light( def compute_light(
@ -43,7 +16,7 @@ def compute_light(
Compute the projection from a given coupling and movement Compute the projection from a given coupling and movement
""" """
frequencies = interpolate_abciss( frequencies = interpolate_abciss(
(movement.sin().psd().sqrt(), coupling[0].abs()) movement.sin().psd().sqrt(), abs(coupling[0])
) )
parts = zeros( parts = zeros(
( (
@ -56,25 +29,28 @@ def compute_light(
phase = (index + 1) * 4 * pi / wavelength phase = (index + 1) * 4 * pi / wavelength
factor_n = (movement * phase).sin().psd().sqrt() factor_n = (movement * phase).sin().psd().sqrt()
coupling_n = coupling[0].abs() coupling_n = abs(coupling[0])
factor_d = (movement * phase).cos().psd().sqrt() 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 = interpolate(
(factor_n, coupling_n, factor_d, coupling_d) factor_n, coupling_n, factor_d, coupling_d
) )
parts[index] = ( parts[index] = opt_compute_light(
sqrt(scatter_factor[index]) scatter_factor[index],
* power_in factor_n,
/ power_out coupling_n,
* (coupling_n * factor_n + coupling_d * factor_d).y factor_d,
coupling_d,
power_in,
power_out,
phase,
) )
return Signal( return Signal(
frequencies, frequencies,
sum(parts), sum(parts),
movement.settings,
) )
@ -86,12 +62,14 @@ def opt_compute_light(
coupling_d: Signal, coupling_d: Signal,
power_in: float, power_in: float,
power_out: float, power_out: float,
phase: float,
) -> NDArray[float64]: ) -> NDArray[float64]:
""" """
Optimized computing of light with pre-computed factor Optimized computing of light with pre-computed factor
""" """
return ( return (
sqrt(scatter_factor) sqrt(scatter_factor)
/ phase
* power_in * power_in
/ power_out / power_out
* (coupling_n * factor_n + coupling_d * factor_d).y * (coupling_n * factor_n + coupling_d * factor_d).y
@ -108,6 +86,7 @@ def fit_compute_light(
power_out: float, power_out: float,
data: Signal, data: Signal,
reference: Signal, reference: Signal,
phase: float,
) -> float: ) -> float:
""" """
Scalar function used to find the right scattering factor Scalar function used to find the right scattering factor
@ -124,13 +103,10 @@ def fit_compute_light(
coupling_d=coupling_d, coupling_d=coupling_d,
power_in=power_in, power_in=power_in,
power_out=power_out, power_out=power_out,
phase=phase,
), ),
factor_n.settings,
) )
+ reference + reference
- data - data
).y ).y
) )
from backscattering_analyzer.signal import Signal # no circular import

131
src/backscattering_analyzer/analyzer.py
View file

@ -1,15 +1,16 @@
from sys import argv from sys import argv
from scipy.interpolate import CubicSpline
from backscattering_analyzer.settings import Settings 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 ( from backscattering_analyzer import (
compute_light, compute_light,
fit_compute_light, fit_compute_light,
opt_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.io.matlab import loadmat
from scipy.optimize import Bounds, minimize
class Analyzer: class Analyzer:
@ -53,9 +54,7 @@ class Analyzer:
self.settings.log("loading bench movement") self.settings.log("loading bench movement")
try: try:
data = loadtxt(file).T data = loadtxt(file).T
self.bench_movement = ( self.bench_movement = Signal(data[0], data[1]) * 1e-6 # um
Signal(data[0], data[1], self.settings) * 1e-6
) # um
except OSError: except OSError:
raise Exception("{file} does not exist".format(file=file)) raise Exception("{file} does not exist".format(file=file))
@ -67,9 +66,7 @@ class Analyzer:
self.settings.log("loading mirror movement") self.settings.log("loading mirror movement")
try: try:
data = loadtxt(file).T data = loadtxt(file).T
self.mirror_movement = ( self.mirror_movement = Signal(data[0], data[1]) * 1e-6
Signal(data[0], data[1], self.settings) * 1e-6
)
except OSError: except OSError:
raise Exception("{file} does not exist".format(file=file)) raise Exception("{file} does not exist".format(file=file))
@ -81,7 +78,7 @@ class Analyzer:
self.settings.log("loading excited h(t)") self.settings.log("loading excited h(t)")
try: try:
data = loadtxt(file).T data = loadtxt(file).T
self.data_signal = Signal(data[0], data[1], self.settings) self.data_signal = Signal(data[0], data[1])
except OSError: except OSError:
raise Exception("{file} does not exist".format(file=file)) raise Exception("{file} does not exist".format(file=file))
@ -93,9 +90,7 @@ class Analyzer:
self.settings.log("loading reference h(t)") self.settings.log("loading reference h(t)")
try: try:
data = loadtxt(file).T data = loadtxt(file).T
self.reference_signal = Signal( self.reference_signal = Signal(data[0], data[1])
data[0], data[1], self.settings
)
except OSError: except OSError:
raise Exception("{file} does not exist".format(file=file)) raise Exception("{file} does not exist".format(file=file))
@ -114,7 +109,6 @@ class Analyzer:
Signal( Signal(
self.modelisation["freq"][0], self.modelisation["freq"][0],
coupling, coupling,
self.settings,
) )
for coupling in coupling_values for coupling in coupling_values
] ]
@ -134,13 +128,14 @@ class Analyzer:
- self.mirror_movement * self.settings.calib_mirror - self.mirror_movement * self.settings.calib_mirror
) )
.detrend("linear") .detrend("linear")
.low_pass_filter(5 * self.settings.vibration_frequency) .filter(end=5 * self.settings.vibration_frequency)
) )
def compute_light(self) -> None: def compute_light(self) -> None:
""" """
Compute projection with current bench excitation Compute projection with current bench excitation
""" """
self.settings.log("computing backscatterd light projection")
self.projection = compute_light( self.projection = compute_light(
scatter_factor=self.settings.scattering_factor, scatter_factor=self.settings.scattering_factor,
coupling=self.coupling, coupling=self.coupling,
@ -159,17 +154,17 @@ class Analyzer:
guess = self.settings.scattering_factor[0] guess = self.settings.scattering_factor[0]
phase = 4 * pi / self.settings.wavelength phase = 4 * pi / self.settings.wavelength
factor_n = (self.movement * phase).sin().psd().sqrt() 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() 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( coupling_d = coupling_d.cut(
10, 200 11,
) # cut signal between 10 and 200 Hz 25,
) # cut signal between 10 and 200 Hz (updated to 11-15 to test)
factor_n, coupling_n, factor_d, coupling_d, data, reference = ( factor_n, coupling_n, factor_d, coupling_d, data, reference = (
interpolate( interpolate(
(
factor_n, factor_n,
coupling_n, coupling_n,
factor_d, factor_d,
@ -178,65 +173,10 @@ class Analyzer:
self.reference_signal.psd().sqrt(), self.reference_signal.psd().sqrt(),
) )
) )
)
bounds = Bounds(1e-30, 1e-3) x = logspace(-15, -3, 1000)
min_result = minimize( y = array(
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
)
)
import matplotlib.pyplot as plt
projection = Signal(
factor_n.x,
opt_compute_light(
scatter_factor=min_result.x,
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,
),
self.settings,
)
_ = plt.loglog(projection.x, projection.y, label="projection")
_ = plt.loglog(reference.x, reference.y, label="référence")
_ = plt.loglog(data.x, data.y, label="data")
_ = plt.loglog(
reference.x, projection.y + reference.y, label="somme"
)
_ = plt.legend()
_ = plt.show()
"""
as the minimization does not work (why ?), visual help
"""
x = logspace(-50, 0, 1000)
y = [
fit_compute_light( fit_compute_light(
x[i], x[i],
factor_n, factor_n,
@ -247,10 +187,39 @@ class Analyzer:
self.settings.power_out, self.settings.power_out,
data, data,
reference, reference,
phase,
) )
for i in range(len(x)) 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=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,
),
)
_ = plt.loglog(x, y) _ = plt.loglog(x, y)
_ = plt.show() _ = plt.show()
return min_result.x _ = plt.loglog(projection.x, projection.y, label="projection")
_ = plt.loglog(reference.x, reference.y, label="référence")
_ = plt.loglog(data.x, data.y, label="data")
_ = plt.loglog(
reference.x, projection.y + reference.y, label="somme"
)
_ = plt.legend()
_ = plt.show()
return factor