%% Parte 1
clear all;
close all;
clc;
% Scelta dei parametri
f0=15;
Fs=400;
A=2;
Tw=15*1/f0;
F=1/Tw;
N=Fs*Tw;
% Creazione dei vettori tempo e segnale
t=[0:N-1]*1/Fs;
x=A*sin(2*pi*f0*t);
figure
plot(t,x);
xlabel('Time (s)');
ylabel('Amplitude (V)');

% Analisi dello spettro
X=fft(x); % calcola i coefficienti complessi di Fourier
X=abs(X)/N; % ne calcoliamo il modulo normalizzato per la somma dei pesi 
            % della finestra rettangolare che è quella di default
X=2*X(1:N/2);  % consideriamo per simmetria solo le prime N/2 componenti, 
               % ma moltiplichiamo per 2 per tenere conto del contributo
               % delle componenti immagine
f=[0:N/2-1]*F; % creiamo un vettore delle frequenze adeguato, in cui i 
               % campioni saranno distanziati di F, di lunghezza coerente
figure    %varie opzioni per visualizzare lo stesso spettro
subplot (3,1,1)
stem(f,X);
xlabel('Frequency (Hz)');
ylabel('Amplitude (V)');
subplot (3,1,2)
plot(f,X);
xlabel('Frequency (Hz)');
ylabel('Amplitude (V)');
subplot (3,1,3)
plot(f,20*log10(X));
xlabel('Frequency (Hz)');
ylabel('Amplitude (dBV)');
%% parte 2

f0=13*F; %valore per cui esista k0 intero per cui f0=k0*F
f0bis=12.6*F;% valore per cui non esista k0 intero per cui f0=k0*F
% N.B. Tw e quindi F resteranno quelli definiti nella parte 1. 
f=[0:N/2-1]*F; % creiamo un vettore delle frequenze adeguato, in cui i 
               % campini saranno distanziati di F, di lunghezza coerente
%creazione dei due segnali
x1=A*sin(2*pi*f0*t);
x1bis=A*sin(2*pi*f0bis*t);

% Per velocizzare il tutto il consiglio è raccogliere le linee di codice
% scritte sopra in una function Matlab da salvare in un file separato
% oppure aggiungere in coda
X1=fft_agile(x1);
X1bis=fft_agile(x1bis);

figure
subplot(2,1,1)
hold on
plot(f,20*log10(X1));
plot(f,20*log10(X1bis),'r');
legend('f0=13F','f0=12.6F');
xlabel('Frequency (Hz)');
ylabel('Amplitude (dBV)');
subplot(2,1,2)
hold on
stem(f,X1);
stem(f,X1bis,'r');
xlabel('Frequency (Hz)');
ylabel('Amplitude (V)');
legend('f0=13F','f0=12.6F');

%% parte 3
w=hann(N);% creare il vettore finestra di hanning 
w1=flattopwin(N); % creare il vettore finestra Flat Top

% Pesare il segnale in cui compariva distorsione per le due finestre 
x1bisw=x1bis.*w';
x1bisw1=x1bis.*w1';

% Confrontare nel dominio del tempo
figure
hold on
plot(t,x1bis);
plot(t,x1bisw,'r');
plot(t,x1bisw1,'g');
xlabel('Time (s)');
ylabel('Amplitude (V)');
legend('Rectangular','Hanning','Flat Top');

% Calcolare la fft dei due segnali finestrati, attenzione a pesare non più
% per N ma per la somma dei campioni; per velocizzare si può creare una
% function fft_finestra in cui si modifica praticamente solo il passaggio
% relativo alla normalizzazione 
X1bisw=fft_finestra(x1bis,w);
X1bisw1=fft_finestra(x1bis,w1);

figure
subplot(2,1,1)
hold on
plot(f,20*log10(X1bis));
plot(f,20*log10(X1bisw),'r');
plot(f,20*log10(X1bisw1),'g');
legend('Rectangular','Hanning','Flat Top');
xlabel('Frequency (Hz)');
ylabel('Amplitude (dBV)');
subplot(2,1,2)
hold on
stem(f,X1bis);
stem(f,X1bisw,'r');
stem(f,X1bisw1,'g');
xlabel('Frequency (Hz)');
ylabel('Amplitude (V)');
legend('Rectangular','Hanning','Flat Top');

%% parte 4
% definisco i parametri (questi sono solo a titolo di pure esempio)
f1=5; 
f2=34; 
f3=127;
Fs=600;
Tw=20.3*1/f1;
F=1/Tw;
N=round(Fs*Tw);
%costruisco il vettore tempo
t=[0:N-1]*1/Fs;
%costruisco i vettori delle varie componenti
y1=sin(2*pi*f1*t);
y2=sin(2*pi*f2*t);
y3=sin(2*pi*f3*t);
%costruisco il segnale composto
signal=y1+y2+y3;

%costruisco i vettori finestre della stessa lunghezza del segnale
w=hann(N);
w1=flattopwin(N);

figure
hold on
plot(t,signal);
plot(t,signal.*w');
plot(t,signal.*w1');
xlabel('Time (s)');
ylabel('Amplitude (V)');
legend('Rectangular','Hanning','Flat Top');

% Calcolo gli spettri senza e con finestra usando le funzioni definite
% prima
Xsig=fft_agile(signal);
Xsigw=fft_finestra(signal,w);
Xsigw1=fft_finestra(signal,w1);

f=[0:N/2-1]*F; % creiamo un vettore delle frequenze adeguato, in cui i 
               % campini saranno distanziati di F, di lunghezza coerente
figure
hold on
plot(f,20*log10(Xsig));
plot(f,20*log10(Xsigw));
plot(f,20*log10(Xsigw1));
legend('Rectangular','Hanning','Flat Top');
xlabel('Frequency (Hz)');
ylabel('Amplitude (dBV)');
figure
hold on
stem(f,Xsig);
stem(f,Xsigw);
stem(f,Xsigw1);
xlabel('Frequency (Hz)');
ylabel('Amplitude (V)');
legend('Rectangular','Hanning','Flat Top');

%% parte 5
% close all;
load('ECG.mat');
load('PPG.mat');
Fs=125;
N=1500;
t=[0:N-1]*1/Fs;
figure
subplot(2,1,1)
plot(t,ECG,'k');
hold on
ECG_corr=ECG-mean(ECG);
plot(t,ECG_corr,'r');
xlabel('Time (s)');
ylabel('Amplitude (mV)');
legend('ECG','ECG corr');
subplot(2,1,2)
plot(t,PPG,'k');
hold on
PPG_corr=PPG-mean(PPG);
plot(t,PPG_corr,'b');
xlabel('Time (s)');
ylabel('Amplitude (mV)');
legend('PPG','PPG corr');

%calcolo dello spettro senza finestra
XECG=fft_agile(ECG_corr);
XPPG=fft_agile(PPG_corr);
%calcolo dello spettro con finestra di Hanning
w=hann(N);
XECGw=fft_finestra(ECG_corr,w);
XPPGw=fft_finestra(PPG_corr,w);

Tw=N*1/Fs;
F=1/Tw;
f=[0:N/2-1]*F;
figure
% subplot(2,1,1)
% plot(f,XECG,'r');
% hold on
% plot(f,XPPG,'b');
% legend('Spettro ECG','Spettro PPG');
% xlabel('Frequency (Hz)');
% ylabel('Amplitude (mV)');
% subplot(2,1,2)
plot(f,XECGw,'r');
hold on
plot(f,XPPGw,'b');
legend('Spettro ECG Hanning windowed','Spettro PPG Hanning windowed');
xlabel('Frequency (Hz)');
ylabel('Amplitude (mV)');
%% METODO 2
[pks,locs] = findpeaks(XECGw);
[pks1,locs1] = findpeaks(XPPGw);

% for PPG signal
[~,indice_max]=max(pks1);
indice_hr_ppg=locs1(indice_max);
HR_ppg_met2=(indice_hr_ppg-1)*F;

% for ECG signal
f_locs=(locs-1)*F;
indici_range=find(f_locs>0.8&f_locs<3.5);
[~,indice_max]=max(pks(indici_range));
indice_hr_ecg=locs(indici_range(indice_max));
HR_ecg_met2=(indice_hr_ecg-1)*F;

figure
plot(f,XECGw,'r');
hold on
plot(f,XPPGw,'b');
xlabel('Frequency (Hz)');
ylabel('Amplitude (mV)');

plot(f(locs),pks,'r*','LineWidth',2);
plot(f(locs1),pks1,'b*','LineWidth',2);
legend('Spettro ECG Hanning windowed','Spettro PPG Hanning windowed','Pks ECG','Pks PPG');
%% METODO 3
th=8;
[pks,locs] = findpeaks(XECGw,'MinPeakHeight',th);
[pks1,locs1] = findpeaks(XPPGw,'MinPeakHeight', th);

% for PPG signal
[~,indice_max]=max(pks1);
indice_hr_ppg=locs1(indice_max);
HR_ppg_met3=(indice_hr_ppg-1)*F;

% for ECG signal
f_locs=(locs-1)*F;
indici_range=find(f_locs>0.8&f_locs<3.5);
[~,indice_max]=max(pks(indici_range));
indice_hr_ecg=locs(indici_range(indice_max));
HR_ecg_met3=(indice_hr_ecg-1)*F;

figure
plot(f,XECGw,'r');
hold on
plot(f,XPPGw,'b');

xlabel('Frequency (Hz)');
ylabel('Amplitude (mV)');
plot(f(locs),pks,'r*','LineWidth',2);
plot(f(locs1),pks1,'b*','LineWidth',2);
legend('Spettro ECG Hanning windowed','Spettro PPG Hanning windowed','Pks ECG','Pks PPG');

%% funzioni

function X=fft_agile(x)
X=fft(x); % calcola i coefficienti complessi di Fourier
N=length(x);
X=abs(X)/N; % ne calcoliamo il modulo normalizzato per la somma dei pesi 
            % della finestra rettangolare che è quella di default
X=2*X(1:N/2);  % consideriamo per simmetria solo le prime N/2 componenti, 
               % ma moltiplichiamo per 2 per tenere conto del contributo
               % delle componenti immagine
end 


function X=fft_finestra(x,w)
X=fft(x.*w'); % calcola i coefficienti complessi di Fourier del segnale 
              %finestrato
N=length(x);
S=sum(w);
X=abs(X)/S; % ne calcoliamo il modulo normalizzato per la somma dei pesi 
            % della finestra rettangolare che è quella di default
X=2*X(1:N/2);  % consideriamo per simmetria solo le prime N/2 componenti, 
               % ma moltiplichiamo per 2 per tenere conto del contributo
               % delle componenti immagine
end