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

# A COMPLETER
# Modification de la fonction somme_sinus

def somme_sinus(parameters,T,fe):
    """
    Construit un signal représentant une somme de sinus
    parameters: liste [frequence, amplitude, phase], contenant les paramètres de chaque sinus à intégrer dans la somme
    T: entier/float, le temps maximum d'évolution du signal
    fe : entier/float, la fréquence d'échantillonnage du signal
    """
    # Création de l'axe temporel
    temps = np.arange(0,T,1/fe)
    
    # Initialisation du signal
    signal = np.zeros_like(temps)
    
    # On itère sur la liste de paramètres
    for p in parameters:
        freq, amp, phase = p
        signal+= amp*np.sin(2*np.pi*freq*temps + phase)
    return temps, signal

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

parameters_s = [[12, 1, -np.pi/2]]

t,s = somme_sinus(parameters_s, 1, 2048)

# Affichage du signal
plt.figure()
plt.plot(t, s)
plt.title("Signal s(t)")
plt.tight_layout()
plt.show()

# A COMPLETER
# Visualisation de s(t) à différentes fréquences d'échantillonnage

K = 8 
fe_init = 2048 # Fréquence d'échantillonnage initiale
T = 1 # Temps d'évolution du signal en secondes

plt.subplots(K//2,2,figsize=(12,16))
for i in range(K):
    fe =  fe_init//(2**i)
    t,s = somme_sinus(parameters_s, T, fe)

    plt.subplot(K//2,2,i+1)
    plt.plot(t, s)
    plt.title(f"Signal s(t), fe = {fe} Hz")
plt.tight_layout()
plt.show()

# A COMPLETER
# Visualisation des spectres d'amplitude de s(t) à différentes fréquences d'échantillonnage
from scipy.signal import find_peaks

K = 8 
fe_init = 2048 # Fréquence d'échantillonnage initiale
T = 1 # Temps d'évolution du signal en secondes

plt.subplots(K//2,2,figsize=(12,16))
for i in range(K):
    fe =  fe_init//(2**i)
    t,s = somme_sinus(parameters_s, T, fe)
    S = np.fft.fftshift(np.fft.fft(s))
    f = np.fft.fftshift(np.fft.fftfreq(len(s),1/fe))
    amp = np.abs(S)*1/len(s)


    plt.subplot(K//2,2,i+1)
    plt.stem(f, amp, markerfmt=' ')
    plt.title(f"Spectre d'amplitude de S(f), fe = {fe} Hz")

    peaks, properties = find_peaks(amp, height=np.max(amp)*0.01)
    print(f'fe = {fe} hz : {[float(f[p]) for p in peaks]}')

plt.tight_layout()
plt.show()

fe = 2048 hz : [-12.0, 12.0]
fe = 1024 hz : [-12.0, 12.0]
fe = 512 hz : [-12.0, 12.0]
fe = 256 hz : [-12.0, 12.0]
fe = 128 hz : [-12.0, 12.0]
fe = 64 hz : [-12.0, 12.0]
fe = 32 hz : [-12.0, 12.0]
fe = 16 hz : [-4.0, 4.0]

# A COMPLETER
# Reconstruction de s(t) à différentes fréquences d'échantillonnage pour obtenir fe' = 2048 Hz

K = 8 
fe_init = 2048 # Fréquence d'échantillonnage initiale
T = 1 # Temps d'évolution du signal en secondes
N_prime = fe_init*T # Nombre de points à avoir après zero-padding

plt.subplots(K//2,2,figsize=(12,16))
for i in range(K):

    fe =  fe_init//(2**i)
    t, s = somme_sinus(parameters_s, T, fe)
    S = np.fft.fftshift(np.fft.fft(s))
    f = np.fft.fftshift(np.fft.fftfreq(len(s),1/fe))
    
    # Zero-padding centré
    N = len(s) # Nombre de points de la TF du signal
    S_pad = np.pad(S,((N_prime-N)//2,(N_prime-N)//2),constant_values=(0,0))
    
    # Correction de l'amplitude
    S_pad *= N_prime/N

    # TF inverse (attention à ne pas oublier ifftshift !)
    s_ = np.fft.ifft(np.fft.ifftshift(S_pad))
    
    # Affichage du résultat
    plt.subplot(K//2,2,i+1)
    # Ici, on définit l'axe temporel t_2048, allant de 0 à 1 secondes avec fe = 2048 Hz
    t_2048 = np.arange(0,1,1/2048)
    plt.plot(t_2048, np.real(s_))
    plt.title(f"Signal s(t) reconstruit à 2048 Hz à partir de fe = {fe} Hz")
plt.tight_layout()
plt.show()

# A COMPLETER
# Construction et affichage du signal x(t)

fe = 2000
T = 1
parameters_x = [[300,1,-np.pi/2],
                [400,1,-np.pi/2],
                [1300,1,-np.pi/2],
                [3600,1,-np.pi/2]]

t,x = somme_sinus(parameters_x, T, fe)

# Affichage du signal
plt.figure()
plt.plot(t, x)
plt.title("Signal x(t)")
plt.tight_layout()
plt.show()

# A COMPLETER
# Transformée de Fourier de x(t) et affichage du spectre d'amplitude
X = np.fft.fftshift(np.fft.fft(x))
f = np.fft.fftshift(np.fft.fftfreq(len(x),1/fe))
amp = np.abs(X)*1/len(x)

plt.stem(f, amp, markerfmt=' ')
plt.title(f"Spectre d'amplitude de X(f), fe = {fe} Hz")
plt.show()

# A COMPLETER
# Détection des pics et de leur amplitude
peaks, properties = find_peaks(amp, height=np.max(amp)*0.01)
print(f'Pics détectés: ')
for peak, height in zip(peaks, properties['peak_heights']):
    print(f"• {f[peak]} Hz - Amplitude {height}")

Pics détectés: 
• -700.0 Hz - Amplitude 0.5000000000000009
• -400.0 Hz - Amplitude 1.000000000000002
• -300.0 Hz - Amplitude 0.5000000000000033
• 300.0 Hz - Amplitude 0.5000000000000033
• 400.0 Hz - Amplitude 1.000000000000002
• 700.0 Hz - Amplitude 0.5000000000000009

# A COMPLETER
# Filtre passe-bas anti-repliement appliqué à x(t), affichage du spectre d'amplitude modifié et détection des pics (fréquence, amplitude)
fc = 900
X[np.abs(f)>fc]=0

amp = np.abs(X)*1/len(x)
plt.stem(f, amp, markerfmt=' ')
plt.title(f"Spectre d'amplitude de X(f), fe = {fe} Hz, fc = {fc} Hz")
plt.show()

peaks, properties = find_peaks(amp, height=np.max(amp)*0.01)
print(f'Pics détectés: ')
for peak, height in zip(peaks, properties['peak_heights']):
    print(f"• {f[peak]} Hz - Amplitude {height}")

Pics détectés: 
• -700.0 Hz - Amplitude 0.5000000000000009
• -400.0 Hz - Amplitude 1.000000000000002
• -300.0 Hz - Amplitude 0.5000000000000033
• 300.0 Hz - Amplitude 0.5000000000000033
• 400.0 Hz - Amplitude 1.000000000000002
• 700.0 Hz - Amplitude 0.5000000000000009

# A COMPLETER
# Mêmes opérations mais avec fe = 2kHz et fc = 1500 Hz

fe = 2000
fc = 1500

T = 1
parameters_x = [[300,1,-np.pi/2],
                [400,1,-np.pi/2],
                [1300,1,-np.pi/2],
                [3600,1,-np.pi/2]]

t,x = somme_sinus(parameters_x, T, fe)

# Affichage du signal
plt.figure()
plt.plot(t, x)
plt.title("Signal x(t)")
plt.tight_layout()
plt.show()

# Transformée de Fourier de x(t) et affichage du spectre d'amplitude
X = np.fft.fftshift(np.fft.fft(x))
f = np.fft.fftshift(np.fft.fftfreq(len(x),1/fe))
amp = np.abs(X)*1/len(x)

plt.stem(f, amp, markerfmt=' ')
plt.title(f"Spectre d'amplitude de X(f), fe = {fe} Hz")
plt.show()

# Filtre passe-bas idéal appliqué à x(t), affichage du spectre d'amplitude modifié et détection des pics (fréquence, amplitude)
X[np.abs(f)>fc]=0

amp = np.abs(X)*1/len(x)
plt.stem(f, amp, markerfmt=' ')
plt.title(f"Spectre d'amplitude de X(f), fe = {fe} Hz, fc = {fc} Hz")
plt.show()

peaks, properties = find_peaks(amp, height=np.max(amp)*0.01)
print(f'Pics détectés: ')
for peak, height in zip(peaks, properties['peak_heights']):
    print(f"• {f[peak]} Hz - Amplitude {height}")

Pics détectés: 
• -700.0 Hz - Amplitude 0.5000000000000009
• -400.0 Hz - Amplitude 1.000000000000002
• -300.0 Hz - Amplitude 0.5000000000000033
• 300.0 Hz - Amplitude 0.5000000000000033
• 400.0 Hz - Amplitude 1.000000000000002
• 700.0 Hz - Amplitude 0.5000000000000009

# A COMPLETER
# Mêmes opérations mais avec fe = 4kHz et fc = 500 Hz

fe = 4000
fc = 500

T = 1
parameters_x = [[300,1,-np.pi/2],
                [400,1,-np.pi/2],
                [1300,1,-np.pi/2],
                [3600,1,-np.pi/2]]

t,x = somme_sinus(parameters_x, T, fe)

# Affichage du signal
plt.figure()
plt.plot(t, x)
plt.title("Signal x(t)")
plt.tight_layout()
plt.show()

# Transformée de Fourier de x(t) et affichage du spectre d'amplitude
X = np.fft.fftshift(np.fft.fft(x))
f = np.fft.fftshift(np.fft.fftfreq(len(x),1/fe))
amp = np.abs(X)*1/len(x)

plt.stem(f, amp, markerfmt=' ')
plt.title(f"Spectre d'amplitude de X(f), fe = {fe} Hz")
plt.show()

# Filtre passe-bas idéal appliqué à x(t), affichage du spectre d'amplitude modifié et détection des pics (fréquence, amplitude)
X[np.abs(f)>fc]=0

amp = np.abs(X)*1/len(x)
plt.stem(f, amp, markerfmt=' ')
plt.title(f"Spectre d'amplitude de X(f), fe = {fe} Hz, fc = {fc} Hz")
plt.show()

peaks, properties = find_peaks(amp, height=np.max(amp)*0.01)
print(f'Pics détectés: ')
for peak, height in zip(peaks, properties['peak_heights']):
    print(f"• {f[peak]} Hz - Amplitude {height}")

Pics détectés: 
• -400.0 Hz - Amplitude 1.000000000000003
• -300.0 Hz - Amplitude 0.5000000000000021
• 300.0 Hz - Amplitude 0.5000000000000021
• 400.0 Hz - Amplitude 1.000000000000003

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

parameters_a = [[5, 2, -np.pi/2],
                [20, 1, -np.pi/3]]

t,a = somme_sinus(parameters_a, 5, 44100)

# Affichage du signal
plt.figure(figsize=(16,4))
plt.plot(t, a)
plt.title("Signal a(t)")
plt.tight_layout()
plt.show()

# A COMPLETER
# Normalisation de a(t) et affichage

amplitude_max_a = np.max(np.abs(a))
a_norm = a/amplitude_max_a

# Affichage du signal normalisé
plt.figure(figsize=(16,4))
plt.plot(t, a_norm)
plt.title("Signal a(t) normalisé")
plt.tight_layout()
plt.show()

# A COMPLETER
# Quantification du signal normalisé en 8 bits et affichage du signal quantifié
a_quant_8_bit = np.round(a_norm*127).astype(np.int8)

# Affichage du signal quantifié
plt.figure(figsize=(16,4))
plt.plot(t, a_quant_8_bit)
plt.title("Signal a(t) quantifié 8 bits")
plt.tight_layout()
plt.show()

# A COMPLETER
# Modification du signal quantifié pour le comparer au signal d'origine
a_prime = (a_quant_8_bit/127)*amplitude_max_a

# Affichage du signal d'origine et du signal quantifié
plt.figure(figsize=(16,4))
plt.plot(t, a, label='Signal original')
plt.plot(t, a_prime, label='Signal quantifié 8 bits signé')
plt.title("Quantification 8 bit de a(t)")
plt.legend()
plt.tight_layout()
plt.show()

# A COMPLETER
# Calcul du nombre de valeurs "uniques" des signaux
print(f"Nombres de valeurs d'amplitudes uniques de a(t) : {len(np.unique(a))}")
print(f"Nombres de valeurs d'amplitudes uniques de a(t) quantifié 8 bits : {len(np.unique(a_prime))}")

Nombres de valeurs d'amplitudes uniques de a(t) : 206400
Nombres de valeurs d'amplitudes uniques de a(t) quantifié 8 bits : 240

# A COMPLETER
# Calcul du rapport signal sur bruit entre les deux signaux
def rapport_sb(signal,signal_compresse):
    return 10 * np.log10(np.sum(signal.astype(float)**2) / np.sum((signal.astype(float) - signal_compresse.astype(float))**2))

print(f"Rapport S/B de la quantification 8 bits : {rapport_sb(a,a_prime)} dB")

Rapport S/B de la quantification 8 bits : 47.38869816491939 dB

# A COMPLETER
# Fonction de quantification N bits
def quantify(signal, N, amp_max):
    """
    Quantifie un signal d'entrée en N bits
    signal : tableau 1D NumPy, signal à quantifier
    N : entier, nombre de bits pour quantification
    amp_max : float, amplitude maximale possible du signal d'entrée
    """

    # Normalisation du signal d'entrée
    signal_norm = np.clip(signal/amp_max,-1,1) # Avec np.clip, on replace les valeurs entre -1 et 1 en définissant le min à -1 et le max à 1

    # Quantification n bits signée
    # On cast en int64 (entier codé sur 64 bits) par défaut, mais on pourrait trouver un type de données moins gourmand en fonction de N
    signal_quant = np.round(signal_norm * 2**(N-1)-1).astype(np.int64)

    return signal_quant

# Quantification de a(t) en 4 bits
a_quant_4_bit = quantify(a,4,amplitude_max_a)

# Affichage du signal a(t) quantifié 4-bits
plt.figure(figsize=(16,4))
plt.plot(t, a_quant_4_bit)
plt.title("Signal a(t) quantifié 4 bits")
plt.tight_layout()
plt.show()

# A COMPLETER
# Fonction de déquantification du signal
def unquantify(signal_quant, N, amp_max):
    """
    Quantifie un signal d'entrée en N bits
    signal_quant : tableau 1D NumPy, signal à quantifier
    N : entier, nombre de bits pour quantification
    amp_max : float, amplitude maximale possible du signal d'entrée
    """
    signal_dequant = (signal_quant/(2**(N-1)-1))*amp_max
    return signal_dequant

# Dé-quantification de a(t) 4 bits
a_dequant_4_bit = unquantify(a_quant_4_bit,4,amplitude_max_a)

# Affichage du signal a(t) déquantifié 4-bits et du signal d'origine
plt.figure(figsize=(16,4))
plt.plot(t, a, label='Signal original')
plt.plot(t, a_dequant_4_bit, label='Signal quantifié 4 bits signé')
plt.title("Quantification 4 bit de a(t)")
plt.tight_layout()
plt.show()

# A COMPLETER
# Quantification N bits de a(t) et affichage du rapport signal sur bruit en fonction de N

amplitude_max_a = np.max(np.abs(a))
N_values = list(range(2,17))
RSB_values = []
for N in N_values:
    a_quant_N_bits = quantify(a, N, amplitude_max_a)
    a_dequant_N_bits = unquantify(a_quant_N_bits, N, amplitude_max_a)
    RSB_values.append(rapport_sb(a,a_dequant_N_bits))

plt.figure()
plt.plot(N_values,RSB_values)
plt.title("Rapport S/B du signal audio quantifié N bits par rapport à l'original")
plt.xlabel("N")
plt.ylabel("RSB")
plt.show()

# A COMPLETER
# Lecture de l'audio, et affichage de la fréquence d'échantillonnage et de la profondeur de quantification
...

Ellipsis

# A COMPLETER
# Ré-échantillonnage et quantification de l'audio HQ, et enregistrement dans différents fichiers WAV
...

Ellipsis

Traitement du Signal - TP6 : Voyage dans le monde numérique !¶

Exercice 1 : Echantillonnage¶

Exercice 2 : Et 1, et 2, et 3 sinus !¶

Exercice 3 : Quantification¶

Exercice 4 : Ca donne quoi sur l'audio ?¶