# -*- coding: utf-8 -*-
"""
Created on Mon Mar 10 16:54:32 2025

@author: AKourgli
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import periodogram, butter, filtfilt, freqz, lfilter

# Paramètres
ts = 0.01        # Pas d'échantillonnage (en s)
Tb = 1           # Temps d'un bit (en s)
Nb = 100         # Nombre de bits
fs = 1 / ts      # Fréquence d'échantillonnage

# Génération du signal NRZ polaire
bits = np.random.choice([-1, 1], Nb)
signal = np.repeat(bits, int(Tb / ts))
t = np.arange(len(signal)) * ts

# Calcul de la DSP du signal original avec le periodogram
f, Pxx = periodogram(signal, fs)

# Filtrage : Passage par le canal
fc = fs / 90  # Fréquence de coupure à fs/4
Ncoef = 6
b, a = butter(Ncoef, fc / (fs / 2), btype='low')
signal_filtre = lfilter(b, a, signal)

# Ajout de bruit
SNR_dB = 20
puissance_signal = np.mean(signal**2)
puissance_bruit = puissance_signal / (10**(SNR_dB/10))
bruit = np.random.normal(0, np.sqrt(puissance_bruit), len(signal))
signal_filtre_bruite = signal_filtre + bruit
#signal_bruite = signal_filtre

# Calcul de la DSP du signal filtré avec le periodogram
f_filt, Pxx_filt = periodogram(signal_filtre, fs)

# Calcul de la réponse en fréquence (DSP) du filtre
w, h = freqz(b, a, worN=1024, fs=fs)
dsp_filter = np.abs(h)**2

#Calcul de la réponse impulsionnelle du filtre
N_imp = 200  # Nombre d'échantillons pour l'impulsion
impulse = np.zeros(N_imp)
impulse[0] = 1
imp_resp = lfilter(b, a, impulse)

# Affichage des résultats
plt.figure()

Subs = int(Tb/ts*20)
# Signal original
plt.subplot(321)
plt.plot(t[:Subs], signal[:Subs], label="Signal NRZ polaire")
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal NRZ polaire")
plt.legend()

# DSP du signal original
plt.subplot(322)
plt.semilogy(f, Pxx, label="DSP originale")
plt.xlabel("Fréquence (Hz)")
plt.ylabel("DSP")
plt.title("DSP du signal original")
plt.ylim(1e-9, 10)
plt.legend()

# DSP du filtre (réponse en fréquence du canal)
plt.subplot(324)
plt.semilogy(w, dsp_filter, label="DSP du canal", color="red")
plt.xlabel("Fréquence (Hz)")
plt.ylabel("Gain")
plt.ylim(1e-9, 10)
plt.title("Densité spectrale de puissance du canal")
plt.legend()

# Signal filtré
plt.subplot(323)
plt.plot(t[:N_imp],imp_resp, label="h canal", color="orange")
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal après passage par le canal en bande limitée sans bruit")
plt.legend()

# Réponse impulsionnelle du filtre
plt.subplot(325)
plt.plot(t[:Subs], signal[:Subs], label="Réponse impulsionnelle du canal")
plt.plot(t[:Subs],signal_filtre_bruite[:Subs], label="Signal reçu")
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal avant et après passage par le canal en bande limitée avec bruit ")
plt.legend()

# DSP du signal filtré
plt.subplot(326)
plt.semilogy(f_filt, Pxx_filt, label="DSP du signal en sortie du canal")
plt.xlabel("Fréquence (Hz)")
plt.ylabel("DSP")
plt.title("DSP du signal filtré")
plt.ylim(1e-9, 10)
plt.legend()

plt.tight_layout()
plt.show()

# Fonction pour tracer le diagramme de l'œil
def plot_eye(sig, ts, Tb, title):
    # Pour une meilleure visualisation, on utilise des segments de 2 périodes de bit
    samples_per_bit = int(Tb / ts)
    samples_per_segment = 2 * samples_per_bit
    n_segments = len(sig) // samples_per_segment
    segments = sig[:n_segments * samples_per_segment].reshape(n_segments, samples_per_segment)
    t_eye = np.linspace(0, 2*Tb, samples_per_segment)
    for seg in segments:
        plt.plot(t_eye, seg, 'navy', alpha=0.1)
    plt.xlabel("Temps (s)")
    plt.ylabel("Amplitude")
    plt.title(title)

# Affichage des diagrammes de l'œil dans une nouvelle figure
plt.figure()

plt.subplot(2, 1, 1)
plot_eye(signal, ts, Tb, "Diagramme de l'œil - Signal émis")

plt.subplot(2, 1, 2)
plot_eye(signal_filtre_bruite, ts, Tb, "Diagramme de l'œil - Signal reçu")

plt.tight_layout()
plt.show()