Source code for cdl.algorithms.signal

# Copyright (c) DataLab Platform Developers, BSD 3-Clause license, see LICENSE file.

"""
.. Signal Processing Algorithms (see parent package :mod:`cdl.algorithms`)
"""

# pylint: disable=invalid-name  # Allows short reference names like x, y, ...

from __future__ import annotations

import abc
import warnings
from typing import Literal

import numpy as np
import scipy.interpolate
import scipy.optimize
import scipy.signal
import scipy.special

# MARK: Level adjustment ---------------------------------------------------------------


[docs] def normalize( yin: np.ndarray, parameter: Literal["maximum", "amplitude", "area", "energy", "rms"] = "maximum", ) -> np.ndarray: """Normalize input array to a given parameter. Args: yin: Input array parameter: Normalization parameter. Defaults to "maximum" Returns: Normalized array """ axis = len(yin.shape) - 1 if parameter == "maximum": maximum = np.max(yin, axis) if axis == 1: maximum = maximum.reshape((len(maximum), 1)) maxarray = np.tile(maximum, yin.shape[axis]).reshape(yin.shape) return yin / maxarray if parameter == "amplitude": ytemp = np.array(yin, copy=True) minimum = np.min(yin, axis) if axis == 1: minimum = minimum.reshape((len(minimum), 1)) ytemp -= minimum return normalize(ytemp, parameter="maximum") if parameter == "area": return yin / yin.sum() if parameter == "energy": return yin / np.sqrt(np.sum(yin * yin.conjugate())) if parameter == "rms": return yin / np.sqrt(np.mean(yin * yin.conjugate())) raise RuntimeError(f"Unsupported parameter {parameter}")
# MARK: Fourier analysis ---------------------------------------------------------------
[docs] def fft1d( x: np.ndarray, y: np.ndarray, shift: bool = True ) -> tuple[np.ndarray, np.ndarray]: """Compute FFT on X,Y data. Args: x: X data y: Y data shift: Shift the zero frequency to the center of the spectrum. Defaults to True. Returns: X data, Y data (tuple) """ y1 = np.fft.fft(y) x1 = np.fft.fftfreq(x.size, d=x[1] - x[0]) if shift: x1 = np.fft.fftshift(x1) y1 = np.fft.fftshift(y1) return x1, y1
[docs] def ifft1d( x: np.ndarray, y: np.ndarray, shift: bool = True ) -> tuple[np.ndarray, np.ndarray]: """Compute iFFT on X,Y data. Args: x: X data y: Y data shift: Shift the zero frequency to the center of the spectrum. Defaults to True. Returns: X data, Y data (tuple) """ if shift: y = np.fft.ifftshift(y) y1 = np.fft.ifft(y) # Recalculate the original time domain array dt = 1.0 / (x[-1] - x[0] + (x[1] - x[0])) x1 = np.arange(y1.size) * dt return x1, y1.real
[docs] def magnitude_spectrum( x: np.ndarray, y: np.ndarray, log_scale: bool = False ) -> tuple[np.ndarray, np.ndarray]: """Compute magnitude spectrum. Args: x: X data y: Y data log_scale: Use log scale. Defaults to False. Returns: Magnitude spectrum (X data, Y data) """ x1, y1 = fft1d(x, y) if log_scale: y_mag = 20 * np.log10(np.abs(y1)) y_mag = np.abs(y1) return x1, y_mag
[docs] def phase_spectrum(x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray]: """Compute phase spectrum. Args: x: X data y: Y data Returns: Phase spectrum in degrees (X data, Y data) """ x1, y1 = fft1d(x, y) y_phase = np.rad2deg(np.angle(y1)) return x1, y_phase
[docs] def psd( x: np.ndarray, y: np.ndarray, log_scale: bool = False ) -> tuple[np.ndarray, np.ndarray]: """Compute Power Spectral Density (PSD), using the Welch method. Args: x: X data y: Y data log_scale: Use log scale. Defaults to False. Returns: Power Spectral Density (PSD): X data, Y data (tuple) """ x1, y1 = scipy.signal.welch(y, fs=sampling_rate(x)) if log_scale: y1 = 10 * np.log10(y1) return x1, y1
[docs] def sort_frequencies(x: np.ndarray, y: np.ndarray) -> np.ndarray: """Sort from X,Y data by computing FFT(y). Args: x: X data y: Y data Returns: Sorted frequencies in ascending order """ freqs, fourier = fft1d(x, y, shift=False) return freqs[np.argsort(fourier)]
# MARK: Peak detection -----------------------------------------------------------------
[docs] def peak_indices( y, thres: float = 0.3, min_dist: int = 1, thres_abs: bool = False ) -> np.ndarray: # Copyright (c) 2014 Lucas Hermann Negri # Unmodified code snippet from PeakUtils 1.3.0 """Peak detection routine. Finds the numeric index of the peaks in *y* by taking its first order difference. By using *thres* and *min_dist* parameters, it is possible to reduce the number of detected peaks. *y* must be signed. Parameters ---------- y : ndarray (signed) 1D amplitude data to search for peaks. thres : float between [0., 1.] Normalized threshold. Only the peaks with amplitude higher than the threshold will be detected. min_dist : int Minimum distance between each detected peak. The peak with the highest amplitude is preferred to satisfy this constraint. thres_abs: boolean If True, the thres value will be interpreted as an absolute value, instead of a normalized threshold. Returns ------- ndarray Array containing the numeric indices of the peaks that were detected """ if isinstance(y, np.ndarray) and np.issubdtype(y.dtype, np.unsignedinteger): raise ValueError("y must be signed") if not thres_abs: thres = thres * (np.max(y) - np.min(y)) + np.min(y) # compute first order difference dy = np.diff(y) # propagate left and right values successively to fill all plateau pixels # (0-value) (zeros,) = np.where(dy == 0) # check if the signal is totally flat if len(zeros) == len(y) - 1: return np.array([]) if len(zeros): # compute first order difference of zero indices zeros_diff = np.diff(zeros) # check when zeros are not chained together (zeros_diff_not_one,) = np.add(np.where(zeros_diff != 1), 1) # make an array of the chained zero indices zero_plateaus = np.split(zeros, zeros_diff_not_one) # fix if leftmost value in dy is zero if zero_plateaus[0][0] == 0: dy[zero_plateaus[0]] = dy[zero_plateaus[0][-1] + 1] zero_plateaus.pop(0) # fix if rightmost value of dy is zero if len(zero_plateaus) > 0 and zero_plateaus[-1][-1] == len(dy) - 1: dy[zero_plateaus[-1]] = dy[zero_plateaus[-1][0] - 1] zero_plateaus.pop(-1) # for each chain of zero indices for plateau in zero_plateaus: median = np.median(plateau) # set leftmost values to leftmost non zero values dy[plateau[plateau < median]] = dy[plateau[0] - 1] # set rightmost and middle values to rightmost non zero values dy[plateau[plateau >= median]] = dy[plateau[-1] + 1] # find the peaks by using the first order difference peaks = np.where( (np.hstack([dy, 0.0]) < 0.0) & (np.hstack([0.0, dy]) > 0.0) & (np.greater(y, thres)) )[0] # handle multiple peaks, respecting the minimum distance if peaks.size > 1 and min_dist > 1: highest = peaks[np.argsort(y[peaks])][::-1] rem = np.ones(y.size, dtype=bool) rem[peaks] = False for peak in highest: if not rem[peak]: sl = slice(max(0, peak - min_dist), peak + min_dist + 1) rem[sl] = True rem[peak] = False peaks = np.arange(y.size)[~rem] return peaks
[docs] def xpeak(x: np.ndarray, y: np.ndarray) -> float: """Return default peak X-position (assuming a single peak). Args: x: X data y: Y data Returns: Peak X-position """ peaks = peak_indices(y) if peaks.size == 1: return x[peaks[0]] return np.average(x, weights=y)
# MARK: Interpolation ------------------------------------------------------------------
[docs] def interpolate( x: np.ndarray, y: np.ndarray, xnew: np.ndarray, method: Literal["linear", "spline", "quadratic", "cubic", "barycentric", "pchip"], fill_value: float | None = None, ) -> np.ndarray: """Interpolate data. Args: x: X data y: Y data xnew: New X data method: Interpolation method fill_value: Fill value. Defaults to None. This value is used to fill in for requested points outside of the X data range. It is only used if the method argument is 'linear', 'cubic' or 'pchip'. Returns: Interpolated Y data """ interpolator_extrap = None if method == "linear": # Linear interpolation using NumPy's interp function: ynew = np.interp(xnew, x, y, left=fill_value, right=fill_value) elif method == "spline": # Spline using 1-D interpolation with SciPy's interpolate package: # pylint: disable=unbalanced-tuple-unpacking knots, coeffs, degree = scipy.interpolate.splrep(x, y, s=0) ynew = scipy.interpolate.splev(xnew, (knots, coeffs, degree), der=0) elif method == "quadratic": # Quadratic interpolation using NumPy's polyval function: coeffs = np.polyfit(x, y, 2) ynew = np.polyval(coeffs, xnew) elif method == "cubic": # Cubic interpolation using SciPy's Akima1DInterpolator class: interpolator_extrap = scipy.interpolate.Akima1DInterpolator(x, y) elif method == "barycentric": # Barycentric interpolation using SciPy's BarycentricInterpolator class: interpolator = scipy.interpolate.BarycentricInterpolator(x, y) ynew = interpolator(xnew) elif method == "pchip": # PCHIP interpolation using SciPy's PchipInterpolator class: interpolator_extrap = scipy.interpolate.PchipInterpolator(x, y) else: raise ValueError(f"Invalid interpolation method {method}") if interpolator_extrap is not None: ynew = interpolator_extrap(xnew, extrapolate=fill_value is None) if fill_value is not None: ynew[xnew < x[0]] = fill_value ynew[xnew > x[-1]] = fill_value return ynew
# MARK: Windowing ----------------------------------------------------------------------
[docs] def windowing( y: np.ndarray, method: Literal[ "barthann", "bartlett", "blackman", "blackman-harris", "bohman", "boxcar", "cosine", "exponential", "flat-top", "hamming", "hanning", "lanczos", "nuttall", "parzen", "rectangular", "taylor", "tukey", "kaiser", "gaussian", ] = "hamming", alpha: float = 0.5, beta: float = 14.0, sigma: float = 7.0, ) -> np.ndarray: """Apply windowing to the input data. Args: x: X data y: Y data method: Windowing function. Defaults to "hamming". alpha: Tukey window parameter. Defaults to 0.5. beta: Kaiser window parameter. Defaults to 14.0. sigma: Gaussian window parameter. Defaults to 7.0. Returns: Windowed Y data """ # Cases without parameters: win_func = { "barthann": scipy.signal.windows.barthann, "bartlett": np.bartlett, "blackman": np.blackman, "blackman-harris": scipy.signal.windows.blackmanharris, "bohman": scipy.signal.windows.bohman, "boxcar": scipy.signal.windows.boxcar, "cosine": scipy.signal.windows.cosine, "exponential": scipy.signal.windows.exponential, "flat-top": scipy.signal.windows.flattop, "hamming": np.hamming, "hanning": np.hanning, "lanczos": scipy.signal.windows.lanczos, "nuttall": scipy.signal.windows.nuttall, "parzen": scipy.signal.windows.parzen, "rectangular": np.ones, "taylor": scipy.signal.windows.taylor, }.get(method) if win_func is not None: return y * win_func(len(y)) # Cases with parameters: if method == "tukey": return y * scipy.signal.windows.tukey(len(y), alpha) if method == "kaiser": return y * np.kaiser(len(y), beta) if method == "gaussian": return y * scipy.signal.windows.gaussian(len(y), sigma) raise ValueError(f"Invalid window type {method}")
# MARK: Curve fitting models -----------------------------------------------------------
[docs] class FitModel(abc.ABC): """Curve fitting model base class"""
[docs] @classmethod @abc.abstractmethod def func(cls, x, amp, sigma, x0, y0): """Return fitting function"""
# pylint: disable=unused-argument
[docs] @classmethod def get_amp_from_amplitude(cls, amplitude, sigma): """Return amp from function amplitude and sigma""" return amplitude
[docs] @classmethod def amplitude(cls, amp, sigma): """Return function amplitude""" return cls.func(0, amp, sigma, 0, 0)
[docs] @classmethod @abc.abstractmethod def fwhm(cls, amp, sigma): """Return function FWHM"""
[docs] @classmethod def half_max_segment(cls, amp, sigma, x0, y0): """Return segment coordinates for y=half-maximum intersection""" hwhm = 0.5 * cls.fwhm(amp, sigma) yhm = 0.5 * cls.amplitude(amp, sigma) + y0 return x0 - hwhm, yhm, x0 + hwhm, yhm
[docs] class GaussianModel(FitModel): """1-dimensional Gaussian fit model"""
[docs] @classmethod def func(cls, x, amp, sigma, x0, y0): """Return fitting function""" return ( amp / (sigma * np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((x - x0) / sigma) ** 2) + y0 )
[docs] @classmethod def get_amp_from_amplitude(cls, amplitude, sigma): """Return amp from function amplitude and sigma""" return amplitude * (sigma * np.sqrt(2 * np.pi))
[docs] @classmethod def amplitude(cls, amp, sigma): """Return function amplitude""" return amp / (sigma * np.sqrt(2 * np.pi))
[docs] @classmethod def fwhm(cls, amp, sigma): """Return function FWHM""" return 2 * sigma * np.sqrt(2 * np.log(2))
[docs] class LorentzianModel(FitModel): """1-dimensional Lorentzian fit model"""
[docs] @classmethod def func(cls, x, amp, sigma, x0, y0): """Return fitting function""" return (amp / (sigma * np.pi)) / (1 + ((x - x0) / sigma) ** 2) + y0
[docs] @classmethod def get_amp_from_amplitude(cls, amplitude, sigma): """Return amp from function amplitude and sigma""" return amplitude * (sigma * np.pi)
[docs] @classmethod def amplitude(cls, amp, sigma): """Return function amplitude""" return amp / (sigma * np.pi)
[docs] @classmethod def fwhm(cls, amp, sigma): """Return function FWHM""" return 2 * sigma
[docs] class VoigtModel(FitModel): """1-dimensional Voigt fit model"""
[docs] @classmethod def func(cls, x, amp, sigma, x0, y0): """Return fitting function""" # pylint: disable=no-member z = (x - x0 + 1j * sigma) / (sigma * np.sqrt(2.0)) return y0 + amp * scipy.special.wofz(z).real / (sigma * np.sqrt(2 * np.pi))
[docs] @classmethod def fwhm(cls, amp, sigma): """Return function FWHM""" wg = GaussianModel.fwhm(amp, sigma) wl = LorentzianModel.fwhm(amp, sigma) return 0.5346 * wl + np.sqrt(0.2166 * wl**2 + wg**2)
# MARK: Misc. analyses -----------------------------------------------------------------
[docs] def find_nearest_zero_point_idx(y: np.ndarray) -> np.ndarray: """Find the x indices where the corresponding y is the closest to zero Args: y: Y data Returns: Indices of the points right before or at zero crossing """ xi = np.where((y[:-1] >= 0) & (y[1:] <= 0) | (y[:-1] <= 0) & (y[1:] >= 0))[0] return xi
[docs] def find_x_at_value(x: np.ndarray, y: np.ndarray, value: float) -> np.ndarray: """Find the x value where the y value is the closest to the given value using linear interpolation to deduce the precise x value. Args: x: X data y: Y data value: Value to find Returns: X value where the Y value is the closest to the given value """ leveled_y = y - value xi_before = find_nearest_zero_point_idx(leveled_y) xi_after = xi_before + 1 if len(xi_before) == 0: return np.array([0.0]) # linear interpolation p = (leveled_y[xi_after] - leveled_y[xi_before]) / (x[xi_after] - x[xi_before]) ori = leveled_y[xi_after] - p * x[xi_after] x0 = -ori / p # where the curve cut the absissa return x0
[docs] def bandwidth(data: np.ndarray, level: float = 3.0) -> float: """Compute the bandwidth of the signal at a given level. Args: data: X,Y data level: Level in dB at which the bandwidth is computed. Defaults to 3.0. Returns: Bandwidth of the signal at the given level: segment coordinates """ x, y = data half_max: float = np.max(y) - level bw = find_x_at_value(x, y, half_max)[0] coords = (x[0], half_max, bw, half_max) return coords
[docs] def contrast(y: np.ndarray) -> float: """Compute contrast Args: y: Input array Returns: Contrast """ max_, min_ = np.max(y), np.min(y) return (max_ - min_) / (max_ + min_)
# MARK: Dynamic parameters -------------------------------------------------------------
[docs] def sinusoidal_model( x: np.ndarray, a: float, f: float, phi: float, offset: float ) -> np.ndarray: """Sinusoidal model function.""" return a * np.sin(2 * np.pi * f * x + phi) + offset
[docs] def sinusoidal_fit( x: np.ndarray, y: np.ndarray ) -> tuple[tuple[float, float, float, float], float]: """Fit a sinusoidal model to the input data. Args: x: X data y: Y data Returns: A tuple containing the fit parameters (amplitude, frequency, phase, offset) and the residuals """ # Initial guess for the parameters # ================================================================================== offset = np.mean(y) amp = (np.max(y) - np.min(y)) / 2 phase_origin = 0 # Search for the maximum of the FFT i_maxfft = np.argmax(np.abs(np.fft.fft(y - offset))) if i_maxfft > len(x) / 2: # If the index is greater than N/2, we are in the mirrored half spectrum # (negative frequencies) i_maxfft = len(x) - i_maxfft freq = i_maxfft / (x[-1] - x[0]) # ================================================================================== def optfunc(fitparams: np.ndarray, x: np.ndarray, y: np.ndarray) -> np.ndarray: """Optimization function.""" return y - sinusoidal_model(x, *fitparams) # Fit the model to the data fitparams = scipy.optimize.leastsq( optfunc, [amp, freq, phase_origin, offset], args=(x, y) )[0] y_th = sinusoidal_model(x, *fitparams) residuals = np.std(y - y_th) return fitparams, residuals
[docs] def sinus_frequency(x: np.ndarray, y: np.ndarray) -> float: """Compute the frequency of a sinusoidal signal. Args: x: x signal data y: y signal data Returns: Frequency of the sinusoidal signal """ fitparams, _residuals = sinusoidal_fit(x, y) return fitparams[1]
[docs] def enob(x: np.ndarray, y: np.ndarray, full_scale: float = 1.0) -> float: """Compute Effective Number of Bits (ENOB). Args: x: x signal data y: y signal data full_scale: Full scale(V). Defaults to 1.0. Returns: Effective Number of Bits (ENOB) """ _fitparams, residuals = sinusoidal_fit(x, y) return -np.log2(residuals * np.sqrt(12) / full_scale)
[docs] def sinad( x: np.ndarray, y: np.ndarray, full_scale: float = 1.0, unit: Literal["dBc", "dBFS"] = "dBc", ) -> float: """Compute Signal-to-Noise and Distortion Ratio (SINAD). Args: x: x signal data y: y signal data full_scale: Full scale(V). Defaults to 1.0. unit: Unit of the input data. Valid values are 'dBc' and 'dBFS'. Defaults to 'dBc'. Returns: Signal-to-Noise and Distortion Ratio (SINAD) """ fitparams, residuals = sinusoidal_fit(x, y) amp = fitparams[0] # Compute the power of the fundamental powf = np.abs(amp / np.sqrt(2)) if unit == "dBc" else full_scale / (2 * np.sqrt(2)) return 20 * np.log10(powf / residuals)
[docs] def thd( x: np.ndarray, y: np.ndarray, full_scale: float = 1.0, unit: Literal["dBc", "dBFS"] = "dBc", nb_harm: int = 5, ) -> float: """Compute Total Harmonic Distortion (THD). Args: x: x signal data y: y signal data full_scale: Full scale(V). Defaults to 1.0. unit: Unit of the input data. Valid values are 'dBc' and 'dBFS'. Defaults to 'dBc'. nb_harm: Number of harmonics to consider. Defaults to 5. Returns: Total Harmonic Distortion (THD) """ fitparams, _residuals = sinusoidal_fit(x, y) offset = np.mean(y) amp, freq = fitparams[:2] ampfft = np.abs(np.fft.fft(y - offset)) # Compute the power of the fundamental if unit == "dBc": powfund = np.max(ampfft[: len(ampfft) // 2]) else: powfund = (full_scale / (2 * np.sqrt(2))) * (len(x) / np.sqrt(2)) sumharm = 0 for i in np.arange(nb_harm + 2)[2:]: a = i * np.ceil(freq * (x[-1] - x[0])) amp = ampfft[int(a - 5) : int(a + 5)] if len(amp) > 0: sumharm += np.max(amp) return 20 * np.log10(sumharm / powfund)
[docs] def sfdr( x: np.ndarray, y: np.ndarray, full_scale: float = 1.0, unit: Literal["dBc", "dBFS"] = "dBc", ) -> float: """Compute Spurious-Free Dynamic Range (SFDR). Args: x: x signal data y: y signal data full_scale: Full scale(V). Defaults to 1.0. unit: Unit of the input data. Valid values are 'dBc' and 'dBFS'. Defaults to 'dBc'. Returns: Spurious-Free Dynamic Range (SFDR) """ fitparams, _residuals = sinusoidal_fit(x, y) # Compute the power of the fundamental if unit == "dBc": powfund = np.max(np.abs(np.fft.fft(y))) else: powfund = (full_scale / (2 * np.sqrt(2))) * (len(x) / np.sqrt(2)) maxspike = np.max(np.abs(np.fft.fft(y - sinusoidal_model(x, *fitparams)))) return 20 * np.log10(powfund / maxspike)
[docs] def snr( x: np.ndarray, y: np.ndarray, full_scale: float = 1.0, unit: Literal["dBc", "dBFS"] = "dBc", ) -> float: """Compute Signal-to-Noise Ratio (SNR). Args: x: x signal data y: y signal data full_scale: Full scale(V). Defaults to 1.0. unit: Unit of the input data. Valid values are 'dBc' and 'dBFS'. Defaults to 'dBc'. Returns: Signal-to-Noise Ratio (SNR) """ fitparams, _residuals = sinusoidal_fit(x, y) # Compute the power of the fundamental if unit == "dBc": powfund = np.max(np.abs(np.fft.fft(y))) else: powfund = (full_scale / (2 * np.sqrt(2))) * (len(x) / np.sqrt(2)) noise = np.sqrt(np.mean((y - sinusoidal_model(x, *fitparams)) ** 2)) return 20 * np.log10(powfund / noise)
[docs] def sampling_period(x: np.ndarray) -> float: """Compute sampling period Args: x: X data Returns: Sampling period """ steps = np.diff(x) if not np.isclose(np.diff(steps).max(), 0, atol=1e-10): warnings.warn("Non-constant sampling signal") return steps[0]
[docs] def sampling_rate(x: np.ndarray) -> float: """Compute mean sampling rate Args: x: X data Returns: Sampling rate """ return 1.0 / sampling_period(x)
# MARK: Pulse analysis -----------------------------------------------------------------
[docs] def fwhm( data: np.ndarray, method: Literal["zero-crossing", "gauss", "lorentz", "voigt"] = "zero-crossing", xmin: float | None = None, xmax: float | None = None, ) -> tuple[float, float, float, float]: """Compute Full Width at Half Maximum (FWHM) of the input data Args: data: X,Y data method: Calculation method. Two types of methods are supported: a zero-crossing method and fitting methods (based on various models: Gauss, Lorentz, Voigt). Defaults to "zero-crossing". xmin: Lower X bound for the fitting. Defaults to None (no lower bound, i.e. the fitting starts from the first point). xmax: Upper X bound for the fitting. Defaults to None (no upper bound, i.e. the fitting ends at the last point) Returns: FWHM segment coordinates """ x, y = data dx, dy, base = np.max(x) - np.min(x), np.max(y) - np.min(y), np.min(y) sigma, mu = dx * 0.1, xpeak(x, y) if isinstance(xmin, float): indices = np.where(x >= xmin)[0] x = x[indices] y = y[indices] if isinstance(xmax, float): indices = np.where(x <= xmax)[0] x = x[indices] y = y[indices] if method == "zero-crossing": hmax = dy * 0.5 + np.min(y) fx = find_x_at_value(x, y, hmax) if fx.size != 2: warnings.warn(f"Ambiguous zero-crossing points (found {fx.size} points)") return fx[0], hmax, fx[-1], hmax try: FitModelClass: type[FitModel] = { "gauss": GaussianModel, "lorentz": LorentzianModel, "voigt": VoigtModel, }[method] except KeyError as exc: raise ValueError(f"Invalid method {method}") from exc def func(params): """Fitting model function""" # pylint: disable=cell-var-from-loop return y - FitModelClass.func(x, *params) amp = FitModelClass.get_amp_from_amplitude(dy, sigma) (amp, sigma, mu, base), _ier = scipy.optimize.leastsq( func, np.array([amp, sigma, mu, base]) ) return FitModelClass.half_max_segment(amp, sigma, mu, base)
[docs] def fw1e2(data: np.ndarray) -> tuple[float, float, float, float]: """Compute Full Width at 1/e� of the input data (using a Gaussian model fitting). Args: data: X,Y data Returns: FW at 1/e² segment coordinates """ x, y = data dx, dy, base = np.max(x) - np.min(x), np.max(y) - np.min(y), np.min(y) sigma, mu = dx * 0.1, xpeak(x, y) amp = GaussianModel.get_amp_from_amplitude(dy, sigma) p_in = np.array([amp, sigma, mu, base]) def func(params): """Fitting model function""" # pylint: disable=cell-var-from-loop return y - GaussianModel.func(x, *params) p_out, _ier = scipy.optimize.leastsq(func, p_in) amp, sigma, mu, base = p_out hw = 2 * sigma yhm = GaussianModel.amplitude(amp, sigma) / np.e**2 + base return mu - hw, yhm, mu + hw, yhm