upd

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