# A COMPLETER
# Import des librairies
import numpy as np
from matplotlib import pyplot as plt

# A COMPLETER
# Développement de la fonction de convolution

# Les fonctions de padding et de shift restent inchangées par rapport à la corrélation
def padding(signal,M):
    """
    Fonction de padding du signal pour l'agrandir à la taille M (si M supérieur à len(signal))
    signal : tableau 1D NumPy, représentant un signal
    M : entier, représentat la taille de padding attendue pour M
    """
    if M<=len(signal):
        return signal
    out = np.zeros((M),dtype=signal.dtype)
    out[:len(signal)] = signal
    return out

def shift(x, S):
    """
    Effectue un décalage de S sur le signal x
    x : tableau 1D NumPy, représentant un signal
    S : entier, représentant le pas de décalage à effectuer
    """
    len_x = len(x)
    y = np.zeros((len_x+abs(S)))
    if S>0:
        y[-len_x:] = x
        return y[:len_x]
    else:
        y[:len_x] = x
        return y[-len_x:]

# Quatre changements à faire dans la fonction de convolution calculées avec un décalage de n cases :
# - Le décalage temporel ne se fait pas sur f mais sur g
# - Le décalage S à effectuer devient S - (len(g)-1)
# - Le signal g est renversé (g -> g[::-1])
# - On n'utilise pas np.conjugate(g) contrairement à la corrélation

def convolve_at_shift_n(f, g, n, Te):
    """
    Calcule l'échantillon n de la convolution de f et g.
    f,g : tableaux 1D NumPy, signaux à convoluer ensemble
    n : entier, pas de décalage à effectuer
    Te : float, temps d'échantillonnage de f et g
    """

    N, M = len(f), len(g)

    # Décalage équivalent pour faire glisser g renversé contre f :
    S = n - (M - 1)

    # On pad comme dans score_correlate pour garder la même logique d’alignement
    len_f, len_g = len(f), len(g)
    f_pad = padding(f, len_g)
    
    # IMPORTANT: convolution -> pas de conjugaison, mais on renverse g
    g_rev = g[::-1]
    g_pad = padding(g_rev, len_f)
    
    # On décale g de S (positif: droite, négatif: gauche)
    g_decale = shift(g_pad, S)

    # Produit scalaire final
    return np.dot(f_pad, g_decale) * Te

# Un changement à faire par rapport au calcul de la corrélation :
# - L'axe temporel évolue de 0 à (len(f)+len(g)-2)*Te secondes
def convolution(f, g, Te):
    """
    Convolution linéaire de f et g.
    f,g : tableaux 1D NumPy, signaux à convoluer ensemble
    Te : float, temps d'échantillonnage de f et g
    """
    N, M = len(f), len(g)
    # Le retour de la convolution est un signal de taille L = len(f)+len(g)-1
    L = N + M - 1
    n_vals = np.arange(L)
    y = np.array([convolve_at_shift_n(f, g, n, Te) for n in n_vals], dtype=np.result_type(f, g, np.float64))
    # L'axe temporel évolue de 0 à L x Te secondes
    t = n_vals * Te
    return t, y

# A COMPLETER
# Création et affichage des deux signaux

# Axe temporel
Fe = 100
Te = 1/Fe
t0 = 2

t = np.arange(0, 5, Te)

# Signaux
x = np.zeros_like(t)
x[np.abs(t-t0)<=0.5] = 1
h = np.exp(-t)                   

# Affichage des signaux
plt.figure()
plt.plot(t,x,label='x(t)')
plt.plot(t,h,label='h(t)')
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Calcul de y(t) = x(t)*h(t) et affichage du résultat

t_y, y = convolution(x,h,Te)

plt.figure()
plt.plot(t_y, y, label="y(t)")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Calcul de la convolution avec NumPy et affichage des résultats de convolution
y_np = np.convolve(x, h, mode='full') * Te # Il ne faut pas oublier de multiplier par Te

plt.figure()
plt.plot(t_y, y, label="y(t) (manuel)")
plt.plot(t_y, y_np, label="y(t) (NumPy)")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Création et affichage du signal triangulaire
Fe = 100
Te = 1/Fe
# Axe temporel
t_x = np.arange(0,2,Te)

# Signal triangulaire
x = 1-np.abs(t_x-1)

# Affichage
plt.plot(t_x,x,label="x(t)")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Création de h(t) et affichage

T = 2

# Axe temporel
t_h = np.arange(0,10.01,Te)

# Signal h(t) peigne de Dirac
h = np.zeros_like(t_h)
h[t_h%T==0] = 1

# Affichage de h(t)
plt.plot(t_h,h,label="h(t)")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Convolution entre x(t) et h(t) et affichage du résultat
t_y, y = convolution(x, h, Te)

plt.plot(t_y,y,label="y(t)")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Modification de h(t) pour corriger y(t)

t_h = np.arange(0,10.01,Te)
h = np.zeros_like(t_h)
h[t_h%T==0] = Fe # Afin de garantir une amplitude correcte en sortie, il faut multiplier le peigne de Dirac par Fe (ou diviser par Te)

plt.plot(t_h,h,label="h(t)")
plt.legend()
plt.tight_layout()
plt.show()

# Convolution entre x(t) et h(t) avec T = 2 secondes
t_y, y = convolution(x, h, Te)

plt.plot(t_y,y,label="y(t)")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Modification de h(t) avec T = 3 secondes
T = 3
t_h = np.arange(0,10.01,Te)
h = np.zeros_like(t_h)
h[t_h%T==0] = Fe

plt.figure()
plt.plot(t_h,h,label="h(t)")
plt.legend()
plt.tight_layout()
plt.show()

# Convolution entre x(t) et h(t) avec T = 3 secondes
t_y, y = convolution(x, h, Te)

plt.figure()
plt.plot(t_y,y,label="y(t)")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Modification de h(t) avec T = 1.5 secondes
T = 1.5
t_h = np.arange(0,10.01,Te)
h = np.zeros_like(t_h)
h[t_h%T==0] = Fe

plt.figure()
plt.plot(t_h,h,label="h(t)")
plt.legend()
plt.tight_layout()
plt.show()

# Convolution entre x(t) et h(t) avec T = 1.5 secondes
t_y, y = convolution(x, h, Te)

plt.figure()
plt.plot(t_y,y,label="y(t)")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Création de x(t) et h(t) et affichage des signaux

Fe = 1000
Te = 1/Fe
f0 = 2
f1 = 5

# Création de l'axe temporel
t = np.arange(0,5,Te)

# Création des signaux
x = np.sin(2*np.pi*f0*t)
h = np.sin(2*np.pi*f1*t)

# Affichage des signaux
plt.subplots(1,2,figsize=(12,4))

plt.subplot(121)
plt.plot(t,x)
plt.title("x(t)")

plt.subplot(122)
plt.plot(t,h)
plt.title("h(t)")

plt.tight_layout()
plt.show()

# A COMPLETER
# FFT réelle de x(t) et h(t) et affichage des spectres d'amplitude

# Calcul de la FFT réelle des signaux
N, M = len(x), len(h)
X = np.fft.rfft(x) 
H = np.fft.rfft(h) 

# Calcul des axes fréquentiels des FFT
freq_X = np.fft.rfftfreq(N,Te)
freq_H = np.fft.rfftfreq(M,Te) 

# Récupération des amplitudes de X et H
amplitude_X = np.abs(X) * 2 / N
amplitude_H = np.abs(H) * 2 / M

# Affichage des spectres d'amplitude
plt.subplots(1,2,figsize=(12,4))

plt.subplot(121)
plt.stem(freq_X,amplitude_X, markerfmt=" ")
plt.title("Spectre d'amplitude de X(f)")

plt.subplot(122)
plt.stem(freq_H,amplitude_H, markerfmt=" ")
plt.title("Spectre d'amplitude de H(f)")

plt.tight_layout()
plt.show()

# A COMPLETER
# Convolution entre x(t) et h(t) et affichage du résultat
t_y, y = convolution(x, h, Te)

plt.plot(t_y,y,label="y(t)")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Vérification du théorème de Plancherel
N, M = len(x), len(h)
L = N + M - 1

# FFT réelles de x et h
X = np.fft.rfft(x, n=L)
H = np.fft.rfft(h, n=L)

# Produit entre X et H
Y = X*H

# FFt réelle inverse de Y
y_freq = np.fft.irfft(Y, n=L) * Te

# Affichage des résultats
plt.plot(t_y,y,label="y(t) via convolution")
plt.plot(t_y,y_freq, label="y(t) via produit des FFT")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Convolution entre x(t) et h(t) avec f1 = [5, 10, 20, 50] Hz
f1_valeurs = [5,10,20,50]

for f1 in f1_valeurs:
    h = np.sin(2*np.pi*f1*t)
    t_y, y = convolution(x, h, Te)
    plt.plot(t_y,y,label=f"y(t) (f1 = {f1} Hz)")

plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Chargement de l'extrait audio
from scipy.io import wavfile
freq_audio, extrait_audio = wavfile.read("audio.wav")
print(f"Fréquence d'échantillonage de l'audio : {freq_audio} Hz")

# Reconstruction de l'axe temporel
t_extrait_audio = np.arange(0, len(extrait_audio)/freq_audio, 1/freq_audio)

# Tracé du signal audio
plt.plot(t_extrait_audio, extrait_audio)
plt.title("Audio")
plt.show()

/tmp/ipykernel_23469/3971863886.py:4: WavFileWarning: Chunk (non-data) not understood, skipping it.
  freq_audio, extrait_audio = wavfile.read("audio.wav")

Fréquence d'échantillonage de l'audio : 44100 Hz

# A COMPLETER
# Chargement de la réponse impulsionnelle
freq_IR, extrait_IR = wavfile.read("eglise.wav")
print(f"Fréquence d'échantillonage de l'audio : {freq_IR} Hz")

# Reconstruction de l'axe temporel
t_extrait_IR = np.arange(0, len(extrait_IR)/freq_IR, 1/freq_IR)

# Tracé de l'IR
plt.plot(t_extrait_IR, extrait_IR)
plt.title("Réponse impulsionnelle")
plt.show()

Fréquence d'échantillonage de l'audio : 44100 Hz

# A COMPLETER
# Normalisation de l'extrait audio et de la réponse impulsionnelle
extrait_audio_normalise = extrait_audio.astype(float)/32768.0
extrait_IR_normalise = extrait_IR.astype(float)/32768.0

# A COMPLETER
# Convolution de l'extrait audio normalisé et de la réponse impulsionnelle normalisée
extrait_audio_reverb = np.convolve(extrait_audio_normalise, extrait_IR_normalise, mode="full")

# A COMPLETER
# Correction de pics trop élevés par normalisation
peak = np.max(np.abs(extrait_audio_reverb))
if peak>0.999:
    extrait_audio_reverb = extrait_audio_reverb/peak*0.999

# A COMPLETER
# Conversion de l'audio en entier 16 bits
extrait_audio_reverb = (extrait_audio_reverb*32767.0).astype(np.int16)

# A COMPLETER
# Affichage du spectre d'audio résultat

# Pour retrouver l'axe temporel, étant donné que notre audio n'a pas changé de fréquences d'échantillonage, il suffit de récupérer la longueur du tableau NumPy de l'audio et de le diviser par la fréquence d'échantillonage
t_audio_reverb = np.arange(0,len(extrait_audio_reverb)/freq_audio,1/freq_audio)
plt.plot(t_audio_reverb, extrait_audio_reverb)
plt.title("Audio avec reverb")
plt.show()

# A COMPLETER
# Enregistrement du son final avec reverb
wavfile.write("audio_reverb.wav",freq_audio,(extrait_audio_reverb).astype(np.int16))

Traitement du Signal - TP4 : Convolution¶

Exercice 1 : Pour faire des convolutions, il faut une fonction...¶

Exercice 2 : J'adore mon signal, je voudrais l'avoir tout le temps !¶

Exercice 3 : Plancherel ?! J'ai déjà entendu cette théorie...¶

Exercice 4 : HEEEEEEE OOOOHHHHH¶

La convolution au service de l'art¶