% Notes and Recommendations:
% - This code calculates the delay between the onset of atrial and ventricular events (AV delay).
% - The program is designed to allow the simultaneous selection of both .mat files, automatically identifying them using the previously saved "ID" variable.
% - The function filterPeaks(Ven, locs, pks, inCua2, FreVi, uhmAm, TiemMin) filters all events and estimates the amplitude of each one.
%   * Requirement: For the function to be detected, it must be located in the same folder as this script.
%   * The variable "uhmAm" is an amplitude threshold that eliminates any event that does not exceed the specified amplitude percentage.
%   * The variable "TiemMin" is a time threshold that eliminates any detected event below this value.
%
% Code Structure:
% The code is divided into the following main stages:
%  1. Simultaneous selection of .mat files (Atria.mat and Ventricles.mat).
%  2. File integrity verification.
%  3. Validation that both .mat files correspond to the same video.
%  4. Inversion of the signals so that the onset of the event represents the highest value.
%  5. Signal normalization.
%  6. Detection of event onsets.
%  7. Filtering of peaks using the "filterPeaks" function.
%  8. Generation of a representative plot of the detection.
%  9. Joint filtering and alignment of the events.
% 10. Calculation of the AV (Atrioventricular) delay.
% 11. Calculation of the mean and standard deviation of the AV delay.


clc
close all
clear

% Prompt user to select atria and ventricles data file and load it
[file,path]=uigetfile('*.mat','Select the Atria.mat and Ventricles.mat','MultiSelect','on');% File generated by the atria and ventricles calculation program
if ~(iscell(file)) % Check if I selected the necessary files.
    error('Select the two files (Atria.mat and ventricles.mat).')
end

file2=char(string(file(1)));
load(fullfile(path, file2)); % Load the first .mat file.
if ID == "Atria" % Check if the first .mat is from the atria or ventricles
    FileA=file3;
elseif ID == "Ventricles"
    FileV=file3;
end

file2=char(string(file(2)));
load(fullfile(path, file2)); % Load the second .mat file.
if ID == "Ventricles" % Check if the second .mat is from the atria or ventricles
    FileV=file3;
elseif ID == "Atria"
    FileA=file3;
end

% Verify that the two .mat files are from the same video and that the files exist.
if ~(string(FileV)==string(FileA)) | ~(exist("FileV")) | ~(exist("FileA"))
    error('The video name is not the same in the two files.')
else
    disp(['Video in... ',file3,' process']);
end

uhmAm=0.1;% Threshold to filter out small peaks
TiemMin=1;% Time threshold to exclude values ​​below this time.
inCua=round(inter*FreVi); % Window adjustable to the stimulation frequency
if rem(inCua,2)==0
    inCua2=(inCua)+1;
else
    inCua2=(inCua);
end

% Invert the signal to detect the start of the event
NorA2=(AFS1+AFS2)/2; % Prom Atria
maxAF=max(NorA2);
NorA2=-(NorA2-maxAF);

maxV=max(VenF(FreVi+1:end));
Ven=-(VenF-maxV);

% Signal normalization
minA=min(NorA2(FreVi+1:end));
maxA=max(NorA2(FreVi+1:end));
Aur=(NorA2-minA)/(maxA-minA);

minV=min(Ven(FreVi+1:end));
maxV=max(Ven(FreVi+1:end));
Ven2=(Ven-minV)/(maxV-minV);

% Signal smoothing
AurF=sgolayfilt(Aur,5,inCua2);
VeF=sgolayfilt(Ven2,5,inCua2);

% Detect peaks for atria and ventricles with minimum peak distance
% Calculation of the peaks of each cardiac event
[Apks,Alocs,Aw,Ap]=findpeaks(AurF,t,'MinPeakDistance',inter*.6);
[weA,~]=filterPeaks(AurF,Alocs,Apks,inCua2,FreVi,uhmAm,TiemMin);
Alocs2=Alocs(weA);
Apks2=Apks(weA);

[Vpks,Vlocs,Vw,Vp]=findpeaks(VeF,t,'MinPeakDistance',inter*.6);
[weV,~]=filterPeaks(VeF,Vlocs,Vpks,inCua2,FreVi,uhmAm,TiemMin);
Vlocs2=Vlocs(weV);
Vpks2=Vpks(weV);

% Plot normalized signals and detected peaks for visual verification
% Graph of the atria, ventricles and detected peaks
figure
plot(t,VeF)
hold on
plot(t,AurF)
hold on
plot(Vlocs2,Vpks2,'o')
hold on
plot(Alocs2,Apks2,'o')

% AV delay calculation
TamV=length(Vpks2); % Number of events
a1=0;
VecAV=[];
Vlocs3=Vlocs2;
Vpks3=Vpks2;
Alocs3=Alocs2;
Apks3=Apks2;
ELV=[];
ELA=[];
% Filter the events to reduce the error
for i2=1:TamV
    if i2==1 % Verify that the first event appearing is an atrial event; if not, delete that ventricular event.
        we5=find(Vlocs2<=Alocs2(i2));
        Lwe5=length(we5);
        if Lwe5>1
            ELV=[ELV,we5(1:end-1)];
        elseif Lwe5==1
            ELV=[ELV,we5];
        else
            we1=find(Alocs2<=Vlocs2(i2));
            Lwe1=length(we1);
            if Lwe1>1
                ELA=[ELA,we1(1:end-1)];
            end
        end
    else
        we3=find(Alocs2>=Vlocs2(i2-1) & Alocs2<=Vlocs2(i2)); % Verify that each ventricular event is preceded by an atrial event; 
        Lwe3=length(we3);
        if Lwe3>1 % if there is more than one ventricular event, only the first is considered and the others are discarded.
            ELA=[ELA,we3(1:end-1)];
        elseif Lwe3==0 % if there is no atrial event, that ventricular event is discarded.
            ELV=[ELV,i2];
        end
    end
end
Vlocs3(ELV)=[];
Vpks3(ELV)=[];
Alocs3(ELA)=[];
Apks3(ELA)=[];

TamV2=length(Vpks3); % Number of events after filtering
for i=1:TamV2 % Calculation of the AV delay for each event
    we3=find(Alocs3<Vlocs3(i));
    AVref=Vlocs3(i)-Alocs3(we3(end));
    VecAV=[VecAV,AVref]; % Vector of AV delay values
end
we4=find(VecAV<=0 | VecAV>=.1); % Only consider delays between 0 - 100 ms
VecAV(we4)=[];
PAV=mean(VecAV) % Average AV Delay
DAV=std(VecAV)  % Standard deviation of AV delay