% Segnali e Sistemi a.a. 2021-22
% Es1 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);
t=[0:N-1]*Tc; %asse dei tempi
figure(1)
plot(t,ecg,'b-')
title('ECG')
xlabel('t(s)')
ylabel('microV')
grid on
axis([0 t(end) 1400 3000])
% zoom
figure(2)
plot(t,ecg,'b-')
title('ECG')
xlabel('t(s)')
ylabel('microV')
grid on
axis([0 5 1400 3000])

%sottrazione della componente continua
s = ecg-mean(ecg);

% calcolo della trasformata di Fourier
S = fftshift(Tc*fft(s));
omega0 = 2*pi/(Tc*N); % passo di campionamento (in frequenza)
%asse delle pulsazioni
omega=(-round((N-1)/2):round(N/2)-1)*omega0; 

% plot del modulo della trasformata
figure
semilogy(omega,abs(S))
xlabel('omega [rad/s]');
axis([0 omega(end) 1e-1 2e3])
grid on
title('Trasformata di Fourier del segnale ECG');

%%%  zoom in scala lineare
figure
plot(omega,abs(S))
xlabel('omega [rad/s]');
axis([0 8*pi 1e0 2e3])
grid on
title('Zoom della Trasformata del Segnale ECG');

%identificazione del primo massimo
index_range=find(omega>=pi & omega<2*pi);
omega_range=omega(index_range);
S_range=S(index_range);

[maxval,pos] = max(abs(S_range));
disp(['Stima periodicità : omega=' num2str(omega_range(pos)) ', Tp=' num2str(2*pi/omega_range(pos))])

 hold on
 plot(omega_range(pos), abs(S_range(pos)),'xr','LineWidth',2)
 hold off
%figure(2)
%hold on
%plot(t,mean(ecg)+maxval*cos(2*pi*f_range(pos)*(t-0.712)),'g--')
%hold off