From dea64a7b7f6e5e31c2e8e0c66b49acb45cd07b51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Boulogne?= Date: Fri, 29 Aug 2025 18:04:30 +0200 Subject: [PATCH] upd --- mysignal/phasefreq.py | 53 ++++++++++++++++++++++++++++--------------- 1 file changed, 35 insertions(+), 18 deletions(-) diff --git a/mysignal/phasefreq.py b/mysignal/phasefreq.py index f07c4d4..eb6bf88 100644 --- a/mysignal/phasefreq.py +++ b/mysignal/phasefreq.py @@ -4,7 +4,7 @@ import pandas as pd import pywt from scipy.fft import fft, fftfreq -from scipy.signal import find_peaks, hilbert, butter, filtfilt, correlate +from scipy.signal import find_peaks, hilbert, butter, filtfilt, correlate, sosfiltfilt from scipy.stats import linregress from scipy.optimize import curve_fit @@ -63,26 +63,43 @@ def filter_signal_by_modes(time, signal, num_modes=1, bandwidth_factor=0.1, nyqu # Pour chaque fréquence dominante, appliquer un filtre passe-bande for freq in dominant_freqs: - # Définir la bande passante (10% de la fréquence centrale de chaque côté) + # Définir la bande passante bandwidth = bandwidth_factor * freq - lowcut = freq - bandwidth + lowcut = max(0.01, freq - bandwidth) # Éviter les fréquences négatives highcut = freq + bandwidth - - if nyquist: - # Éviter les fréquences négatives ou supérieures à la fréquence de Nyquist - nyq = 0.5 * fs - if lowcut <= 0: - lowcut = 0.01 - if highcut >= nyq: - highcut = nyq - 0.01 - # Créer le filtre Butterworth - lowcut = lowcut / nyq - highcut = highcut / nyq - b, a = butter(4, [lowcut, highcut], btype='band') - - # Appliquer le filtre (filtfilt pour éviter le déphasage) - filtered_signal = filtfilt(b, a, signal_clean) + if nyquist: + # Éviter les fréquences supérieures à la fréquence de Nyquist + nyq = 0.5 * fs + highcut = min(highcut, nyq * 0.99) # Laisser une marge + + # Normaliser les fréquences + lowcut_norm = lowcut / nyq + highcut_norm = highcut / nyq + else: + lowcut_norm = lowcut + highcut_norm = highcut + + print(f"Filtrage à {freq:.2f} Hz: [{lowcut_norm:.3f}, {highcut_norm:.3f}] (normalisé)") + + # Vérifier si le signal est suffisamment long + period = 1 / lowcut if lowcut > 0 else 0 + required_cycles = 40 # Vous avez mentionné 30 cycles + min_length = required_cycles * period * fs + + if len(signal_clean) < min_length: + print(f"Avertissement: Signal court pour {freq:.2f} Hz. " + f"Requis: {min_length:.0f} échantillons, Disponible: {len(signal_clean)}") + + # Méthode alternative: utiliser un filtre SOS (Second-Order Sections) + # plus stable pour les signaux courts + sos = butter(4, [lowcut_norm, highcut_norm], btype='band', output='sos') + filtered_signal = sosfiltfilt(sos, signal_clean) + else: + # Méthode standard + b, a = butter(4, [lowcut_norm, highcut_norm], btype='band') + filtered_signal = filtfilt(b, a, signal_clean) + filtered_signals.append(filtered_signal) # Le temps reste le même après filtrage