clear all; close all; clc;
 
 
 
%% Traitement du signal | Matlab #3: Synthèse du filtre Butterworth
%% Cours & Projets: www.Electronique-Mixte.fr


% 1. Savoir l'importance du filtre 
% 2. Savoir tracer la fonction du transfert (FT) (Gain & Phase) 
%    d'un filtre quelconque
% 3. Exemple du filtrage d'un signal ECG bruité 
% 4. Limitations du filtre 
% 5. Etc.





%% Paramètres du signal ECG 
 
    
T=1; %% Période normalisée Toujours=1!!!!!!!!!!!! 
n=500;
ecg_dc=0.14; 
    
p_wave   =[0.1 0.1*T 0.4*T];    % 3
    % a_pwav=p_wave(1);   % Amplitude
    % d_pwav=p_wave(2);   % Duration
    % t_pwav=p_wave(3);   % p-r interval
 
q_wave   =[0.2 0.1*T ];         % -  
qrs_wave =[1 0.1*T];
s_wave   =[0.3 0.2*T ];         % -
t_wave   =[0.2 0.2*T 0.2*T];    % 3
u_wave   =[0.1 0.1*T];
 
% Génération du signal 
ECG_s = ecg(n, p_wave, q_wave, qrs_wave, s_wave, t_wave, u_wave);
ECG_s=ECG_s-ecg_dc; % Suppression de la composante DC

 
 
%% Périodisation d'un signal 
 
% Génération d'une distribution peine de Dirac
f0=1; t0=1/f0; n0=10; N=n0*n; 
t=linspace(0,n0*t0/2,N);
A=1; s_n=A*sin(2*pi*f0*t);

s_b=sign(s_n); 
s_diff=abs(conv(s_b,[-1 1], 'same'));
s_diff=s_diff/2;

% bords
s_diff(1)=0;
s_diff(end)=0;


%% Périodisation 

ecg_p=conv(s_diff,ECG_s, 'same');
ecg_p=ecg_p+ecg_dc; 

%% Filtrage de Butterworth 


% Params du filtre 
ts=t(2)-t(1); fs=1/ts;  % Fréquence Fs 

fc=(fs/2)/4;            % Fréquence de coupure

fn=fc/(fs/2);           % Fréquence normalisée (ou pulsation)

n_fil=12;                % Ordre du filtre 

% Génération des coefs 
[b, a]=butter(n_fil,fn, 'low'); % 'high', 'stop', 'low', 'bandpass'

% Affichage de la FT du filtre 
figure(1); freqz(b,a); % Ou bien fvtool(b,a);




%% Ajout du bruit 

% Bruit
sig=0.5e-1; 
b_n=sig*randn(1,N);

% Ajout du bruit 
s_b_n=b_n+ecg_p; 

% Filtrage 
s_fil=filter(b,a,s_b_n);


%% Affichage 

figure(2)
subplot(211); plot(t, ecg_p, 'linewidth',2); grid on;
xlabel('Temps(s)'); ylabel('Amplitude(V)'); legend('Original'); 
xlim([t(1) t(end)]); 

subplot(212); plot(t, s_b_n, 'linewidth',2); grid on;
xlabel('Temps(s)'); ylabel('Amplitude(V)'); legend('Bruité'); 
xlim([t(1) t(end)]); 




figure(3)
subplot(311); plot(t, ecg_p, 'linewidth',2); grid on;
xlabel('Temps(s)'); ylabel('Amplitude(V)'); legend('Original'); 
xlim([t(1) t(end)]); 

subplot(312); plot(t, s_b_n, 'linewidth',2); grid on; hold on; 
xlabel('Temps(s)'); ylabel('Amplitude(V)'); legend('Bruité'); 
xlim([t(1) t(end)]); 


subplot(313); plot(t, s_fil, 'r', 'linewidth',2); grid on; hold off;
xlabel('Temps(s)'); ylabel('Amplitude(V)'); legend('Filtré'); 
xlim([t(1) t(end)]);  %title(num2str(i));