function [] = ekganalysis(filename, f, year, month, day, hour, minute, second)
% A function to open and interpret the binary data file generated from the
% embedded ECG logger. The date funtionality makes Matlab lag like there's
% no tomorrow, so that's currently disabled.
%
% filename = binary file created by the logger located in
% Matlab's seach path.
%
% f = sampling frequency. 256Hz or 512Hz.
%
% The rest is self explainatory.
%
% example:
%
% >> ekganalysis('ekgdata8.bin', 512, 2014, 07, 27, 22, 61, 0);
%
% Where 2014-07-22 22:61:00 is the date the user started his logging session.
buffer = 0;
timestamps = 0;
datevector = [0 0 0 0 0 0];
fid = fopen(filename,'rb'); % Open the binary file for reading with file pointer fid
[data,count] = fread(fid, 'uint16'); % Scan the data into a vector, in this case called data
%Start the seach for saved timestamps and store them in buffer.
%Binary files with timestamps will contain 0xFFFF near the end followed by
%32 bits timestamps in seconds since the logging session started.
%The timestamps are for convenience stored as two 16-bits values.
for i = length(data) : -1: 1
if(data(i) == 65535)
for j = i : length(data)
buffer(end+1) = data(i);
data(i) = [];
end
break;
end
end
count = length(data);
%Check if file contains timestamps
if( length(buffer) > 1)
buffer(1:2) = [];
%Merge the 16-bits values into 32-bits values.
for i = 1:2: length(buffer)
x = uint16(buffer(i));
y = uint16(buffer(i+1));
bytepack = uint32(x);
bytepack = bitshift(bytepack, 16);
timestamps(end+1)= bitor(bytepack, uint32(y));
end
end
timestamps(1) = [];
timestamps = fliplr(timestamps);
timestamps = timestamps * 512;
%Produce a vector of seconds instead of samples
n = linspace(0, count/f, count);
%seconds = floor(count/f);
%minutes = floor(seconds / 60);
%hours = floor(minutes / 60);
%for i = 0: seconds
% datevector(end+1,:) = [year, month, day, hour, minute, second+i]
%end
%datevector(1,:) = [];
%datenumber = datenum(datevector)
%nnn = linspace(datenumber(1), datenumber(seconds), count);
%Plot the data
subplot(2,1,1);
plot(n, data);
xlabel('seconds');
%datetick('x', 0);
%Plot the timestamps
hold on;
for i = 1 : length(timestamps)
x=[timestamps(i)/f,timestamps(i)/f];
y=[0,2^12*1.1];
plot(x,y, 'red');
end
hold off;
%Analyse the frequency spectrum
XS = fft(data,f);
%Plotting the FFT. The 0-component of the FFT is the DC component which
%will be massive. With Matlab's 1 indexing (the index of arrays start at 1
%instead of 0), we can remove this component from the plot by starting the
%plot from position 2. Since the FFT is repeating and symmetrical, we also
%just need the first half of the spectrum, 1-f/2 Hz.
subplot(2,1,2);
stem(abs(XS(2:f/2)), 'filled');
xlabel('Hz');