diff --git a/pyproject.toml b/pyproject.toml index ee3605a..eea42b9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -36,6 +36,7 @@ show-report = true [tool.ruff] line-length = 72 +builtins = ["_"] [tool.basedpyright] venvPath = ".env" diff --git a/src/science_signal/__init__.py b/src/science_signal/__init__.py index 0d01e6d..22a0a4f 100644 --- a/src/science_signal/__init__.py +++ b/src/science_signal/__init__.py @@ -25,6 +25,7 @@ def interpolate(*signals: "Signal") -> tuple["Signal", ...]: Return each signal with a common frequency/time range (smallest range and highest rate) """ + logger.debug(_("interpolate {n} signals").format(n=len(signals))) splines = [signal.spline() for signal in signals] x_range = interpolate_abciss(*signals) diff --git a/src/science_signal/generator.py b/src/science_signal/generator.py index ac4bc64..560ca64 100644 --- a/src/science_signal/generator.py +++ b/src/science_signal/generator.py @@ -1,4 +1,4 @@ -from numpy import arange, pi, sin as np_sin +from numpy import arange, pi, sign, sin as np_sin from numpy.random import default_rng from science_signal.signal import Signal @@ -8,6 +8,7 @@ def sin( rate: float = 10, frequency: float = 1, amplitude: float = 1, + offset: float = 0, phase: float = 0, ) -> Signal: """ @@ -16,10 +17,11 @@ def sin( duration: duration of the signal in second rate: rate of the signal in Hz (s-1) frequency: frequency of the wanted sinus + offset: mean of the wanted sinus phase: between 0 and 2 pi (if more, congruate) """ x = arange(0, duration, 1 / rate) - y = amplitude * np_sin(2 * pi * frequency * x + phase) + y = offset + amplitude * np_sin(2 * pi * frequency * x + phase) return Signal( x, @@ -32,6 +34,7 @@ def cos( rate: float = 10, frequency: float = 1, amplitude: float = 1, + offset: float = 0, phase: float = 0, ) -> Signal: """ @@ -42,7 +45,36 @@ def cos( frequency: frequency of the wanted sinus phase: between 0 and 2 pi (if more, congruate) """ - return sin(duration, rate, frequency, amplitude, phase + pi / 2) + return sin( + duration, rate, frequency, amplitude, offset, phase + pi / 2 + ) + + +def square( + duration: float = 10, + rate: float = 10, + frequency: float = 1, + amplitude: float = 1, + offset: float = 0, + phase: float = 0, +) -> Signal: + """ + Generate square signal + + duration: duration of the signal in second + rate: rate of the signal in Hz (s-1) + frequency: frequency of the wanted square signal + phase: between 0 and 2 pi (if more, congruate) + """ + x = arange(0, duration, 1 / rate) + y = offset + amplitude * sign( + np_sin(2 * pi * frequency * x + phase) + ) + + return Signal( + x, + y, + ) def gaussian_noise( @@ -56,12 +88,36 @@ def gaussian_noise( Generate a gaussian noise signal duration: duration of the signal in second - rate: rate of the signal in hz (s-1) + rate: rate of the signal in Hz (s-1) sigma: standard deviation of the wanted distribution mu: mean of the distribution + seed: optional, allowd to specify a seed for testing purpose """ x = arange(0, duration, 1 / rate) return Signal( x, default_rng(seed).normal(mu, sigma, len(x)), ) + + +def white_noise( + duration: float = 10, + rate: float = 10, + minimum: float = -1, + maximum: float = 1, + seed: None | int = None, +) -> Signal: + """ + Generate a white noise signal + + duration: duration of the signal in seconds + rate: rate of the signal in hz (s-1) + minimum: minimum value of the noise + maximum: maximum value of the noise + seed: optional, allowd to specify a seed for testing purpose + """ + x = arange(0, duration, 1 / rate) + return Signal( + x, + default_rng(seed).uniform(minimum, maximum, len(x)), + ) diff --git a/src/science_signal/signal.py b/src/science_signal/signal.py index 721c586..169bcec 100644 --- a/src/science_signal/signal.py +++ b/src/science_signal/signal.py @@ -1,4 +1,4 @@ -from scipy.interpolate import CubicSpline # pyright: ignore[reportMissingTypeStubs] +from scipy.interpolate import CubicSpline # type: ignore[reportMissingTypeStubs] from numpy.typing import ArrayLike, NDArray from numpy import ( append, @@ -7,14 +7,13 @@ from numpy import ( linspace, array, logical_and, - logical_or, sin, where, arange, zeros, ) -from scipy.signal import detrend, welch -from scipy.fft import rfftfreq, rfft, irfft +from scipy.signal import detrend, welch # type: ignore[reportMissingTypeStubs] +from scipy.fft import rfftfreq, rfft, irfft # type: ignore[reportUnknownVariableType] from science_signal import interpolate @@ -77,7 +76,7 @@ class Signal: start = min(self.x) if end is None: end = max(self.x) - indexes = where(logical_and(self.x >= start, self.x <= end)) + indexes = where(logical_and(self.x >= start, self.x <= end)) # type: ignore[reportOperatorIssue] return Signal( self.x[indexes], self.y[indexes], @@ -114,30 +113,49 @@ class Signal: def filter( self, start: None | float = None, end: None | float = None ) -> "Signal": - freq_x = rfftfreq(len(self), self.sampling) + freq_x = rfftfreq(len(self), self.sampling) # type: ignore[reportUnknownVariableType] freq_y = rfft(self.y) + rate = 1 / (freq_x[1] - freq_x[0]) # type: ignore[reportUnknownVariableType] + if start is None: start = min(abs(freq_x)) if end is None: end = max(abs(freq_x)) - index_to_remove = where( - logical_or(abs(freq_x) <= start, abs(freq_x) >= end) - ) - freq_y[index_to_remove] = 0 + index_start = where(abs(freq_x) < start) # type: ignore[reportOperatorIssue] + index_end = where(abs(freq_x) > end) # type: ignore[reportOperatorIssue] + + if len(index_start) != 0: + """ + damping_start = 10 ** -( # type: ignore[reportUnknownVariableType] + arange(1, len(index_start[0]) + 1, 1, dtype=float64) + / rate + ) + freq_y[index_start] = freq_y[index_start] * damping_start # type: ignore[reportIndexIssue] + """ + freq_y[index_start] = 0 + if len(index_end) != 0: + """ + damping_end = 10 ** -( # type: ignore[reportUnknownVariableType] + arange(1, len(index_end[0]) + 1, 1, dtype=float64) + / rate + ) + freq_y[index_end] = freq_y[index_end] * damping_end # type: ignore[reportAny] + """ + freq_y[index_end] = 0 y = irfft(freq_y) return Signal( self.x[: len(y)], - y[: len(self.x)], + y[: len(self.x)], # type: ignore[reportArgumentType] ) def psd(self, fft_length: int = 10) -> "Signal": """ Compute psd of a given signal """ - freq, psd = welch( + freq, psd = welch( # type: ignore[reportUnknownVariableType] self.y, self.rate, nperseg=fft_length * self.rate, @@ -192,7 +210,7 @@ class Signal: # suppose same range but different rates return self.operator_signal( Signal( - linspace(self.x[0], self.x[-1], len(others)), + linspace(self.x[0], self.x[-1], len(other)), array(other), ), operator,