Vérification des performances d'un filtre¶

import numpy as np
import matplotlib.pyplot as plt

# Création du signal bruité à filtrer

# Fréquence d'échantillonnage
fe = 110e3  # Hz

# Signal bruité
T = 0.05  # Durée du signal
nEchantillons = int(T * fe)  # Nombre d'échantillons
t = np.linspace(0, T, nEchantillons, endpoint=False)
a = 0.02  # Amplitude du signal à retrouver
f0 = 600.0  # Fréquence du signal à retrouver
x = 0.1 * np.sin(2 * np.pi * 1.5 * np.sqrt(t))
x += 0.01 * np.cos(2 * np.pi * 314 * t + 0.1)
x += a * np.cos(2 * np.pi * f0 * t + .12)
x += 0.03 * np.cos(2 * np.pi * 2000 * t)

plt.plot(t, x, label="Signal bruité")
plt.xlabel("Temps(s)")
plt.ylabel("Amplitude")
plt.title("Signal bruité")
plt.grid()
plt.show()

from numpy.fft import fft, fftfreq

# Calcul FFT
X = fft(x)  # Transformée de fourier
freq = fftfreq(x.size, d=1/fe)  # Fréquences de la transformée de Fourier

# Calcul du nombre d'échantillon
N = x.size

# On prend la valeur absolue de l'amplitude uniquement pour les fréquences positives et normalisation
X_abs = np.abs(X[:N//2])*2.0/N
# On garde uniquement les fréquences positives
freq_pos = freq[:N//2]

plt.plot(freq_pos, X_abs, label="Amplitude absolue")
plt.xlim(0, 2500)  # On réduit la plage des fréquences à la zone utile
plt.grid()
plt.title("Analyse fréquentielle du signal")
plt.xlabel(r"Fréquence (Hz)")
plt.ylabel(r"Amplitude $|X(f)|$"),
plt.show()

from scipy import signal

# Fonction de création du filtre passe-bande de Butterworth
def butter_passeBande(basseFreq, hauteFreq, fe, ordre=3):
    f_nyq = 0.5 * fe  # Fréquence de nyquist
    f_basse = basseFreq / f_nyq
    f_haute = hauteFreq / f_nyq
    b, a = signal.butter(ordre, [f_basse, f_haute], btype='band', analog=False)
    return b, a


# Fréquence de coupure
fc1 = 400  # Hz
fc2 = 800  # Hz

# Test pour différents ordres du filtre de butterworth et tracé de la réponse en fréquence du filtre
for ordre in [2, 3, 4, 5, 6]:
    b, a = butter_passeBande(fc1, fc2, fe, ordre=ordre)
    w, h = signal.freqz(b, a, worN=2000)
    plt.plot((fe * 0.5 / np.pi) * w, abs(h),
             label="ordre = {:d}".format(ordre))

plt.plot([0, 0.5 * fe], [np.sqrt(0.5), np.sqrt(0.5)],
         '--', label='Valeur limite')
plt.xlabel('Fréquence (Hz)')
plt.xlim(0, 2*fc2)
plt.ylabel('Gain')
plt.grid(True)
plt.legend(loc='best')
plt.title("Réponse en fréquence du filtre pour différents ordres")
plt.show()

# Application du filtre passe-bande en mode causal
b, a = butter_passeBande(fc1, fc2, fe, ordre=4)
s_but_bp = signal.lfilter(b, a, x)

# Application du filtre passe-bande en mode double passe, acausal
s_but_bp2 = signal.filtfilt(b, a, x)

plt.plot(t, x, color='silver', label='Signal brut')
plt.plot(t, s_but_bp, label='Signal filtré causal')
plt.plot(t, s_but_bp2, label='Signal filtré double passe', alpha=0.7)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Filtrage du signal")
plt.legend()
plt.show()

def fir_passeBande(basseFreq, hauteFreq, fe, largeur=10, attenuation=20.0):
    f_nyq = 0.5 * fe  # Fréquence de nyquist

    # La largeur relative à la fréquence du passage de l'état passant à bloqué en Hz
    largeur = largeur/f_nyq

    # Calcul l'ordre du filtre FIR.
    # Attenuation en dB du filtre
    N, beta = signal.kaiserord(attenuation, largeur)
    #plt.plot(taps, 'bo-', linewidth=2)
    #plt.title('Coefficient du filtre (%d taps)' % N)
    # plt.grid(True)
    # plt.show()

    # Fréquences de coupure
    f_basse = basseFreq / f_nyq
    f_haute = hauteFreq / f_nyq

    # Création du filtre
    if N % 2 == 0:  # N doit être impair pour un passe-bande
        N += 1
    c = signal.firwin(N, [f_basse, f_haute],
                      pass_zero='bandpass', window=('kaiser', 5))

    return c

# Fréquence de coupure
fc1 = 400  # Hz
fc2 = 1000  # Hz

# Création du filtre passe-bande
c = fir_passeBande(fc1, fc2, fe, largeur=50)

# Réponse en fréquence du filtre
w, h = signal.freqz(c, worN=4000)

plt.plot((w/np.pi)*fe/2, abs(h), linewidth=2)
plt.xlabel('Fréquence (Hz)')
plt.ylabel('Gain')
plt.title('Réponse en fréquence')
plt.ylim(-0.05, 1.05)
plt.xlim(0, 2*fc2)
plt.grid(True)

# Application du filtre passe-bande en mode causal
s_fir_bp = signal.lfilter(c, 1.0, x)

plt.plot(t, x, color='silver', label='Signal brut')
plt.plot(t, s_fir_bp, label='Signal filtré FIR')
plt.plot(t, s_but_bp, label='Signal filtré Butterworth', alpha=0.5)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Filtrage du signal")
plt.legend()
plt.show()