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

% carica i dati ecg
load data_ecg.mat 
fc = 125; % Hz
wc = 2*pi*fc; % rad/s 
Tc = 1/fc; % s

% sottrazione della componente continua
% serve solo a non avere un picco in omega=0
s = ecg - mean(ecg);
N = length(s);
t = [0:N-1]*Tc; % asse dei tempi

%% mostro il segnale nel tempo (con e senza zoom)
figure(1)
subplot(2,2,1)
plot(t,s,'b-')
title('ECG')
xlabel('t [s]')
ylabel('s(t) [microV]')
grid on
axis([0 t(end) ylim])
subplot(2,2,3)
plot(t,s,'b-')
title('Zoom di ECG')
xlabel('t [s]')
ylabel('s(t) [microV]')
grid on
axis([0 5 ylim])

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

% mostro la trasformata (con e senza zoom)
subplot(2,2,2)
semilogy(omega,abs(S))
xlabel('\omega [rad/s]');
ylabel('S(\omega) [microV/Hz]')
axis([0 omega(end) 1e-1 2e3])
grid on
title('Trasformata di Fourier di ECG');
subplot(2,2,4)
plot(omega,abs(S))
xlabel('\omega [rad/s]');
ylabel('S(\omega) [microV/Hz]')
axis([0 8*pi 1e0 2e3])
grid on
title('Zoom della trasformata di ECG');

%% cerco il massimo
% cerco la sua posizione nel vettore
pos = find((abs(S)==max(abs(S))) & (omega>0))
% uso 'pos' per mostrare nel plot dove sta il massimo 
hold on
plot(omega(pos), abs(S(pos)),'xr','LineWidth',2)
hold off
% riconosco che sta nella second armonica, quindi
% omega(pos) = 2 om0 = 4 pi / Tp
% uso la formula per trovare Tp
disp(['Stima periodicità : Tp=' num2str(4*pi/omega(pos))])