%hr = heartRateEstimation(imgDirectory, fps, windowDuration, windowShift)
%%%%%%%%%%%%%%%%%
% hr = heartRateEstimation(imgDirectory, fps, windowDuration, windowShift)
%
% Task: Finding the heartrate of person from selcted images
%
% Inputs:
%	-imgDirectory
% -fps
%
% Output:
%	-heartrate:
%
% Author: Thomas Périn
% Date: 17/02/2023
%
% Note: images were extracted with ffmpeg -i "video.mkv" ../img/gg_%04d.png
%%%%%%%%%%%%%%%%%

clc
clear all

% temp var
fps = 15;
windowDuration = 30;
roi = [280 580 600 850];

% global
redChannel = [];
greenChannel = [];
blueChannel = [];

%list the images with i the same directory
listImg = dir(['*.png']);

% iteration for every images available
for l_img=1:windowDuration*fps
	% load the current image
  image_original = imread([listImg(l_img).name]);

	% crop the current image around the face
	image_face = image_original(roi(1):roi(2), roi(3):roi(4), :);

	% spatial average for r, g, b channels
	r_mean = mean(mean(image_face(:,:,1)));
	g_mean = mean(mean(image_face(:,:,2)));
	b_mean = mean(mean(image_face(:,:,3)));

	% store the current "average" pixel in global arrays
	redChannel = [redChannel r_mean];
	greenChannel = [greenChannel g_mean];
	blueChannel = [blueChannel b_mean];
end

% estimate temporal average and standard deviation
redChannel_avg = mean(redChannel);
redChannel_std = std(redChannel);
greenChannel_avg = mean(greenChannel);
greenChannel_std = std(greenChannel);
blueChannel_avg = mean(blueChannel);
blueChannel_std = std(blueChannel);

% normalize your data
% first define length of our data
L = columns(greenChannel);
for i=1:L
  greenChannel_normalized(i) = (greenChannel(i) - greenChannel_avg)/greenChannel_std;
end

% select only greenChannel for simplicity (no ICA)
greenChannel_fft = fft(greenChannel_normalized);

% power spectrum using link(https://www.mathworks.com/help/matlab/ref/fft.html)
P2 = abs(greenChannel_fft/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = 15*(0:(L/2))/L;
plot(f,P1)

% find the peak in the range ([0.75 4] Hz)
% define range and round to get intergers
low_index = round(0.75*L/15);
high_index = round(2*L/15);
% find max value and index inn new range
[a, ia] = max(P1(low_index:high_index));
% adapt the index for full range
heart_frequency = f(ia+low_index)

% determine heart rate
% convert the index from Hz to bpm
heart_rate = heart_frequency*60