%%  Comparison of my peak detect alfgorithm vs. the one proposed by Physionet
% The comparison is on physionet database - AFIB patients

%%  Cleaning
clc
clear all
close all

%% 
Path='Physionet_AFIB\';
Record_name='p44c';
%% Load annotations
 fileID = fopen([Path 'an' Record_name '.txt'],'r');
 formatSpec = '%*s%d%*s%*d%*d%*d\n';
 Annotations = fscanf(fileID,formatSpec)
 fclose(fileID);
 
 %%  Load signals
   fileID = fopen([Path Record_name '.txt'],'r');
 formatSpec = '%d:%*f%*f%*f\n';
 [Ts_min] = fscanf(fileID,formatSpec);
  fclose(fileID);
  fileID = fopen([Path Record_name '.txt'],'r');
 formatSpec = '%*d:%f%*f%*f\n';
 [Ts] = fscanf(fileID,formatSpec);
  fclose(fileID);
  fileID = fopen([Path  Record_name '.txt'],'r'); 
 formatSpec = '%*d:%*f\t%f%*f\n';
 [Signal] = fscanf(fileID,formatSpec);
 fclose(fileID);
 Ts=Ts+Ts_min*60;
 Annotations(Annotations>length(Signal))=[];
 %%  Place a Physionet annotation
 plot(Ts,Signal)
 hold all
 plot(Ts(Annotations),Signal(Annotations),'ro');
 
 %%  Apply my algorithm on the signal
         %% Method suggested by me
        handles=[];
        sample_rate=1/(Ts(2)-Ts(1)); % Hz
        peakInds=peakfinder(Signal(:,1)',sample_rate,handles);    
        %figure(2)
        %plot(tm,sig(:,1));hold on;grid on
        plot(Ts(peakInds),Signal(peakInds,1),'rx')
        xlabel('[Sec]')
        ylabel('Voltage [NU]');
        %title('Method suggested by me')
        legend('ECG','Method suggested by Physionet','Method suggested by me')
        drawnow
        
%%  Statistics of compatibility with Physionet
%%  Variables
        false_positive=0;
        true_positive=0;
        false_negative=0;
        true_negative=0;
        ann=Annotations;
%%  Definitions
allowed_delta=0.02; %[sec]
allowed_samples_ann_diff=round(sample_rate*allowed_delta);

        for i=1:length(ann) % Taking Physweb as gold standard
           I = find((peakInds==ann(i))|(peakInds==ann(i)-1)|(peakInds==ann(i)+1)|(peakInds-1==ann(i))|(peakInds+1==ann(i)));
           if length(I)>0
               true_positive=true_positive+1;
           else
               false_negative=false_negative+1;
               ann_i=ann(i)
           end
        end
      for i=1:length(peakInds) % Taking Physweb as gold standard
           I = find((ann==peakInds(i))|(ann==peakInds(i)-1)|(ann==peakInds(i)+1)|(ann-1==peakInds(i))|(ann+1==peakInds(i))|(ann==peakInds(i)-2)|(ann==peakInds(i)+2));
           if length(I)==0
               false_positive=false_positive+1;
              peakInds_i= peakInds(i);
           end
     end

        true_negative=length(Annotations)-false_negative-false_positive-true_positive;
 
        R_peaks_number=length(Annotations);
        
      
 