From f60f22d5230bb935b6906e68feea4be07ca27b61 Mon Sep 17 00:00:00 2001 From: brunozc Date: Fri, 4 Jul 2025 10:01:22 +0200 Subject: [PATCH 1/3] fix: improve time padding logic in TimeSignalProcessing and add tests for time interpolation --- SignalProcessingTools/time_signal.py | 2 +- tests/test_time_signal.py | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/SignalProcessingTools/time_signal.py b/SignalProcessingTools/time_signal.py index fdd961c..d9283bd 100644 --- a/SignalProcessingTools/time_signal.py +++ b/SignalProcessingTools/time_signal.py @@ -110,7 +110,7 @@ def __init__(self, # pad signal at the end if necessary to get full windows if signal_length % window_size != 0: self.signal = np.append(self.signal, np.zeros(window_size - (signal_length % window_size))) - self.time = np.append(self.time, np.zeros(window_size - (signal_length % window_size))) + self.time = np.append(self.time, self.time[-1]+np.cumsum(np.ones(window_size - (signal_length % window_size))*(1/Fs))) self.operations.append(f"Signal padded with zeros (original length: {signal_length}, new length: {len(self.signal)})") def __str__(self) -> str: diff --git a/tests/test_time_signal.py b/tests/test_time_signal.py index 15b9cac..b9bbe90 100644 --- a/tests/test_time_signal.py +++ b/tests/test_time_signal.py @@ -4,6 +4,7 @@ from SignalProcessingTools.time_signal import TimeSignalProcessing, IntegrationRules, Windows TOL = 3e-3 +FLOAT_TOL = 1e-12 FREQ = 6 AMP = 1.75 @@ -215,6 +216,11 @@ def test_psd(test_data): np.testing.assert_almost_equal(sig.frequency_Pxx[np.argmax(sig.Pxx)], FREQ, 2) assert (np.abs((np.max(sig.Pxx) - peak_psd) / peak_psd) < 0.035) + # test the time interpolation + assert len(sig.time) == np.ceil(50001/4000)*4000 + assert ((np.diff(sig.time)-1/500) < FLOAT_TOL).all() + + def test_v_eff(): """ Test the v_eff function From d527ee41fb341f9f478b4bb7430b1c30e43c9ea2 Mon Sep 17 00:00:00 2001 From: brunozc Date: Fri, 4 Jul 2025 10:02:24 +0200 Subject: [PATCH 2/3] feat: add pre-commit configuration and update requirements for pre-commit package --- .pre-commit-config.yaml | 10 ++++++++++ requirements.txt | 3 ++- setup.cfg | 3 ++- 3 files changed, 14 insertions(+), 2 deletions(-) create mode 100644 .pre-commit-config.yaml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..bfd963e --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,10 @@ +repos: + - repo: https://github.com/google/yapf + rev: v0.43.0 + hooks: + - id: yapf + name: yapf + language: python + entry: yapf + args: [-i] + types: [python] diff --git a/requirements.txt b/requirements.txt index 906b331..7e435cc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,4 +2,5 @@ numpy==2.2.5 scipy==1.15.0 matplotlib==3.10.0 pytest==8.0.2 -tox==4.13.0 \ No newline at end of file +tox==4.13.0 +pre_commit==4.2.0 \ No newline at end of file diff --git a/setup.cfg b/setup.cfg index f6ab469..47b6a59 100644 --- a/setup.cfg +++ b/setup.cfg @@ -23,4 +23,5 @@ python_requires = >=3.10 [options.extras_require] testing = pytest>=8.0 - tox>=4.13 \ No newline at end of file + tox>=4.13 + pre_commit>=4.2.0 \ No newline at end of file From 771e469b9892bb43b0b399169349809397e8b15e Mon Sep 17 00:00:00 2001 From: brunozc Date: Fri, 4 Jul 2025 10:02:47 +0200 Subject: [PATCH 3/3] applied linter --- SignalProcessingTools/space_signal.py | 97 ++++++++---- SignalProcessingTools/time_signal.py | 204 ++++++++++++++++++-------- example_space_signal.py | 5 +- example_time_signal.py | 38 +++-- tests/test_space_signal.py | 35 +++-- tests/test_time_signal.py | 144 ++++++++++++------ 6 files changed, 349 insertions(+), 174 deletions(-) diff --git a/SignalProcessingTools/space_signal.py b/SignalProcessingTools/space_signal.py index 0d96cd2..138124e 100644 --- a/SignalProcessingTools/space_signal.py +++ b/SignalProcessingTools/space_signal.py @@ -8,7 +8,11 @@ class SpaceSignalProcessing: """ SignalProcessing class for processing signals in space. """ - def __init__(self, x: npt.NDArray[np.float64], values: npt.NDArray[np.float64], Fs: Optional[float] = None): + + def __init__(self, + x: npt.NDArray[np.float64], + values: npt.NDArray[np.float64], + Fs: Optional[float] = None): """ Initializes the ProcessSignal object. @@ -40,7 +44,6 @@ def __init__(self, x: npt.NDArray[np.float64], values: npt.NDArray[np.float64], self.max_fast = None self.max_fast_Dx = None - def compute_track_longitudinal_levels(self): """ Computes the track longitudinal levels, following EN 13848-1:2006. @@ -53,20 +56,34 @@ def compute_track_longitudinal_levels(self): - D3: 70m < lambda <= 150m (1/150 Hz < f <= 1/70 Hz) """ - sig = TimeSignalProcessing(self.coordinates, self.signal_raw, Fs=self.fs) - sig.filter([1/5., 1.], 4, type_filter="bandpass", filter_design=FilterDesign.BUTTERWORTH) + sig = TimeSignalProcessing(self.coordinates, + self.signal_raw, + Fs=self.fs) + sig.filter([1 / 5., 1.], + 4, + type_filter="bandpass", + filter_design=FilterDesign.BUTTERWORTH) self.d0 = sig.signal sig.reset() - sig.filter([1/25., 1/3.], 4, type_filter="bandpass", filter_design=FilterDesign.BUTTERWORTH) + sig.filter([1 / 25., 1 / 3.], + 4, + type_filter="bandpass", + filter_design=FilterDesign.BUTTERWORTH) self.d1 = sig.signal sig.reset() - sig.filter([1/70., 1/25.], 4, type_filter="bandpass", filter_design=FilterDesign.BUTTERWORTH) + sig.filter([1 / 70., 1 / 25.], + 4, + type_filter="bandpass", + filter_design=FilterDesign.BUTTERWORTH) self.d2 = sig.signal sig.reset() - sig.filter([1/150., 1/70.], 4, type_filter="bandpass", filter_design=FilterDesign.BUTTERWORTH) + sig.filter([1 / 150., 1 / 70.], + 4, + type_filter="bandpass", + filter_design=FilterDesign.BUTTERWORTH) self.d3 = sig.signal sig.reset() @@ -82,22 +99,23 @@ def compute_Hmax(self, convert_m2mm: bool = True): """ # octave bands used for the processing - one_third_octave_bands = [[.08, .10], - [.10, .126], - [.126, .16], - [.16, .20], - [.20, .253], - [.253, .32], - [.32, .40], - [.40, .50], - [.50, .63], - ] - + one_third_octave_bands = [ + [.08, .10], + [.10, .126], + [.126, .16], + [.16, .20], + [.20, .253], + [.253, .32], + [.32, .40], + [.40, .50], + [.50, .63], + ] # setting for the processing self.DXmaxFast = 1 nb_fft_min = 256 # minimum number of samples for the power spectral density - derivative = [0, 0, 0, 2, 2, 2, 2, 2, 2] # number of times that each frequency band is derived + derivative = [0, 0, 0, 2, 2, 2, 2, 2, + 2] # number of times that each frequency band is derived # RMS of the square root of the power spectral density self.rms_bands = np.zeros(len(one_third_octave_bands)) @@ -111,26 +129,34 @@ def compute_Hmax(self, convert_m2mm: bool = True): self.signal = self.signal_raw * 1000 # compute the power spectral density - n_fft = int(np.max([2 ** (np.ceil(np.log2(len(self.signal)))), nb_fft_min])) + n_fft = int( + np.max([2**(np.ceil(np.log2(len(self.signal)))), nb_fft_min])) # if signal is odd length, add a zero to make it even if len(self.signal) % 2 != 0: signal = np.append(self.signal, 0) - coordinates = np.append(self.coordinates, self.coordinates[-1] + (self.coordinates[1] - self.coordinates[0])) + coordinates = np.append( + self.coordinates, self.coordinates[-1] + + (self.coordinates[1] - self.coordinates[0])) else: signal = self.signal coordinates = self.coordinates - sig = TimeSignalProcessing(coordinates, signal, Fs=self.fs, window=Windows.HAMMING, + sig = TimeSignalProcessing(coordinates, + signal, + Fs=self.fs, + window=Windows.HAMMING, window_size=len(signal)) sig.psd(nb_points=n_fft, detrend=False) # compute the rsm psd - self.__rms_effective(sig.frequency_Pxx, sig.Pxx, one_third_octave_bands, derivative) + self.__rms_effective(sig.frequency_Pxx, sig.Pxx, + one_third_octave_bands, derivative) # compute the effective values self.__effective_values(one_third_octave_bands, derivative) - - def __rms_effective(self, frequency: npt.NDArray[np.float64], Pxx: npt.NDArray[np.float64], - one_third_octave_bands: List[Tuple[float, float]], derivative: List[int]): + def __rms_effective(self, frequency: npt.NDArray[np.float64], + Pxx: npt.NDArray[np.float64], + one_third_octave_bands: List[Tuple[float, float]], + derivative: List[int]): """ Computes RMS square root of power spectral density @@ -148,10 +174,13 @@ def __rms_effective(self, frequency: npt.NDArray[np.float64], Pxx: npt.NDArray[n for i, band in enumerate(one_third_octave_bands): # find indexes where the bands exist idx = np.where((frequency >= band[0]) & (frequency < band[1]))[0] - Pxx[idx] = (2 * np.pi * frequency[idx]) ** (2 * derivative[i]) * Pxx[idx] + Pxx[idx] = (2 * np.pi * frequency[idx])**(2 * + derivative[i]) * Pxx[idx] self.rms_bands[i] = np.sqrt(np.sum(Pxx[idx] * delta_f)) - def __effective_values(self, one_third_octave_bands: List[Tuple[float, float]], derivative: List[int]): + def __effective_values(self, one_third_octave_bands: List[Tuple[float, + float]], + derivative: List[int]): """ Computes the effective values of the signal @@ -170,8 +199,13 @@ def __effective_values(self, one_third_octave_bands: List[Tuple[float, float]], for i, band in enumerate(one_third_octave_bands): derivative_value = derivative[i] - sig = TimeSignalProcessing(self.coordinates, self.signal, Fs=self.fs) - sig.filter(np.array(band), N=3, type_filter="bandpass", filter_design=FilterDesign.BUTTERWORTH) + sig = TimeSignalProcessing(self.coordinates, + self.signal, + Fs=self.fs) + sig.filter(np.array(band), + N=3, + type_filter="bandpass", + filter_design=FilterDesign.BUTTERWORTH) new_signal = sig.signal while derivative_value != 0: @@ -181,7 +215,8 @@ def __effective_values(self, one_third_octave_bands: List[Tuple[float, float]], ksi = np.linspace(0, n * tau, int(n * tau / dx + 1)) g = fout * np.exp(-ksi / tau) - convoluted_signal = np.sqrt(np.convolve(new_signal**2, g) * dx / tau) + convoluted_signal = np.sqrt( + np.convolve(new_signal**2, g) * dx / tau) self.max_fast[i] = np.max(convoluted_signal) idx = np.floor((len(self.signal) - np.floor(self.DXmaxFast / dx)) / 2) + \ np.linspace(0, np.floor(self.DXmaxFast / dx)-1, int(np.floor(self.DXmaxFast / dx))) diff --git a/SignalProcessingTools/time_signal.py b/SignalProcessingTools/time_signal.py index d9283bd..52bedf3 100644 --- a/SignalProcessingTools/time_signal.py +++ b/SignalProcessingTools/time_signal.py @@ -14,6 +14,7 @@ class FilterDesign(Enum): CHEBYSHEV = 2 ELLIPTIC = 3 + class IntegrationRules(Enum): """ Integration rules @@ -21,6 +22,7 @@ class IntegrationRules(Enum): TRAPEZOID = 1 SIMPSON = 2 + class Windows(Enum): """ Windows types @@ -34,10 +36,12 @@ class Windows(Enum): RECTANGULAR = 'boxcar' TRIANG = 'triang' + class TimeSignalProcessing: """ Signal processing class for time signals """ + def __init__(self, time: npt.NDArray[np.float64], signal: npt.NDArray[np.float64], @@ -71,8 +75,7 @@ def __init__(self, self.Sxx = None self.frequency_Sxx = None self.time_Sxx = None - self.fft_settings = {"nb_points": None, - "half_representation": False} + self.fft_settings = {"nb_points": None, "half_representation": False} # Track operations performed on the signal self.operations = [] @@ -93,25 +96,38 @@ def __init__(self, self.use_window = False else: if window_size == 0: - raise ValueError("When using a window the `window_size` must be specified") + raise ValueError( + "When using a window the `window_size` must be specified") if window_size % 2 != 0: raise ValueError("Window length must be even") if window not in Windows: - raise ValueError(f"Window type {window} not supported. Available types: {list(Windows)}") + raise ValueError( + f"Window type {window} not supported. Available types: {list(Windows)}" + ) if window_size > signal_length: - raise ValueError(f"Window length ({window_size}) cannot be greater than signal length ({signal_length}).") + raise ValueError( + f"Window length ({window_size}) cannot be greater than signal length ({signal_length})." + ) self.window = self.__create_window(window, window_size) self.window_size = window_size self.window_type = window - self.nb_windows = int(np.ceil((signal_length / window_size) * 2 - 2)) + self.nb_windows = int( + np.ceil((signal_length / window_size) * 2 - 2)) self.use_window = True # pad signal at the end if necessary to get full windows if signal_length % window_size != 0: - self.signal = np.append(self.signal, np.zeros(window_size - (signal_length % window_size))) - self.time = np.append(self.time, self.time[-1]+np.cumsum(np.ones(window_size - (signal_length % window_size))*(1/Fs))) - self.operations.append(f"Signal padded with zeros (original length: {signal_length}, new length: {len(self.signal)})") + self.signal = np.append( + self.signal, + np.zeros(window_size - (signal_length % window_size))) + self.time = np.append( + self.time, self.time[-1] + np.cumsum( + np.ones(window_size - (signal_length % window_size)) * + (1 / Fs))) + self.operations.append( + f"Signal padded with zeros (original length: {signal_length}, new length: {len(self.signal)})" + ) def __str__(self) -> str: """ @@ -148,7 +164,8 @@ def __str__(self) -> str: return "\n".join(info) @staticmethod - def __create_window(window_type: Windows, size: int) -> npt.NDArray[np.float64]: + def __create_window(window_type: Windows, + size: int) -> npt.NDArray[np.float64]: """ Create a window array of specified type and size @@ -234,10 +251,10 @@ def fft(self, # peak amplitude of stationary sinusoids. spectrum_w[:, w] = np.fft.fft(signal_w, nfft) / normalise_fct - self.amplitude = np.mean(np.abs(spectrum_w), axis=1) # self.phase = np.unwrap(np.angle(np.mean(spectrum_w, axis=1))) - self.phase = np.angle(np.mean(np.exp(1j * np.angle(spectrum_w)), axis=1)) + self.phase = np.angle( + np.mean(np.exp(1j * np.angle(spectrum_w)), axis=1)) # compute frequency self.frequency = np.linspace(0, 1, nfft) * self.Fs @@ -249,9 +266,11 @@ def fft(self, self.phase = self.phase[:int(nfft / 2)] # FFT settings: needed to perform inverse FFT - self.fft_settings = {"nb_points": nfft, - "half_representation": half_representation, - "odd_length": odd_length} + self.fft_settings = { + "nb_points": nfft, + "half_representation": half_representation, + "odd_length": odd_length + } # Add to operations list op_info = f"FFT (points: {nfft}, half representation: {half_representation})" @@ -275,7 +294,8 @@ def inv_fft(self): "Please compute FFT with full representation.") if self.use_window: - raise ValueError("Cannot perform inverse FFT on the windowed signal.") + raise ValueError( + "Cannot perform inverse FFT on the windowed signal.") # get FFT settings odd_length = self.fft_settings["odd_length"] @@ -290,7 +310,8 @@ def inv_fft(self): # inverse of the FFT signal self.signal_inv = np.real(spectrum_inv) * len(spectrum) # time from frequency - self.time_inv = np.cumsum(np.ones(len(spectrum)) * 1 / self.Fs) - 1 / self.Fs + self.time_inv = np.cumsum( + np.ones(len(spectrum)) * 1 / self.Fs) - 1 / self.Fs if odd_length: # remove last sample @@ -298,11 +319,17 @@ def inv_fft(self): self.time_inv = self.time_inv[:-1] # Add to operations list - self.operations.append("Inverse FFT" + (" with windowing" if self.use_window else "")) - - def integrate(self, rule: IntegrationRules = IntegrationRules.TRAPEZOID, - baseline: bool = False, moving: bool = False, hp: bool = False, ini_cond: float = 0., - fpass: float = 0.5, n: int = 6): + self.operations.append("Inverse FFT" + + (" with windowing" if self.use_window else "")) + + def integrate(self, + rule: IntegrationRules = IntegrationRules.TRAPEZOID, + baseline: bool = False, + moving: bool = False, + hp: bool = False, + ini_cond: float = 0., + fpass: float = 0.5, + n: int = 6): """ Numerical integration of signal @@ -322,9 +349,13 @@ def integrate(self, rule: IntegrationRules = IntegrationRules.TRAPEZOID, # integration rule if rule == IntegrationRules.TRAPEZOID: - self.signal = integrate.cumulative_trapezoid(self.signal, self.time, initial=ini_cond) + self.signal = integrate.cumulative_trapezoid(self.signal, + self.time, + initial=ini_cond) elif rule == IntegrationRules.SIMPSON: - self.signal = integrate.cumulative_simpson(self.signal, x=self.time, initial=ini_cond) + self.signal = integrate.cumulative_simpson(self.signal, + x=self.time, + initial=ini_cond) else: sys.exit("Integration rule not supported") @@ -346,13 +377,18 @@ def integrate(self, rule: IntegrationRules = IntegrationRules.TRAPEZOID, if moving: op_details.append("moving average correction") if hp: - op_details.append(f"highpass filter (cutoff: {fpass} Hz, order: {n})") + op_details.append( + f"highpass filter (cutoff: {fpass} Hz, order: {n})") self.operations.append(f"Integration ({', '.join(op_details)})") - - def filter(self, Fpass: float, N: int, filter_design: FilterDesign = FilterDesign.ELLIPTIC, - type_filter: str = "lowpass", rp: float = 0.01, rs: int = 60): + def filter(self, + Fpass: float, + N: int, + filter_design: FilterDesign = FilterDesign.ELLIPTIC, + type_filter: str = "lowpass", + rp: float = 0.01, + rs: int = 60): """ Filter signal @@ -377,11 +413,23 @@ def filter(self, Fpass: float, N: int, filter_design: FilterDesign = FilterDesig # design filter if filter_design == FilterDesign.ELLIPTIC: - z, p, k = signal.ellip(N, rp, rs, np.array(Fpass) / (self.Fs / 2), btype=type_filter, output='zpk') + z, p, k = signal.ellip(N, + rp, + rs, + np.array(Fpass) / (self.Fs / 2), + btype=type_filter, + output='zpk') elif filter_design == FilterDesign.BUTTERWORTH: - z, p, k = signal.butter(N, np.array(Fpass) / (self.Fs / 2), btype=type_filter, output='zpk') + z, p, k = signal.butter(N, + np.array(Fpass) / (self.Fs / 2), + btype=type_filter, + output='zpk') elif filter_design == FilterDesign.CHEBYSHEV: - z, p, k = signal.cheby1(N, rp, np.array(Fpass) / (self.Fs / 2), btype=type_filter, output='zpk') + z, p, k = signal.cheby1(N, + rp, + np.array(Fpass) / (self.Fs / 2), + btype=type_filter, + output='zpk') sos = signal.zpk2sos(z, p, k) @@ -389,8 +437,8 @@ def filter(self, Fpass: float, N: int, filter_design: FilterDesign = FilterDesig self.signal = signal.sosfiltfilt(sos, self.signal) # Add to operations list - self.operations.append(f"Filter ({type_filter}, cutoff: {Fpass} Hz, order: {N})") - + self.operations.append( + f"Filter ({type_filter}, cutoff: {Fpass} Hz, order: {N})") def psd(self, detrend: str = "linear", nb_points: Optional[int] = None): """ @@ -403,11 +451,15 @@ def psd(self, detrend: str = "linear", nb_points: Optional[int] = None): """ if detrend not in ["linear", False]: - raise ValueError("Detrend method not supported. Available methods: ['linear', False]") + raise ValueError( + "Detrend method not supported. Available methods: ['linear', False]" + ) # check if window is initialized if not self.use_window: - raise ValueError("No window defined. Please define a window when initialising SignalProcessing.") + raise ValueError( + "No window defined. Please define a window when initialising SignalProcessing." + ) # if nb_points is None: nb_points is window length if nb_points is None: @@ -416,12 +468,18 @@ def psd(self, detrend: str = "linear", nb_points: Optional[int] = None): nfft = nb_points # compute PSD using Welch method - self.frequency_Pxx, self.Pxx = signal.welch(self.signal, fs=self.Fs, nperseg=self.window_size, nfft=nfft, - window=self.window_type.value, scaling='density', detrend=detrend) + self.frequency_Pxx, self.Pxx = signal.welch( + self.signal, + fs=self.Fs, + nperseg=self.window_size, + nfft=nfft, + window=self.window_type.value, + scaling='density', + detrend=detrend) # Add to operations list - self.operations.append(f"PSD (window: {self.window_type.name}, size: {self.window_size})") - + self.operations.append( + f"PSD (window: {self.window_type.name}, size: {self.window_size})") def v_eff_SBR(self, n: int = 4, tau: float = 0.125): """ @@ -438,13 +496,12 @@ def v_eff_SBR(self, n: int = 4, tau: float = 0.125): qsi = np.linspace(0, n * tau, int(n * tau * self.Fs + 1)) g = fout * np.exp(-qsi / tau) - # Frequency weighting parameters v0 = 1 / 1000 # Reference velocity [m/s] - f0 = 5.6 # Reference frequency [Hz] + f0 = 5.6 # Reference frequency [Hz] # Handle even/odd signal length for FFT - if self.signal.shape[0] % 2 != 0: + if self.signal.shape[0] % 2 != 0: nv1 = int(self.signal.shape[0] / 2 + 0.5) nv2 = int(self.signal.shape[0] / 2 - 0.5) else: @@ -456,25 +513,25 @@ def v_eff_SBR(self, n: int = 4, tau: float = 0.125): freq = np.arange(df, (nv1 + 1) * df, df) # Create high-pass weighting filter (human perception curve) - Hv = (1 / v0) * 1 / (np.sqrt(1 + (f0 / freq) ** 2)) + Hv = (1 / v0) * 1 / (np.sqrt(1 + (f0 / freq)**2)) Hv = np.append(0, Hv) # Add DC component # Create low-pass filter with 50 Hz cutoff cut_off_number = int(np.ceil(50 / df)) if cut_off_number < nv1: Hv2 = np.zeros(Hv.shape[0]) - Hv2[:cut_off_number+1] = 1 + Hv2[:cut_off_number + 1] = 1 else: Hv2 = np.ones(Hv.shape[0]) # Applies the frequency weighting functions Fv = np.fft.fft(self.signal) - Fhv = Hv2 * Hv * Fv[:nv1+1] + Fhv = Hv2 * Hv * Fv[:nv1 + 1] Fv = np.append(Fhv, np.flipud(np.conj(Fhv[1:nv2]))) v_eff = np.real(np.fft.ifft(Fv)) # moving root-mean-square through convolution with the exponential decay function `g` - v_eff = np.sqrt( np.convolve(v_eff**2, g) * (1 / self.Fs) /tau) + v_eff = np.sqrt(np.convolve(v_eff**2, g) * (1 / self.Fs) / tau) self.v_eff = v_eff[:self.signal.shape[0]] @@ -514,14 +571,19 @@ def spectrogram(self): Compute spectrogram of signal """ # compute spectrogram - f, t, Sxx = signal.spectrogram(self.signal, fs=self.Fs, window=self.window_type.value, - nperseg=self.window_size, noverlap=self.window_size // 8) + f, t, Sxx = signal.spectrogram(self.signal, + fs=self.Fs, + window=self.window_type.value, + nperseg=self.window_size, + noverlap=self.window_size // 8) self.Sxx = Sxx self.frequency_Sxx = f self.time_Sxx = t # Add to operations list - self.operations.append(f"Spectrogram (nperseg: {self.window_size}, noverlap: {self.window_size // 8})") + self.operations.append( + f"Spectrogram (nperseg: {self.window_size}, noverlap: {self.window_size // 8})" + ) def one_third_octave_bands(self): """ @@ -534,24 +596,31 @@ def one_third_octave_bands(self): # https://en.wikipedia.org/wiki/Octave_band initial_frequency_band_number = -20 final_frequency_band_number = 33 - names = ("10", "12.5", "16", "20", "25", "31.5", "40", "50", "63", "80", "100", "125", "160", "200", "250", "315", "400", - "500", "630", "800", "1000", "1250", "1600", "2000", "2500", "3.150", "4000", "5000", "6300", "8000", - "10000", "12500", "16000", "20000") + names = ("10", "12.5", "16", "20", "25", "31.5", "40", "50", "63", + "80", "100", "125", "160", "200", "250", "315", "400", "500", + "630", "800", "1000", "1250", "1600", "2000", "2500", "3.150", + "4000", "5000", "6300", "8000", "10000", "12500", "16000", + "20000") # compute centre frequencies of the bands - f_centre = 1000 * (2 ** (np.arange(initial_frequency_band_number, final_frequency_band_number) / 3)) - f_upper = f_centre * (2 ** (1 / 6)) - f_lower = f_centre / (2 ** (1 / 6)) + f_centre = 1000 * (2**(np.arange(initial_frequency_band_number, + final_frequency_band_number) / 3)) + f_upper = f_centre * (2**(1 / 6)) + f_lower = f_centre / (2**(1 / 6)) # sum the signal for the bands if (self.Pxx is None) and (self.amplitude is None): - raise ValueError("No PSD nor FFT computed. Please compute either first.") + raise ValueError( + "No PSD nor FFT computed. Please compute either first.") if self.Pxx is not None: # determine the frequency bands - idx = np.where((f_lower < np.max(self.frequency_Pxx)) & (f_upper > np.min(self.frequency_Pxx)))[0] + idx = np.where((f_lower < np.max(self.frequency_Pxx)) + & (f_upper > np.min(self.frequency_Pxx)))[0] if len(idx) == 0: - raise ValueError("No frequency bands found in the PSD. Please check the frequency bands.") + raise ValueError( + "No frequency bands found in the PSD. Please check the frequency bands." + ) delta_f = self.frequency_Pxx[1] - self.frequency_Pxx[0] @@ -561,14 +630,19 @@ def one_third_octave_bands(self): for i, val in enumerate(idx): self.octave_bands_Pxx[i] = float(names[val]) - mask = (self.frequency_Pxx >= f_lower[val]) & (self.frequency_Pxx < f_upper[val]) - self.octave_bands_Pxx_power[i] = np.sum(self.Pxx[mask] * delta_f) + mask = (self.frequency_Pxx + >= f_lower[val]) & (self.frequency_Pxx < f_upper[val]) + self.octave_bands_Pxx_power[i] = np.sum(self.Pxx[mask] * + delta_f) if self.amplitude is not None: # determine the frequency bands - idx = np.where((f_lower < np.max(self.frequency)) & (f_upper > np.min(self.frequency)))[0] + idx = np.where((f_lower < np.max(self.frequency)) + & (f_upper > np.min(self.frequency)))[0] if len(idx) == 0: - raise ValueError("No frequency bands found in the FFT. Please check the frequency bands.") + raise ValueError( + "No frequency bands found in the FFT. Please check the frequency bands." + ) # compute the FFT for the bands self.octave_bands_fft = np.zeros(len(idx)) @@ -576,5 +650,7 @@ def one_third_octave_bands(self): for i, val in enumerate(idx): self.octave_bands_fft[i] = float(names[val]) - mask = (self.frequency >= f_lower[val]) & (self.frequency < f_upper[val]) - self.octave_bands_fft_power[i] = np.sum(self.amplitude[mask]**2) + mask = (self.frequency >= f_lower[val]) & (self.frequency + < f_upper[val]) + self.octave_bands_fft_power[i] = np.sum( + self.amplitude[mask]**2) diff --git a/example_space_signal.py b/example_space_signal.py index 2aae46b..7bb469e 100644 --- a/example_space_signal.py +++ b/example_space_signal.py @@ -3,14 +3,12 @@ from SignalProcessingTools.space_signal import SpaceSignalProcessing import pickle - # Create test data x = np.linspace(0, 100, 50001) omega = 2 * np.pi * 6 y = 1.75 * np.sin(omega * x) y_noise = y + 0.01 * np.sin(120 * x) - # Create a SpaceSignalProcessing object and demonstrate basic functionality print("----------------------------------------------") print("EXAMPLE 1: Track longitudinal level processing") @@ -63,7 +61,6 @@ plt.tight_layout() plt.show() - # Example 2: Computing Hmax parameters print("------------------------------------") print("EXAMPLE 2: Computing Hmax parameters") @@ -75,7 +72,7 @@ 0.002 * np.sin(2 * np.pi * 0.1 * x_track) + # Long wavelength component 0.001 * np.sin(2 * np.pi * 0.2 * x_track) + # Medium wavelength component 0.0005 * np.sin(2 * np.pi * 0.4 * x_track) + # Short wavelength component - 0.0002 * np.random.randn(len(x_track)) # Random noise + 0.0002 * np.random.randn(len(x_track)) # Random noise ) sig_hmax = SpaceSignalProcessing(x_track, track_irregularity) diff --git a/example_time_signal.py b/example_time_signal.py index b0b9bbf..a707e90 100644 --- a/example_time_signal.py +++ b/example_time_signal.py @@ -2,14 +2,12 @@ import matplotlib.pylab as plt from SignalProcessingTools.time_signal import TimeSignalProcessing, IntegrationRules, Windows - # Create test data x = np.linspace(0, 100, 50001) omega = 2 * np.pi * 6 y = 1.75 * np.sin(omega * x) y_noise = y + 0.01 * np.sin(120 * x) - # Create a SignalProcessing object and demonstrate basic functionality print("------------------------------------") print("EXAMPLE 1: Basic FFT and integration") @@ -47,12 +45,14 @@ sig.reset() print(sig) - # Example 2: Working with windowed signals and PSD print("---------------------------------------") print("EXAMPLE 2: Windowed processing and PSD") print("---------------------------------------") -sig_window = TimeSignalProcessing(x, y_noise, window=Windows.HAMMING, window_size=4096) +sig_window = TimeSignalProcessing(x, + y_noise, + window=Windows.HAMMING, + window_size=4096) print(sig_window) # Calculate and plot PSD @@ -69,8 +69,10 @@ # Calculate and plot spectrogram sig_window.spectrogram() plt.figure(figsize=(10, 6)) -plt.pcolormesh(sig_window.time_Sxx, sig_window.frequency_Sxx, - 10 * np.log10(sig_window.Sxx), shading='gouraud') +plt.pcolormesh(sig_window.time_Sxx, + sig_window.frequency_Sxx, + 10 * np.log10(sig_window.Sxx), + shading='gouraud') plt.ylim(0, 20) plt.ylabel('Frequency [Hz]') plt.xlabel('Time [s]') @@ -80,7 +82,6 @@ print(sig_window) - # Example 3: Signal filtering print("---------------------------") print("EXAMPLE 3: Signal filtering") @@ -102,7 +103,9 @@ sig_filter.filter(10, 4, type_filter="lowpass") plt.figure(figsize=(10, 6)) -plt.plot(sig_filter.time[:500], sig_filter.signal[:500], label='Filtered signal') +plt.plot(sig_filter.time[:500], + sig_filter.signal[:500], + label='Filtered signal') plt.plot(x[:500], y[:500], '--', label='Original clean signal') plt.xlabel('Time [s]') plt.ylabel('Amplitude') @@ -113,7 +116,6 @@ print(sig_filter) - # Example 4: Inverse FFT print("----------------------") print("EXAMPLE 4: Inverse FFT") @@ -137,7 +139,10 @@ plt.figure(figsize=(10, 6)) plt.plot(sig_ifft.time[:500], sig_ifft.signal[:500], label='Original signal') -plt.plot(sig_ifft.time_inv[:500], sig_ifft.signal_inv[:500], '--', label='Reconstructed signal') +plt.plot(sig_ifft.time_inv[:500], + sig_ifft.signal_inv[:500], + '--', + label='Reconstructed signal') plt.xlabel('Time [s]') plt.ylabel('Amplitude') plt.title('Signal Reconstruction with Inverse FFT') @@ -146,11 +151,11 @@ plt.show() # Calculate RMSE between original and reconstructed signals -rmse = np.sqrt(np.sum((sig_ifft.signal[:500] - sig_ifft.signal_inv[:500]) ** 2) / 500) +rmse = np.sqrt( + np.sum((sig_ifft.signal[:500] - sig_ifft.signal_inv[:500])**2) / 500) print(f"RMSE between original and reconstructed signals: {rmse:.6f}") print(sig_ifft) - # Example 5: Effective velocity calculation using SBR method print("-------------------------------------------") print("EXAMPLE 5: Effective velocity (SBR method)") @@ -158,7 +163,8 @@ # Create a vibration signal (using more complex multi-frequency signal) t = np.linspace(0, 10, 5001) -vib_signal = 0.5 * np.sin(2 * np.pi * 2 * t) + 0.3 * np.sin(2 * np.pi * 8 * t) + 0.2 * np.sin(2 * np.pi * 15 * t) +vib_signal = 0.5 * np.sin(2 * np.pi * 2 * t) + 0.3 * np.sin( + 2 * np.pi * 8 * t) + 0.2 * np.sin(2 * np.pi * 15 * t) sig_veff = TimeSignalProcessing(t, vib_signal) # Calculate effective velocity @@ -176,7 +182,6 @@ print(sig_veff) - # Example 6: Different integration rules print("---------------------------------------") print("EXAMPLE 6: Integration rules comparison") @@ -196,7 +201,8 @@ # Analytical integration for comparison (integral of sin(2πt) is -cos(2πt)/(2π)) vel_analytical = -np.cos(2 * np.pi * 1 * t_acc) / (2 * np.pi) -vel_analytical = vel_analytical - np.mean(vel_analytical) # baseline correction for comparison +vel_analytical = vel_analytical - np.mean( + vel_analytical) # baseline correction for comparison plt.figure(figsize=(10, 6)) plt.plot(t_acc, vel_analytical, 'k-', label='Analytical solution') @@ -212,4 +218,4 @@ print("Trapezoid integration:") print(sig_trap) print("\nSimpson integration:") -print(sig_simp) \ No newline at end of file +print(sig_simp) diff --git a/tests/test_space_signal.py b/tests/test_space_signal.py index 732f430..5d8b8fc 100644 --- a/tests/test_space_signal.py +++ b/tests/test_space_signal.py @@ -32,12 +32,16 @@ def test_track_quality_index(test_data): with open("./tests/data/track_alignment_results.txt", "r") as fi: data = json.load(fi) - assert np.allclose(sig.coordinates, data["coordinates"], rtol=1e-5, atol=1e-8) + assert np.allclose(sig.coordinates, + data["coordinates"], + rtol=1e-5, + atol=1e-8) assert np.allclose(sig.d0, data["D0"], rtol=1e-5, atol=1e-8) assert np.allclose(sig.d1, data["D1"], rtol=1e-5, atol=1e-8) assert np.allclose(sig.d2, data["D2"], rtol=1e-5, atol=1e-8) assert np.allclose(sig.d3, data["D3"], rtol=1e-5, atol=1e-8) + def test_Hmax(test_data): """ Test the Hmax function @@ -48,19 +52,22 @@ def test_Hmax(test_data): sig.compute_Hmax() - - rms_band_matlab = np.array([2087.55705457531, 1139.29553877343, 793.047457091564, 548.095561097181, - 648.015015656438, 521.608760233929, 790.563948013129, 886.097705321285, - 1342.44544888507]) - h_max_matlab = np.array([6557.20346546820, 3764.49040821894, 2505.54201229720, 1727.21064286318, - 1208.73808675389, 1182.77849807824, 1996.50964604940, 2682.54554936051, - 3681.41028314498]) - h_max_dx_matlab = np.array([293.263231128877, 783.641358934000, 1300.63422139907, 524.680132894043, - 139.978582308885, 583.726452878631, 662.410999843909, 1240.10908108192, - 1127.56942541745]) - - + rms_band_matlab = np.array([ + 2087.55705457531, 1139.29553877343, 793.047457091564, 548.095561097181, + 648.015015656438, 521.608760233929, 790.563948013129, 886.097705321285, + 1342.44544888507 + ]) + h_max_matlab = np.array([ + 6557.20346546820, 3764.49040821894, 2505.54201229720, 1727.21064286318, + 1208.73808675389, 1182.77849807824, 1996.50964604940, 2682.54554936051, + 3681.41028314498 + ]) + h_max_dx_matlab = np.array([ + 293.263231128877, 783.641358934000, 1300.63422139907, 524.680132894043, + 139.978582308885, 583.726452878631, 662.410999843909, 1240.10908108192, + 1127.56942541745 + ]) assert np.allclose(sig.rms_bands, rms_band_matlab, rtol=1e-3, atol=1e-8) assert np.allclose(sig.max_fast, h_max_matlab, rtol=1e-3, atol=1e-8) - assert np.allclose(sig.max_fast_Dx, h_max_dx_matlab, rtol=1e-3, atol=1e-8) \ No newline at end of file + assert np.allclose(sig.max_fast_Dx, h_max_dx_matlab, rtol=1e-3, atol=1e-8) diff --git a/tests/test_time_signal.py b/tests/test_time_signal.py index b9bbe90..60a649b 100644 --- a/tests/test_time_signal.py +++ b/tests/test_time_signal.py @@ -8,6 +8,7 @@ FREQ = 6 AMP = 1.75 + @pytest.fixture def test_data(): """ @@ -19,6 +20,7 @@ def test_data(): y_noise = y + 0.01 * np.sin(120 * x) return x, y, y_noise + def test_fft(test_data): """ Test the fft function @@ -33,10 +35,10 @@ def test_fft(test_data): assert len(sig.signal) == 50001 assert len(sig.time) == 50001 - np.testing.assert_almost_equal(sig.frequency[np.argmax(sig.amplitude)], FREQ, 2) + np.testing.assert_almost_equal(sig.frequency[np.argmax(sig.amplitude)], + FREQ, 2) np.testing.assert_almost_equal(np.max(sig.amplitude), AMP, 2) - # results full representation sig.fft(half_representation=False) @@ -44,16 +46,23 @@ def test_fft(test_data): assert len(sig.signal) == 50001 assert len(sig.time) == 50001 - np.testing.assert_almost_equal(sig.frequency[np.argmax(sig.amplitude[:int(len(sig.amplitude) / 2)])], FREQ, 2) - np.testing.assert_almost_equal(np.max(sig.amplitude[:int(len(sig.amplitude) / 2)]), AMP / 2, 2) + np.testing.assert_almost_equal( + sig.frequency[np.argmax(sig.amplitude[:int(len(sig.amplitude) / 2)])], + FREQ, 2) + np.testing.assert_almost_equal( + np.max(sig.amplitude[:int(len(sig.amplitude) / 2)]), AMP / 2, 2) # example with spectral leakage y = 1.75 * np.sin(2.675 * 2 * np.pi * x) sig = TimeSignalProcessing(x, y) sig.fft() - np.testing.assert_almost_equal(sig.frequency[np.argmax(sig.amplitude[:int(len(sig.amplitude) / 2)])], 2.675, 2) - np.testing.assert_almost_equal(np.max(sig.amplitude[:int(len(sig.amplitude) / 2)]), 1.137, 2) + np.testing.assert_almost_equal( + sig.frequency[np.argmax(sig.amplitude[:int(len(sig.amplitude) / 2)])], + 2.675, 2) + np.testing.assert_almost_equal( + np.max(sig.amplitude[:int(len(sig.amplitude) / 2)]), 1.137, 2) + def test_fft_nb_points(test_data): """ @@ -66,10 +75,11 @@ def test_fft_nb_points(test_data): sig.fft(nb_points=2**18) # check if signal lenght has been adapted to window size - assert len(sig.amplitude) == (2**18)/2 - assert len(sig.frequency) == (2**18)/2 + assert len(sig.amplitude) == (2**18) / 2 + assert len(sig.frequency) == (2**18) / 2 - np.testing.assert_almost_equal(sig.frequency[np.argmax(sig.amplitude)], FREQ, 2) + np.testing.assert_almost_equal(sig.frequency[np.argmax(sig.amplitude)], + FREQ, 2) np.testing.assert_almost_equal(np.max(sig.amplitude), AMP, 2) # results full representation @@ -79,8 +89,11 @@ def test_fft_nb_points(test_data): assert len(sig.amplitude) == 2**18 assert len(sig.frequency) == 2**18 - np.testing.assert_almost_equal(sig.frequency[np.argmax(sig.amplitude[:int(len(sig.amplitude) / 2)])], FREQ, 2) - np.testing.assert_almost_equal(np.max(sig.amplitude[:int(len(sig.amplitude) / 2)]), AMP / 2, 2) + np.testing.assert_almost_equal( + sig.frequency[np.argmax(sig.amplitude[:int(len(sig.amplitude) / 2)])], + FREQ, 2) + np.testing.assert_almost_equal( + np.max(sig.amplitude[:int(len(sig.amplitude) / 2)]), AMP / 2, 2) def test_fft_window(test_data): @@ -90,7 +103,9 @@ def test_fft_window(test_data): x, y, _ = test_data # assert that sig raises a Value error - with pytest.raises(ValueError, match="When using a window the `window_size` must be specified"): + with pytest.raises( + ValueError, + match="When using a window the `window_size` must be specified"): sig = TimeSignalProcessing(x, y, window=Windows.HAMMING) # test with window - half representation @@ -101,8 +116,11 @@ def test_fft_window(test_data): assert len(sig.signal) == 54000 assert len(sig.time) == 54000 - np.testing.assert_almost_equal(sig.frequency[np.argmax(sig.amplitude[:int(len(sig.amplitude))])], FREQ, 2) - np.testing.assert_almost_equal(np.max(sig.amplitude[:int(len(sig.amplitude))]), AMP, 2) + np.testing.assert_almost_equal( + sig.frequency[np.argmax(sig.amplitude[:int(len(sig.amplitude))])], + FREQ, 2) + np.testing.assert_almost_equal( + np.max(sig.amplitude[:int(len(sig.amplitude))]), AMP, 2) # full representation sig.fft(half_representation=False) @@ -111,8 +129,11 @@ def test_fft_window(test_data): assert len(sig.signal) == 54000 assert len(sig.time) == 54000 - np.testing.assert_almost_equal(sig.frequency[np.argmax(sig.amplitude[:int(len(sig.amplitude))])], FREQ, 2) - np.testing.assert_almost_equal(np.max(sig.amplitude[:int(len(sig.amplitude))]), AMP / 2, 2) + np.testing.assert_almost_equal( + sig.frequency[np.argmax(sig.amplitude[:int(len(sig.amplitude))])], + FREQ, 2) + np.testing.assert_almost_equal( + np.max(sig.amplitude[:int(len(sig.amplitude))]), AMP / 2, 2) def test_ifft(test_data): @@ -126,7 +147,11 @@ def test_ifft(test_data): sig.fft(half_representation=True) # assert that sig raises a Value error - with pytest.raises(NotImplementedError, match="Half representation not supported for inverse FFT. Please compute FFT with full representation."): + with pytest.raises( + NotImplementedError, + match= + "Half representation not supported for inverse FFT. Please compute FFT with full representation." + ): sig.inv_fft() # test with window - full representation @@ -138,8 +163,9 @@ def test_ifft(test_data): # check if signal lenght has been adapted to window size assert len(sig.signal) == len(sig.signal_inv) - rmse = np.sqrt(np.sum((sig.signal - sig.signal_inv) ** 2) / len(y)) - assert(rmse < TOL) + rmse = np.sqrt(np.sum((sig.signal - sig.signal_inv)**2) / len(y)) + assert (rmse < TOL) + def test_ifft_window(test_data): """ @@ -152,7 +178,11 @@ def test_ifft_window(test_data): sig.fft(half_representation=True) # assert that sig raises a Value error - with pytest.raises(NotImplementedError, match="Half representation not supported for inverse FFT. Please compute FFT with full representation."): + with pytest.raises( + NotImplementedError, + match= + "Half representation not supported for inverse FFT. Please compute FFT with full representation." + ): sig.inv_fft() # test with window - half representation @@ -160,7 +190,9 @@ def test_ifft_window(test_data): sig.fft(half_representation=False) # assert that sig raises a Value error - with pytest.raises(ValueError, match="Cannot perform inverse FFT on the windowed signal."): + with pytest.raises( + ValueError, + match="Cannot perform inverse FFT on the windowed signal."): sig.inv_fft() @@ -175,13 +207,18 @@ def test_int(test_data): omega = 2 * np.pi * FREQ int_sig = -AMP * np.cos(omega * x) / omega - rmse = np.sqrt(np.sum((sig.signal - int_sig) ** 2) / len(int_sig)) - assert(rmse < TOL) + rmse = np.sqrt(np.sum((sig.signal - int_sig)**2) / len(int_sig)) + assert (rmse < TOL) sig = TimeSignalProcessing(x, y) - sig.integrate(baseline=True, hp=True, rule=IntegrationRules.SIMPSON, fpass=1, n=6) - rmse = np.sqrt(np.sum((sig.signal - int_sig) ** 2) / len(int_sig)) - assert(rmse < TOL) + sig.integrate(baseline=True, + hp=True, + rule=IntegrationRules.SIMPSON, + fpass=1, + n=6) + rmse = np.sqrt(np.sum((sig.signal - int_sig)**2) / len(int_sig)) + assert (rmse < TOL) + def test_filter(test_data): """ @@ -192,15 +229,21 @@ def test_filter(test_data): sig.filter(10, 4, type_filter="lowpass") # compare between 200 and -200 to avoid edge effects - rmse = np.sqrt(np.sum((sig.signal[200:-200] - y[200:-200]) ** 2) / len(y[200:-200])) - assert(rmse < TOL) + rmse = np.sqrt( + np.sum((sig.signal[200:-200] - y[200:-200])**2) / len(y[200:-200])) + assert (rmse < TOL) + def test_psd(test_data): """ Test the psd function """ x, y, _ = test_data - with pytest.raises(ValueError, match="No window defined. Please define a window when initialising SignalProcessing."): + with pytest.raises( + ValueError, + match= + "No window defined. Please define a window when initialising SignalProcessing." + ): sig = TimeSignalProcessing(x, y) sig.psd() @@ -208,17 +251,18 @@ def test_psd(test_data): sig.psd() # power - power_sinus_wave = AMP ** 2 / 2 + power_sinus_wave = AMP**2 / 2 bin_width = sig.Fs / sig.window_size - ENBW = np.sum(sig.window ** 2) / (np.sum(sig.window)**2) * sig.window_size + ENBW = np.sum(sig.window**2) / (np.sum(sig.window)**2) * sig.window_size peak_psd = power_sinus_wave / (ENBW * bin_width) - np.testing.assert_almost_equal(sig.frequency_Pxx[np.argmax(sig.Pxx)], FREQ, 2) + np.testing.assert_almost_equal(sig.frequency_Pxx[np.argmax(sig.Pxx)], FREQ, + 2) assert (np.abs((np.max(sig.Pxx) - peak_psd) / peak_psd) < 0.035) # test the time interpolation - assert len(sig.time) == np.ceil(50001/4000)*4000 - assert ((np.diff(sig.time)-1/500) < FLOAT_TOL).all() + assert len(sig.time) == np.ceil(50001 / 4000) * 4000 + assert ((np.diff(sig.time) - 1 / 500) < FLOAT_TOL).all() def test_v_eff(): @@ -235,7 +279,7 @@ def test_v_eff(): v_eff = fi.read().splitlines() v_eff = np.array([list(map(float, r.split(";"))) for r in v_eff]) - time = np.linspace(0, (raw.shape[0]-1) / 500, raw.shape[0]) + time = np.linspace(0, (raw.shape[0] - 1) / 500, raw.shape[0]) # compute veff for i in range(raw.shape[1]): @@ -243,6 +287,7 @@ def test_v_eff(): sig.v_eff_SBR() np.testing.assert_almost_equal(sig.v_eff, np.array(v_eff)[:, i], 2) + def test_str_representation(test_data): """ Test the __str__ method to verify operations are tracked correctly @@ -283,6 +328,7 @@ def test_str_representation(test_data): assert any("Filter" in op for op in sig.operations) assert any("PSD" in op for op in sig.operations) + def test_spectrogram(test_data): """ Test the spectrogram function @@ -293,10 +339,11 @@ def test_spectrogram(test_data): # check if signal lenght has been adapted to window size assert sig.Sxx.shape == (301, 95) - assert sig.time_Sxx.shape == (95,) - assert sig.frequency_Sxx.shape == (301,) + assert sig.time_Sxx.shape == (95, ) + assert sig.frequency_Sxx.shape == (301, ) - np.testing.assert_almost_equal(sig.frequency_Sxx[np.where(sig.Sxx==np.max(sig.Sxx))[0][0]], FREQ, 0) + np.testing.assert_almost_equal( + sig.frequency_Sxx[np.where(sig.Sxx == np.max(sig.Sxx))[0][0]], FREQ, 0) np.testing.assert_almost_equal(np.max(sig.Sxx), 1.27, 2) # # plot spectrogram @@ -355,7 +402,10 @@ def test_reset(test_data): assert sig.time_Sxx is None # Verify FFT settings are reset - assert sig.fft_settings == {"nb_points": None, "half_representation": False} + assert sig.fft_settings == { + "nb_points": None, + "half_representation": False + } # Verify operations history is cleared assert len(sig.operations) == 0 @@ -363,6 +413,7 @@ def test_reset(test_data): # Verify string representation shows no operations assert "No operations performed yet" in str(sig) + def test_one_third_octave(test_data): """ Test the one third octave function @@ -370,11 +421,13 @@ def test_one_third_octave(test_data): x, y, _ = test_data # definition of frequencies used in the test - freqs_used = [10, 12.5, 16, 20, 25, 31.5, 40, 50, 63, 80, 100, 125, 160, 200, 250] + freqs_used = [ + 10, 12.5, 16, 20, 25, 31.5, 40, 50, 63, 80, 100, 125, 160, 200, 250 + ] # compute the power at 20 Hz - f_center = 1000 * (2 ** ((-17) / 3)) - f_max = f_center * (2 ** (1 / 6)) - f_min = f_center / (2 ** (1 / 6)) + f_center = 1000 * (2**((-17) / 3)) + f_max = f_center * (2**(1 / 6)) + f_min = f_center / (2**(1 / 6)) # test FFT sig = TimeSignalProcessing(x, y) @@ -385,7 +438,7 @@ def test_one_third_octave(test_data): idx = np.where((sig.frequency >= f_min) & (sig.frequency < f_max))[0] - assert sig.octave_bands_fft_power[3] == np.sum(sig.amplitude[idx] ** 2) + assert sig.octave_bands_fft_power[3] == np.sum(sig.amplitude[idx]**2) assert sig.octave_bands_fft[3] == 20 # test PSDF @@ -394,7 +447,8 @@ def test_one_third_octave(test_data): sig.one_third_octave_bands() assert all(sig.octave_bands_Pxx == freqs_used) - idx = np.where((sig.frequency_Pxx >= f_min) & (sig.frequency_Pxx < f_max))[0] + idx = np.where((sig.frequency_Pxx >= f_min) + & (sig.frequency_Pxx < f_max))[0] delta_freq = sig.frequency_Pxx[1] - sig.frequency_Pxx[0] assert sig.octave_bands_Pxx_power[3] == np.sum(sig.Pxx[idx] * delta_freq) assert sig.octave_bands_Pxx[3] == 20