// ================================
// Paramètres utilisateur
// ================================
directory = "C:\Users\Fourcade\Documents\Logiciel PC\Audio\HornResp\Export\";
filename = directory+"50Horn.txt";  //  fichier d'entrée
fc       = 55;                         // fréquence de coupure en Hz
N        = 2;                          // ordre Butterworth utilisé (LR sera d'ordre 2N)

// Type de filtre : "HP" (passe-haut) ou "LP" (passe-bas)
filterType = "HP";

// Application d'un brickwall à fc
brickwall = 0 ;

// ================================
// Lecture du fichier 
// ================================
data = fscanfMat(filename);

f = data(:, 1);   // fréquences en Hz
X = data(:, 7);   // élongation non filtrée

// Gestion éventuelle du point f = 0 pour le passe-haut
f_eps = f;
eps_val = 1.0e-9;
f_eps(f_eps == 0) = eps_val;

// ================================
// Fonctions utilitaires : module du filtre LR
// ================================

// Butterworth passe-bas (amplitude) d'ordre N
function Hb = butter_lp_mag(f, fc, N)
    // f, fc > 0 (Hz)
    Hb = 1 ./ sqrt(1 + (f ./ fc).^(2*N));
endfunction

// Butterworth passe-haut (amplitude) d'ordre N
function Hb = butter_hp_mag(f, fc, N)
    // f, fc > 0 (Hz)
    Hb = 1 ./ sqrt(1 + (fc ./ f).^(2*N));
endfunction

// Linkwitz–Riley passe-bas (amplitude) d'ordre 2N
function Hlr = lr_lp_mag(f, fc, N)
    Hb = butter_lp_mag(f, fc, N);
    Hlr = Hb.^2; // cascade de deux Butterworth
endfunction

// Linkwitz–Riley passe-haut (amplitude) d'ordre 2N
function Hlr = lr_hp_mag(f, fc, N)
    Hb = butter_hp_mag(f, fc, N);
    Hlr = Hb.^2;
endfunction

// ================================
// Calcul de la réponse du filtre
// ================================
select filterType
case "HP" then
    H = lr_hp_mag(f_eps, fc, N);
case "LP" then
    H = lr_lp_mag(f, fc, N);
else
    error("filterType doit valoir ""HP"" ou ""LP"".");
end

// ================================
// Application du filtre sur l'élongation
// (domaine fréquentiel : Xfiltrée(f) = X(f) * |H(f)| )
// ================================
X_filt = X .* H;

// Application du brickwall
if brickwall==1 then
    // Recherche indice fc
    [ecart_min, indice] = min(abs(f - fc)) ;
    // Mise à zero valeur
    X_filt(1:indice) = 0 ;
end 
// ================================
// Tracés
// ================================
// Color : noir 1, bleu=2, green=3, bleu clair=4, red=5,  purple=6, 
scf(0);

// Courbe d'élongation avant / après
//subplot(2,1,1);
plot2d(f, X);
plot2d(f, X_filt,style=5);
h=gca();
h.grid=[1,1]
h.log_flags="lnn"; // x:log,y=log,z=lin
h.sub_ticks=[0,0];

// X ticks
ticks = h.x_ticks
ticks.labels = ["10";"20";"40";"80";"120";"180"];
ticks.locations = [10;20;40;80;120;180];
h.x_ticks = ticks;
xlabel("Fréquence (Hz)");
ylabel("Élongation (m ou unité)");
title("Élongation de la membrane - Avant / Après filtre LR");
//legends(['LR2';'LR4';'LR8';'LR16'],[2,3,4,5],opt="ur")
legends(['Mixte LR2+Brickwall'],[5],opt="ur")
xgrid();

