HR estimation using camera-based technique and compare it with ground truth ECG

script_infrared_driving.m

-%Infrared_driving_face
+%Infrared_driving_face_ROI
 clc; clear all; close all;
 addpath(genpath('E:\TU_braunschweig\work2\examp\infrared_driving\signal'));
 cd 'E:\TU_braunschweig\work2\examp\infrared_driving\mugam'
@@ -16,11 +16,9 @@ overlap_length = 10; % overlapping 10 sec signal
 
 [b,a] = butter(4,[0.75 2]/(fs/2)); % Heart rate % bidirectional filter bandpass butterworth filter
  
-for j = 4
-    %length(files_cam)
+for j = 1:length(files_cam)
     out(j) = load(files_cam(j).name);
-    x = mean(out(j).signal,2); % averaging of the signal- 3 channels
-%     x = out(j).signal(:,3);
+    x = mean(out(j).signal,2); % averaging of all channels
     t = 1:length(x);
     time = (t-1)/(fs*60); % Time in minute
 %     figure,
@@ -31,6 +29,7 @@ for j = 4
     % Preprocessing
     face_sig1 = detrend(x); % detrend the signal
     face_sig = filtfilt(b,a,face_sig1); % Filtered signal in the range of [0.75 2 Hz]
+    
     % 60 sec signal with overlapping window of 10 sec signal
     face = overlap_window_signal(x,fs,window_length,overlap_length);
     face_sig_overlap = overlap_window_signal(face_sig,fs,window_length,overlap_length);
@@ -38,9 +37,10 @@ for j = 4
     
     filename_gt = remove_extension(files_gt(j).name);
     filename_camera = remove_extension(files_cam(j).name);
-    files_gt(j).name
-    files_cam(j).name
+
     if filename_gt == filename_camera
+        
+        % Ground truth
         out_gt{j} = dlmread(files_gt(j).name);
         Mx = out_gt{1, j}(:,2); % ECG signal - Eindhoven II lead
         x1 = out_gt{1, j}(:,1); 
@@ -52,9 +52,8 @@ for j = 4
         time1_overlap = overlap_window_signal(time1,sf,window_length,overlap_length);
     
         % power spectrum
-        for k = 34
-            %size(face_sig_overlap,2) 
-            %34
+        for k = 1:size(face_sig_overlap,2) 
+      
             face_overlap = face_sig_overlap(:,k);
             [pxx,f,peak_f,max_f] = freq_spectrum(face_overlap,fs); % periodogram
             HR_infrared_driving(j,k) = round((peak_f*60),1); % HR in face signal
@@ -64,39 +63,39 @@ for j = 4
             
             
             % Visualization
-            figure,plot(time_overlap(:,k),face(:,k),'k','LineWidth',2);
-            xlabel('Time (min)','FontWeight','bold','FontSize',20)
-            ylabel('Amplitude(a.u.)','FontWeight','bold','FontSize',20)
+            %figure,plot(time_overlap(:,k),face(:,k),'k','LineWidth',2);
+            %xlabel('Time (min)','FontWeight','bold','FontSize',20)
+            %ylabel('Amplitude(a.u.)','FontWeight','bold','FontSize',20)
 %             ylim([110 140])
-            set(gca,'LooseInset',get(gca,'TightInset'));
-            title('Raw signal','FontWeight','bold','FontSize',15)
+            %set(gca,'LooseInset',get(gca,'TightInset'));
+            %title('Raw signal','FontWeight','bold','FontSize',15)
             
             
-            figure,plot(time_overlap(:,k), face_overlap,'k','LineWidth',2)
-            xlabel('Time (min)','FontWeight','bold','FontSize',20)
-            ylabel('Amplitude(a.u.)','FontWeight','bold','FontSize',20)
-            set(gca,'LooseInset',get(gca,'TightInset'));
-            title('Preprocessed signal','FontWeight','bold','FontSize',15)
+            %figure,plot(time_overlap(:,k), face_overlap,'k','LineWidth',2)
+            %xlabel('Time (min)','FontWeight','bold','FontSize',20)
+            %ylabel('Amplitude(a.u.)','FontWeight','bold','FontSize',20)
+            %set(gca,'LooseInset',get(gca,'TightInset'));
+            %title('Preprocessed signal','FontWeight','bold','FontSize',15)
             
             
-            figure,plot(f,pxx,'k','LineWidth',2)
+            %figure,plot(f,pxx,'k','LineWidth',2)
 %             legend
-            hold on
-            plot(peak_f,max_f,'rp','MarkerSize', 10);
-            xlim([0.4 2.1])
-            xlabel('Frequency (Hz)','FontWeight','bold','FontSize',20)
-            ylabel('Amplitude(a.u.)','FontWeight','bold','FontSize',20)
-            set(gca,'LooseInset',get(gca,'TightInset'));
-            title(['Camera: HR is ',num2str(HR_infrared_driving(j,k)), 'bpm'],'FontWeight','bold','FontSize',15)
-            hold off
+            %hold on
+            %plot(peak_f,max_f,'rp','MarkerSize', 10);
+            %xlim([0.4 2.1])
+            %xlabel('Frequency (Hz)','FontWeight','bold','FontSize',20)
+            %ylabel('Amplitude(a.u.)','FontWeight','bold','FontSize',20)
+            %set(gca,'LooseInset',get(gca,'TightInset'));
+            %title(['Camera: HR is ',num2str(HR_infrared_driving(j,k)), 'bpm'],'FontWeight','bold','FontSize',15)
+            %hold off
 
-            figure,plot(time1_overlap(:,k),groundtruth(:,k),'k','LineWidth',2); hold on;
-            offset = 0.01;
-            plot((R_peaks_final+(k*fs*sf)-(sf*fs))/(sf*60), max(groundtruth(:,k))+offset*2, 'rp', 'MarkerSize', 10)
-            xlabel('Time (min)','FontWeight','bold','FontSize',20)
-            ylabel('Amplitude(a.u.)','FontWeight','bold','FontSize',20)
-            set(gca,'LooseInset',get(gca,'TightInset'));
-            title(['Ground truth: HR is ',num2str(Heart_rate), 'bpm'],'FontWeight','bold','FontSize',15)
+            %figure,plot(time1_overlap(:,k),groundtruth(:,k),'k','LineWidth',2); hold on;
+            %offset = 0.01;
+            %plot((R_peaks_final+(k*fs*sf)-(sf*fs))/(sf*60), max(groundtruth(:,k))+offset*2, 'rp', 'MarkerSize', 10)
+            %xlabel('Time (min)','FontWeight','bold','FontSize',20)
+            %ylabel('Amplitude(a.u.)','FontWeight','bold','FontSize',20)
+            %set(gca,'LooseInset',get(gca,'TightInset'));
+            %title(['Ground truth: HR is ',num2str(Heart_rate), 'bpm'],'FontWeight','bold','FontSize',15)
 
         end 
     end
@@ -113,17 +112,10 @@ end
 % figure, plot(Mean_camera,'o','DisplayName','camera')
 % hold on, plot(Mean_gt,'o','DisplayName','GT')
 
-% %% Bland altman plot
+% %% Bland altman plot & correlation plot
 % data1 = reshape(HR_infrared,[],1);
 % data2 = reshape(HR_gt,[],1);
 % MakeBlandAltmanPlots(data1,data2,'HR','HR',TitleForPlots)
-% bland_altman( data1,data2)
 
 
-% for j = 1:length(files_gt) 
-%     out_gt{j} = dlmread(files_gt(j).name);
-%         Mx = out_gt{1, j}(:,2); % ECG signal - Eindhoven II lead
-%         x1 = out_gt{1, j}(:,1); 
-%         time1 = (x1-1)/(sf*60); % time in minutes
-%         figure,plot(time1,Mx)
-% end
\ No newline at end of file
+% end