Upload New File

script_infrared_driving.m

+%Infrared_driving_face
+clc; clear all; close all;
+addpath(genpath('E:\TU_braunschweig\work2\examp\infrared_driving\signal'));
+cd 'E:\TU_braunschweig\work2\examp\infrared_driving\mugam'
+addpath('E:\TU_braunschweig\work2\examp\infrared_driving\mugam')
+files_cam = dir('*.mat');
+fs = 10; %frame rate 
+
+% Ground truth
+cd 'E:\TU_braunschweig\work2\examp\infrared_driving\signal'
+files_gt = dir('*.txt');
+sf = 500; % sampling frequency of ECG
+
+window_length = 60; % window length 60 sec of the signal
+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)
+    out(j) = load(files_cam(j).name);
+    x = mean(out(j).signal,2); % averaging of the signal- 3 channels
+%     x = out(j).signal(:,3);
+    t = 1:length(x);
+    time = (t-1)/(fs*60); % Time in minute
+%     figure,
+%     plot(time, x,'k','LineWidth',2)
+%     xlabel('Time (min)')
+%     ylabel ('Amplitude')
+
+    % 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);
+    time_overlap = overlap_window_signal(time,fs,window_length,overlap_length);
+    
+    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
+        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
+        % Conversion from raw data to millivolts signal
+        M1 = convert_mv(Mx);
+        
+        groundtruth = overlap_window_signal(M1,sf,window_length,overlap_length);
+        time1_overlap = overlap_window_signal(time1,sf,window_length,overlap_length);
+    
+        % power spectrum
+        for k = 34
+            %size(face_sig_overlap,2) 
+            %34
+            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
+            
+            [R_peaks_final, Heart_rate] = script_infrared_driving_ECG(groundtruth(:,k),sf);
+            HR_gt(j,k) = Heart_rate;
+            
+            
+            % 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)
+%             ylim([110 140])
+            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(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
+
+            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
+    Mean_camera(j,:) = mean(HR_infrared_driving(j,:));
+    Mean_gt(j,:) = mean(HR_gt(j,:));
+    Std_camera(j,:) = std(HR_infrared_driving(j,:));
+    Std_gt(j,:) = std(HR_gt(j,:));
+    RMSE(j,:) = sqrt(mean((HR_gt(j,:) - HR_infrared_driving(j,:)).^2));
+    mean_abs_error(j,:) = mean(abs((HR_gt(j,:) - HR_infrared_driving(j,:))));
+end
+
+% figure, plot(RMSE,'o','DisplayName','RMSE')
+% figure, plot(mean_abs_error,'o','DisplayName','MAE')
+% figure, plot(Mean_camera,'o','DisplayName','camera')
+% hold on, plot(Mean_gt,'o','DisplayName','GT')
+
+% %% Bland altman 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