% Segnali e Sistemi a.a. 2021-22
% Es2 Lab 6, ECG
clear all
clc
close all

% CARICA I DATI 
load data_ecg.mat %carica i dati ecg ed ecg_dist
fc=125; %Hz 
wc=2*pi*fc; %rad/s 
Tc=1/fc; %s
N=length(ecg_dist);
t=[0:N-1]*Tc; %asse dei tempi

figure
subplot(2,1,1)
plot(t,ecg_dist,'b-')
xlabel('t(s)')
ylabel('microV')
grid on
axis([0 t(end) 1400 3000])
title('ECG disturbato')

subplot(2,1,2)
plot(t,ecg,'r-')
title('ECG')
xlabel('t(s)')
ylabel('microV')
grid on
axis([0 t(end) 1400 3000])

sdist=ecg_dist-mean(ecg_dist);
Tf_sdist=fftshift(Tc*fft(sdist));
omega0 = 2*pi/(Tc*N); % passo di campionamento (in frequenza)
%asse delle pulsazioni
omega=(-round((N-1)/2):round(N/2)-1)*omega0; 

% Crea un filtro PA idale con cutoff 0.4 pi rad/s
wcut=0.4*pi;
ind=find(omega<=wcut & omega>=-wcut);

Tf_sfilt=Tf_sdist;
Tf_sfilt(ind)=zeros(length(Tf_sfilt(ind)),1);

ecg_filt=1/Tc*ifft(ifftshift(Tf_sfilt))+mean(ecg_dist);


figure
subplot(2,1,1)
plot(t,ecg_filt,'b-')
xlabel('t(s)')
ylabel('microV')
grid on
axis([0 t(end) 1400 3000])

title('ECG filtrato')
subplot(2,1,2)
plot(t,ecg,'r-')
title('ECG')
xlabel('t(s)')
ylabel('microV')
grid on
axis([0 t(end) 1400 3000])