518 lines
17 KiB
Python
518 lines
17 KiB
Python
import warnings
|
|
import numpy as np
|
|
import pandas as pd
|
|
import pywt
|
|
|
|
from scipy.fft import fft, fftfreq
|
|
from scipy.signal import find_peaks, hilbert, butter, filtfilt, correlate, sosfiltfilt
|
|
from scipy.stats import linregress
|
|
from scipy.optimize import curve_fit
|
|
|
|
|
|
def filter_signal_by_modes(time, signal, num_modes=1, bandwidth_factor=0.1):
|
|
"""
|
|
Filtre un signal pour extraire ses composantes fréquentielles dominantes.
|
|
|
|
Parameters:
|
|
time (array-like): Array contenant les données temporelles.
|
|
signal (array-like): Array contenant les données du signal.
|
|
num_modes (int): Nombre de modes fréquentiels à extraire.
|
|
|
|
Returns:
|
|
tuple: (filtered_signals, frequencies, time_filtered)
|
|
filtered_signals: Liste des signaux filtrés pour chaque mode
|
|
frequencies: Fréquences correspondantes à chaque mode
|
|
time_filtered: Temps après filtrage (même que l'entrée)
|
|
"""
|
|
# Convertir en arrays numpy
|
|
time = np.array(time)
|
|
signal = np.array(signal)
|
|
|
|
# Calculer la période d'échantillonnage et la fréquence
|
|
dt = np.mean(np.diff(time))
|
|
fs = 1 / dt
|
|
|
|
# Supprimer la composante DC
|
|
signal_clean = signal - np.mean(signal)
|
|
|
|
# Calculer la FFT
|
|
n = len(time)
|
|
fft_signal = fft(signal_clean)
|
|
freqs = fftfreq(n, dt)
|
|
|
|
# Prendre seulement les fréquences positives
|
|
positive_mask = freqs > 0
|
|
freqs_pos = freqs[positive_mask]
|
|
fft_pos = fft_signal[positive_mask]
|
|
|
|
# Trouver les pics dans le spectre
|
|
magnitude = np.abs(fft_pos)
|
|
peaks, _ = find_peaks(magnitude, height=0.1*np.max(magnitude))
|
|
peak_freqs = freqs_pos[peaks]
|
|
peak_mags = magnitude[peaks]
|
|
|
|
# Trier les pics par magnitude et sélectionner les num_modes plus importants
|
|
idx = np.argsort(peak_mags)[-num_modes:]
|
|
dominant_freqs = peak_freqs[idx]
|
|
|
|
# Trier les fréquences par ordre croissant
|
|
dominant_freqs.sort()
|
|
|
|
# Initialiser la liste des signaux filtrés
|
|
filtered_signals = []
|
|
|
|
# Pour chaque fréquence dominante, appliquer un filtre passe-bande
|
|
for freq in dominant_freqs:
|
|
# Définir la bande passante
|
|
bandwidth = bandwidth_factor * freq
|
|
lowcut = max(0.01, freq - bandwidth) # Éviter les fréquences négatives
|
|
highcut = freq + bandwidth
|
|
|
|
# É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
|
|
|
|
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
|
|
time_filtered = time
|
|
|
|
return filtered_signals, dominant_freqs, time_filtered
|
|
|
|
|
|
def analyze_signal_Hilbert(time, e_signal, s_signal, freq_rtol=0.01):
|
|
"""
|
|
Analyzes the dominant frequencies, phase shift, time shift, and periods between two signals e_signal(t) and s_signal(t) from a DataFrame.
|
|
Uses Hilbert transform method for phase shift calculation (single mode only).
|
|
"""
|
|
# Extract signals and time
|
|
#time = df[time_column].values
|
|
#e_signal = df[e_signal_column].values
|
|
#s_signal = df[s_signal_column].values
|
|
|
|
# Calculate sampling period and rate
|
|
dt = np.mean(np.diff(time))
|
|
fs = 1 / dt
|
|
|
|
# Remove mean to eliminate DC component
|
|
e_signal_clean = e_signal - np.mean(e_signal)
|
|
s_signal_clean = s_signal - np.mean(s_signal)
|
|
|
|
# Compute FFT
|
|
n = len(time)
|
|
e_fft = fft(e_signal_clean)
|
|
s_fft = fft(s_signal_clean)
|
|
freqs = fftfreq(n, dt)
|
|
|
|
# Get positive frequencies only
|
|
positive_mask = freqs > 0
|
|
freqs_pos = freqs[positive_mask]
|
|
e_fft_pos = e_fft[positive_mask]
|
|
s_fft_pos = s_fft[positive_mask]
|
|
|
|
# Find dominant frequencies using peak detection
|
|
e_magnitude = np.abs(e_fft_pos)
|
|
s_magnitude = np.abs(s_fft_pos)
|
|
|
|
# Find peaks for e signal
|
|
e_peaks, _ = find_peaks(e_magnitude, height=0.1*np.max(e_magnitude))
|
|
e_peak_freqs = freqs_pos[e_peaks]
|
|
e_peak_mags = e_magnitude[e_peaks]
|
|
|
|
# Find peaks for s signal
|
|
s_peaks, _ = find_peaks(s_magnitude, height=0.1*np.max(s_magnitude))
|
|
s_peak_freqs = freqs_pos[s_peaks]
|
|
s_peak_mags = s_magnitude[s_peaks]
|
|
|
|
# Sort peaks by magnitude and select the top one
|
|
e_idx = np.argsort(e_peak_mags)[-1]
|
|
s_idx = np.argsort(s_peak_mags)[-1]
|
|
|
|
freq_e = e_peak_freqs[e_idx]
|
|
freq_s = s_peak_freqs[s_idx]
|
|
|
|
# Check if frequencies match within 1%
|
|
if np.abs(freq_e - freq_s) / ((freq_e + freq_s)/2) > freq_rtol:
|
|
raise ValueError(f"Frequency difference exceeds {freq_rtol*100}%. Phase shift cannot be determined.")
|
|
|
|
# Calculate mean frequency for phase shift calculation
|
|
mean_freq = (freq_e + freq_s) / 2
|
|
|
|
# Create a bandpass filter around the mean frequency
|
|
nyq = 0.5 * fs
|
|
low = (mean_freq - 0.1 * mean_freq) / nyq
|
|
high = (mean_freq + 0.1 * mean_freq) / nyq
|
|
|
|
if low <= 0:
|
|
low = 0.01
|
|
if high >= 1:
|
|
high = 0.99
|
|
|
|
b, a = butter(4, [low, high], btype='band')
|
|
|
|
# Filter both signals
|
|
e_filtered = filtfilt(b, a, e_signal_clean)
|
|
s_filtered = filtfilt(b, a, s_signal_clean)
|
|
|
|
# Apply Hilbert transform to filtered signals
|
|
e_analytic_filtered = hilbert(e_filtered)
|
|
s_analytic_filtered = hilbert(s_filtered)
|
|
|
|
# Extract phase
|
|
e_phase_filtered = np.unwrap(np.angle(e_analytic_filtered))
|
|
s_phase_filtered = np.unwrap(np.angle(s_analytic_filtered))
|
|
|
|
# Calculate phase difference
|
|
phase_diff = s_phase_filtered - e_phase_filtered
|
|
|
|
# Average phase difference over the stable part of the signal
|
|
n_stable = len(phase_diff) // 4
|
|
stable_phase_diff = phase_diff[n_stable:-n_stable]
|
|
|
|
# Calculate mean phase difference
|
|
mean_phase_diff = np.mean(stable_phase_diff)
|
|
|
|
# Normalize to [-π, π]
|
|
mean_phase_diff = np.arctan2(np.sin(mean_phase_diff), np.cos(mean_phase_diff))
|
|
|
|
# Convert to degrees
|
|
phase_diff_deg = np.degrees(mean_phase_diff)
|
|
|
|
# Calculate time shift
|
|
time_shift = mean_phase_diff / (2 * np.pi * mean_freq)
|
|
|
|
# Calculate periods
|
|
period_e = 1 / freq_e
|
|
period_s = 1 / freq_s
|
|
|
|
print(f'phase = {mean_phase_diff}')
|
|
return (period_e, period_s, freq_e, freq_s,
|
|
mean_phase_diff, phase_diff_deg, time_shift)
|
|
|
|
def analyze_signal_sinfit(time, e_signal, s_signal, freq_rtol=0.01):
|
|
"""
|
|
Analyzes the dominant frequencies, phase shift, time shift, and periods between two signals e_signal(t) and s_signal(t) from a DataFrame.
|
|
Uses curve fitting method for phase shift calculation (single mode only).
|
|
"""
|
|
# Extract signals and time
|
|
#time = df[time_column].values
|
|
#e_signal = df[e_signal_column].values
|
|
#s_signal = df[s_signal_column].values
|
|
|
|
# Calculate sampling period and rate
|
|
dt = np.mean(np.diff(time))
|
|
fs = 1 / dt
|
|
|
|
# Remove mean to eliminate DC component
|
|
e_signal_clean = e_signal - np.mean(e_signal)
|
|
s_signal_clean = s_signal - np.mean(s_signal)
|
|
|
|
# Compute FFT
|
|
n = len(time)
|
|
e_fft = fft(e_signal_clean)
|
|
s_fft = fft(s_signal_clean)
|
|
freqs = fftfreq(n, dt)
|
|
|
|
# Get positive frequencies only
|
|
positive_mask = freqs > 0
|
|
freqs_pos = freqs[positive_mask]
|
|
e_fft_pos = e_fft[positive_mask]
|
|
s_fft_pos = s_fft[positive_mask]
|
|
|
|
# Find dominant frequencies using peak detection
|
|
e_magnitude = np.abs(e_fft_pos)
|
|
s_magnitude = np.abs(s_fft_pos)
|
|
|
|
# Find peaks for e signal
|
|
e_peaks, _ = find_peaks(e_magnitude, height=0.1*np.max(e_magnitude))
|
|
e_peak_freqs = freqs_pos[e_peaks]
|
|
e_peak_mags = e_magnitude[e_peaks]
|
|
|
|
# Find peaks for s signal
|
|
s_peaks, _ = find_peaks(s_magnitude, height=0.1*np.max(s_magnitude))
|
|
s_peak_freqs = freqs_pos[s_peaks]
|
|
s_peak_mags = s_magnitude[s_peaks]
|
|
|
|
# Sort peaks by magnitude and select the top one
|
|
e_idx = np.argsort(e_peak_mags)[-1]
|
|
s_idx = np.argsort(s_peak_mags)[-1]
|
|
|
|
freq_e = e_peak_freqs[e_idx]
|
|
freq_s = s_peak_freqs[s_idx]
|
|
|
|
# Check if frequencies match within 1%
|
|
if np.abs(freq_e - freq_s) / ((freq_e + freq_s)/2) > freq_rtol:
|
|
raise ValueError(f"Frequency difference exceeds {freq_rtol*100}%. Phase shift cannot be determined.")
|
|
|
|
# Calculate mean frequency for phase shift calculation
|
|
mean_freq = (freq_e + freq_s) / 2
|
|
|
|
# Define sine function for fitting
|
|
def sine_func(t, A, phi, offset):
|
|
return A * np.sin(2 * np.pi * mean_freq * t + phi) + offset
|
|
|
|
# Fit sine wave to e signal
|
|
try:
|
|
popt_e, _ = curve_fit(sine_func, time, e_signal,
|
|
p0=[np.std(e_signal), 0, np.mean(e_signal)],
|
|
bounds=([0, -np.pi, -np.inf], [np.inf, np.pi, np.inf]))
|
|
except:
|
|
popt_e = [np.std(e_signal), 0, np.mean(e_signal)]
|
|
|
|
# Fit sine wave to s signal
|
|
try:
|
|
popt_s, _ = curve_fit(sine_func, time, s_signal,
|
|
p0=[np.std(s_signal), 0, np.mean(s_signal)],
|
|
bounds=([0, -np.pi, -np.inf], [np.inf, np.pi, np.inf]))
|
|
except:
|
|
popt_s = [np.std(s_signal), 0, np.mean(s_signal)]
|
|
|
|
# Calculate phase difference
|
|
phase_diff = popt_s[1] - popt_e[1]
|
|
|
|
# Normalize to [-π, π]
|
|
phase_diff = np.arctan2(np.sin(phase_diff), np.cos(phase_diff))
|
|
|
|
# Convert to degrees
|
|
phase_diff_deg = np.degrees(phase_diff)
|
|
|
|
# Calculate time shift
|
|
time_shift = phase_diff / (2 * np.pi * mean_freq)
|
|
|
|
# Calculate periods
|
|
period_e = 1 / freq_e
|
|
period_s = 1 / freq_s
|
|
|
|
print(f'phase = {phase_diff}')
|
|
return (period_e, period_s, freq_e, freq_s,
|
|
phase_diff, phase_diff_deg, time_shift)
|
|
|
|
|
|
|
|
def analyze_signal_cross_correlation(time, e_signal, s_signal, freq_rtol=0.01):
|
|
"""
|
|
Analyzes signals using cross-correlation method for phase shift calculation (single mode only).
|
|
"""
|
|
# Extract signals and time
|
|
#time = df[time_column].values
|
|
#e_signal = df[e_signal_column].values
|
|
#s_signal = df[s_signal_column].values
|
|
|
|
# Calculate sampling period and rate
|
|
dt = np.mean(np.diff(time))
|
|
fs = 1 / dt
|
|
|
|
# Remove mean to eliminate DC component
|
|
e_signal_clean = e_signal - np.mean(e_signal)
|
|
s_signal_clean = s_signal - np.mean(s_signal)
|
|
|
|
# Compute FFT
|
|
n = len(time)
|
|
e_fft = fft(e_signal_clean)
|
|
s_fft = fft(s_signal_clean)
|
|
freqs = fftfreq(n, dt)
|
|
|
|
# Get positive frequencies only
|
|
positive_mask = freqs > 0
|
|
freqs_pos = freqs[positive_mask]
|
|
e_fft_pos = e_fft[positive_mask]
|
|
s_fft_pos = s_fft[positive_mask]
|
|
|
|
# Find dominant frequencies using peak detection
|
|
e_magnitude = np.abs(e_fft_pos)
|
|
s_magnitude = np.abs(s_fft_pos)
|
|
|
|
# Find peaks for e signal
|
|
e_peaks, _ = find_peaks(e_magnitude, height=0.1*np.max(e_magnitude))
|
|
e_peak_freqs = freqs_pos[e_peaks]
|
|
e_peak_mags = e_magnitude[e_peaks]
|
|
|
|
# Find peaks for s signal
|
|
s_peaks, _ = find_peaks(s_magnitude, height=0.1*np.max(s_magnitude))
|
|
s_peak_freqs = freqs_pos[s_peaks]
|
|
s_peak_mags = s_magnitude[s_peaks]
|
|
|
|
# Sort peaks by magnitude and select the top one
|
|
e_idx = np.argsort(e_peak_mags)[-1]
|
|
s_idx = np.argsort(s_peak_mags)[-1]
|
|
|
|
freq_e = e_peak_freqs[e_idx]
|
|
freq_s = s_peak_freqs[s_idx]
|
|
|
|
# Check if frequencies match within 1%
|
|
if np.abs(freq_e - freq_s) / ((freq_e + freq_s)/2) > freq_rtol:
|
|
raise ValueError(f"Frequency difference exceeds {freq_rtol*100}%. Phase shift cannot be determined.")
|
|
|
|
# Calculate mean frequency for phase shift calculation
|
|
mean_freq = (freq_e + freq_s) / 2
|
|
|
|
# Create a bandpass filter around the mean frequency
|
|
nyq = 0.5 * fs
|
|
low = (mean_freq - 0.1 * mean_freq) / nyq
|
|
high = (mean_freq + 0.1 * mean_freq) / nyq
|
|
|
|
if low <= 0:
|
|
low = 0.01
|
|
if high >= 1:
|
|
high = 0.99
|
|
|
|
b, a = butter(4, [low, high], btype='band')
|
|
|
|
# Filter both signals
|
|
e_filtered = filtfilt(b, a, e_signal_clean)
|
|
s_filtered = filtfilt(b, a, s_signal_clean)
|
|
|
|
# Compute cross-correlation
|
|
correlation = correlate(s_filtered, e_filtered, mode='full')
|
|
lags = np.arange(-len(e_filtered) + 1, len(e_filtered))
|
|
|
|
# Find the lag with maximum correlation
|
|
max_lag = lags[np.argmax(correlation)]
|
|
|
|
# Calculate time shift
|
|
time_shift = max_lag * dt
|
|
|
|
# Calculate phase shift
|
|
phase_diff = 2 * np.pi * mean_freq * time_shift
|
|
|
|
# Normalize to [-π, π]
|
|
phase_diff = np.arctan2(np.sin(phase_diff), np.cos(phase_diff))
|
|
|
|
# Convert to degrees
|
|
phase_diff_deg = np.degrees(phase_diff)
|
|
|
|
# Calculate periods
|
|
period_e = 1 / freq_e
|
|
period_s = 1 / freq_s
|
|
|
|
print(f'phase = {phase_diff}')
|
|
return (period_e, period_s, freq_e, freq_s,
|
|
phase_diff, phase_diff_deg, time_shift)
|
|
|
|
|
|
def analyze_signal_wavelet(time, e_signal, s_signal, freq_rtol=0.01):
|
|
"""
|
|
Analyzes signals using wavelet transform for phase shift calculation (single mode only).
|
|
"""
|
|
# Extract signals and time
|
|
#time = df[time_column].values
|
|
#e_signal = df[e_signal_column].values
|
|
#s_signal = df[s_signal_column].values
|
|
|
|
# Calculate sampling period and rate
|
|
dt = np.mean(np.diff(time))
|
|
fs = 1 / dt
|
|
|
|
# Remove mean to eliminate DC component
|
|
e_signal_clean = e_signal - np.mean(e_signal)
|
|
s_signal_clean = s_signal - np.mean(s_signal)
|
|
|
|
# Compute FFT
|
|
n = len(time)
|
|
e_fft = fft(e_signal_clean)
|
|
s_fft = fft(s_signal_clean)
|
|
freqs = fftfreq(n, dt)
|
|
|
|
# Get positive frequencies only
|
|
positive_mask = freqs > 0
|
|
freqs_pos = freqs[positive_mask]
|
|
e_fft_pos = e_fft[positive_mask]
|
|
s_fft_pos = s_fft[positive_mask]
|
|
|
|
# Find dominant frequencies using peak detection
|
|
e_magnitude = np.abs(e_fft_pos)
|
|
s_magnitude = np.abs(s_fft_pos)
|
|
|
|
# Find peaks for e signal
|
|
e_peaks, _ = find_peaks(e_magnitude, height=0.1*np.max(e_magnitude))
|
|
e_peak_freqs = freqs_pos[e_peaks]
|
|
e_peak_mags = e_magnitude[e_peaks]
|
|
|
|
# Find peaks for s signal
|
|
s_peaks, _ = find_peaks(s_magnitude, height=0.1*np.max(s_magnitude))
|
|
s_peak_freqs = freqs_pos[s_peaks]
|
|
s_peak_mags = s_magnitude[s_peaks]
|
|
|
|
# Sort peaks by magnitude and select the top one
|
|
e_idx = np.argsort(e_peak_mags)[-1]
|
|
s_idx = np.argsort(s_peak_mags)[-1]
|
|
|
|
freq_e = e_peak_freqs[e_idx]
|
|
freq_s = s_peak_freqs[s_idx]
|
|
|
|
# Check if frequencies match within 1%
|
|
if np.abs(freq_e - freq_s) / ((freq_e + freq_s)/2) > freq_rtol:
|
|
raise ValueError(f"Frequency difference exceeds {freq_rtol*100}%. Phase shift cannot be determined.")
|
|
|
|
# Calculate mean frequency for phase shift calculation
|
|
mean_freq = (freq_e + freq_s) / 2
|
|
|
|
# Perform continuous wavelet transform
|
|
scales = np.arange(1, 128)
|
|
coefficients_e, frequencies_e = pywt.cwt(e_signal_clean, scales, 'morl', sampling_period=dt)
|
|
coefficients_s, frequencies_s = pywt.cwt(s_signal_clean, scales, 'morl', sampling_period=dt)
|
|
|
|
# Find the scale index closest to the mean frequency
|
|
scale_idx = np.argmin(np.abs(frequencies_e - mean_freq))
|
|
|
|
# Extract the wavelet coefficients at this scale
|
|
coeffs_e = coefficients_e[scale_idx, :]
|
|
coeffs_s = coefficients_s[scale_idx, :]
|
|
|
|
# Calculate the instantaneous phase using the Hilbert transform
|
|
analytic_e = hilbert(coeffs_e)
|
|
analytic_s = hilbert(coeffs_s)
|
|
|
|
phase_e = np.unwrap(np.angle(analytic_e))
|
|
phase_s = np.unwrap(np.angle(analytic_s))
|
|
|
|
# Calculate phase difference
|
|
phase_diff = phase_s - phase_e
|
|
|
|
# Average phase difference over the stable part of the signal
|
|
n_stable = len(phase_diff) // 4
|
|
stable_phase_diff = phase_diff[n_stable:-n_stable]
|
|
|
|
# Calculate mean phase difference
|
|
mean_phase_diff = np.mean(stable_phase_diff)
|
|
|
|
# Normalize to [-π, π]
|
|
mean_phase_diff = np.arctan2(np.sin(mean_phase_diff), np.cos(mean_phase_diff))
|
|
|
|
# Convert to degrees
|
|
phase_diff_deg = np.degrees(mean_phase_diff)
|
|
|
|
# Calculate time shift
|
|
time_shift = mean_phase_diff / (2 * np.pi * mean_freq)
|
|
|
|
# Calculate periods
|
|
period_e = 1 / freq_e
|
|
period_s = 1 / freq_s
|
|
|
|
print(f'phase = {mean_phase_diff}')
|
|
return (period_e, period_s, freq_e, freq_s,
|
|
mean_phase_diff, phase_diff_deg, time_shift)
|