# -*- coding: utf-8 -*-
"""
Created on Tue Mar  4 21:22:43 2025

@author: AKourgli
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import upfirdn, welch

# Paramètres
sps = 128                             # Échantillons par symbole
num_symbols = 5000
alphas = [0, 1]
span = 10                             # Durée du filtre
snr_db = 30

# Génération des symboles BPSK
symbols = np.random.choice([-1, 1], num_symbols)
samples_per_eye = 2 * sps
t_eye = (np.arange(2*sps)) / sps

# Upsampling
upsampled = np.zeros(num_symbols * sps)
upsampled[::sps] = symbols

# Configuration des graphiques (3 colonnes : Œil, DSP, Réponse impulsionnelle)
fig, axes = plt.subplots(len(alphas), 4, figsize=(15, 2 * len(alphas)))

# Fonction PSD théorique
def raised_cosine_psd(f, alpha):
    f_abs = np.abs(f)
    psd = np.zeros_like(f)
    if alpha == 0:
        psd[f_abs <= 0.5] = 1
    else:
        mask1 = f_abs <= (1 - alpha)/2
        mask2 = (f_abs > (1 - alpha)/2) & (f_abs <= (1 + alpha)/2)
        psd[mask1] = 1
        psd[mask2] = 0.5 * (1 + np.cos( (np.pi / alpha) * (f_abs[mask2] - (1 - alpha)/2 )))
    return psd

for idx, alpha in enumerate(alphas):
    # Conception du filtre RRC
    t = np.linspace(-span/2, span/2, span * sps + 1)
    h = np.zeros_like(t)
    
    for i, tt in enumerate(t):
        tt += 1e-9  # Éviter tt=0
        if alpha == 0:
            h[i] = np.sinc(tt)
        elif alpha == 1:
            if np.isclose(abs(tt), 0.25, atol=1e-6):
                h[i] = (1 / np.sqrt(2)) * ((1 + 2/np.pi) * np.sin(np.pi/4) + (1 - 2/np.pi) * np.cos(np.pi/4))
            elif tt == 0:
                h[i] = 1 + 4 / np.pi
            else:
                num = 4 * tt * np.cos(2 * np.pi * tt)  # Numérateur corrigé
                den = np.pi * tt * (1 - (4 * tt)**2)
                h[i] = num / den
        else:
            if np.isclose(abs(tt), 1/(4 * alpha), atol=1e-6):
                h[i] = (alpha / np.sqrt(2)) * ((1 + 2/np.pi) * np.sin(np.pi/(4*alpha)) + (1 - 2/np.pi) * np.cos(np.pi/(4*alpha)))
            elif tt == 0:
                h[i] = 1 - alpha + 4 * alpha / np.pi
            else:
                num = np.sin(np.pi * tt * (1 - alpha)) + 4 * alpha * tt * np.cos(np.pi * tt * (1 + alpha))
                den = np.pi * tt * (1 - (4 * alpha * tt)**2)
                h[i] = num / den
    
    h /= np.linalg.norm(h)  # Normalisation

    # Filtrage
    filtered = upfirdn(h, upsampled, up=1, down=1)
    filtered = filtered[span*sps//2 : -span*sps//2]

    # # Tracé du diagramme de l'œil
    offset = span * sps // 2
    for i in range(100):
        start = offset + i * sps
        end = start + 2 * sps
        if end < len(filtered):
            axes[idx, 3].plot(t_eye,filtered[start:end], 'navy', alpha=0.1)
    axes[idx, 3].set_title(f"Œil (α={alpha})")
    axes[idx, 3].grid(True)

    # Tracé avec bruit (optionnel)
    signal_power = np.mean(filtered**2)
    noise_power = signal_power / (10 ** (snr_db / 10))
    noise = np.random.normal(0, np.sqrt(noise_power), len(filtered))
    filtered_noisy = filtered + noise

    for i in range(100):
        start = offset + i * sps
        end = start + 2 * sps
        if end < len(filtered_noisy):
            axes[idx, 0].plot(t_eye,filtered_noisy[start:end], 'navy', alpha=0.1)
    axes[idx, 0].set_title(f"RRC (α={alpha}) - SNR={snr_db}dB")
    axes[idx, 0].set_xlabel("Temps (T)")
    axes[idx, 0].grid()
    
    # Tracé de la DSP
    f_theo = np.linspace(-1, 1, 1000)
    psd_theo = raised_cosine_psd(f_theo, alpha)
    psd_theo_db = 10 * np.log10(psd_theo + 1e-9)
    axes[idx, 2].plot(f_theo, psd_theo_db, 'k--', label='Théorique')
    axes[idx, 2].set_title(f"DSP (α={alpha})")
    axes[idx, 2].set_xlim(-1.2, 1.2)
    axes[idx, 2].set_ylim(-60, 5)
    axes[idx, 2].grid(True)
    axes[idx, 2].legend()

    # Tracé de la réponse impulsionnelle h
    t_symbols = t  # Temps en périodes symbole
    axes[idx, 1].plot(t_symbols, h, 'g')
    axes[idx, 1].set_title(f"Réponse impulsionnelle (α={alpha})")
    axes[idx, 1].set_xlabel("Temps (T)")
    axes[idx, 1].set_ylabel("Amplitude")
    axes[idx, 1].grid(True)
    

plt.tight_layout()
plt.show()