%We are using a threshold-based version for bpm estimation %Only the z axis of the acc is used %NOTE: depending on the measurement device we have a highly different sample rate. the smartwatches are not capable of providing a constant sample rate. the xsens on the other hand is able to do this. %load file provided by the sensor readout app % SMARTWATCH LG WEAR ------> 100 hz - 1000hz %measurements = dlmread('../../measurements/lgWear/PR_recording_80bpm_4-4_177596720.csv', ';'); % %measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176527527.csv', ';'); %measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176606785.csv', ';'); %measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176696356.csv', ';'); %measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176820066.csv', ';'); % %measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176931941.csv', ';'); %double %measurements = dlmread('../../measurements/lgWear/recording_72bpm_4-4_176381633.csv', ';'); %measurements = dlmread('../../measurements/lgWear/recording_72bpm_4-4_176453327.csv', ';'); % %measurements = dlmread('../../measurements/lgWear/recording_100bpm_4-4_176073767.csv', ';'); %* %measurements = dlmread('../../measurements/lgWear/recording_100bpm_4-4_176165357.csv', ';'); %measurements = dlmread('../../measurements/lgWear/recording_100bpm_4-4_176230146.csv', ';'); %measurements = dlmread('../../measurements/lgWear/recording_100bpm_4-4_176284687.csv', ';'); % %measurements = dlmread('../../measurements/lgWear/recording_100bpm_4-4_177368860.csv', ';'); %(besonders) %measurements = dlmread('../../measurements/lgWear/recording_180bpm_4-4_177011641.csv', ';'); % %measurements = dlmread('../../measurements/lgWear/recording_180bpm_4-4_177064915.csv', ';'); % % SMARTWATCH G WATCH WEAR R ----> 100hz - 250hz %measurements = dlmread('../measurements/wearR/PR_recording_80bpm_4-4_177596720.csv', ';'); %* %measurements = dlmread('../measurements/wearR/recording_48bpm_4-4_176527527.csv', ';'); %measurements = dlmread('../measurements/wearR/recording_48bpm_4-4_176606785.csv', ';'); %measurements = dlmread('../../measurements/wearR/recording_48bpm_4-4_176696356.csv', ';'); %* %measurements = dlmread('../measurements/wearR/recording_48bpm_4-4_176820066.csv', ';'); %measurements = dlmread('../measurements/wearR/recording_48bpm_4-4_176931941.csv', ';'); %double %measurements = dlmread('../measurements/wearR/recording_72bpm_4-4_176381633.csv', ';'); %measurements = dlmread('../measurements/wearR/recording_72bpm_4-4_176453327.csv', ';'); %measurements = dlmread('../measurements/wearR/recording_100bpm_4-4_176073767.csv', ';'); * %measurements = dlmread('../measurements/wearR/recording_100bpm_4-4_176165357.csv', ';'); %measurements = dlmread('../measurements/wearR/recording_100bpm_4-4_176230146.csv', ';'); %* 72? %measurements = dlmread('../measurements/wearR/recording_100bpm_4-4_176284687.csv', ';'); %*48? %measurements = dlmread('../measurements/wearR/recording_100bpm_4-4_177368860.csv', ';'); %(besonders) %measurements = dlmread('../measurements/wearR/recording_180bpm_4-4_177011641.csv', ';'); * %measurements = dlmread('../measurements/wearR/recording_180bpm_4-4_177064915.csv', ';'); * %files = dir(fullfile('../../measurements/lgWear/', '*.csv')); %files = dir(fullfile('../../measurements/wearR/', '*.csv')); files = dir(fullfile('../../measurements/peter_failed/', '*.csv')); for file = files' filename = [file.folder '/' file.name]; measurements = dlmread(filename, ';'); %draw the raw acc data m_idx = []; m_idx = (measurements(:,2)==10); m = measurements(m_idx, :); %Interpolate to generate a constant sample rate to 250hz (4ms per sample) sample_rate_ms = 4;%ms [~, m_unique_idx] = unique(m(:,1)); %matlab requirs unique timestamps for interp m = m(m_unique_idx, :); t = m(:,1); %timestamps t_interp = t(1):sample_rate_ms:t(length(t)); m_interp = interp1(t,m(:,3:5),t_interp); %put all together again m = [t_interp', t_interp', m_interp]; figure(1); plot(m(:,1),m(:,3)) %x legend("x", "location", "eastoutside"); figure(2); plot(m(:,1),m(:,4)) %y legend("y", "location", "eastoutside"); figure(3); plot(m(:,1),m(:,5)) %z legend("z", "location", "eastoutside"); waitforbuttonpress(); %save timestamps timestamps = m(:,1); data = m(:,3); %only z %TODO: Different window sizes for periods under 16.3 s window_size = 2048; %about 2 seconds using 2000hz, 16.3 s using 250hz overlap = 256; bpm_per_window_ms = []; bpm_per_window = []; for i = window_size+1:length(data) %wait until window is filled with new data if(mod(i,overlap) == 0) %measure periodicity of window and use axis with best periodicity [corr_x, lag_x] = xcov(m(i-window_size:i,3), (window_size/4), "coeff"); [corr_y, lag_y] = xcov(m(i-window_size:i,4), (window_size/4), "coeff"); [corr_z, lag_z] = xcov(m(i-window_size:i,5), (window_size/4), "coeff"); corr_x_pos = corr_x; corr_y_pos = corr_y; corr_z_pos = corr_z; corr_x_pos(corr_x_pos<0)=0; corr_y_pos(corr_y_pos<0)=0; corr_z_pos(corr_z_pos<0)=0; [peak_x, idx_x_raw] = findpeaks(corr_x_pos, 'MinPeakHeight', 0.1,'MinPeakDistance', 50, 'MinPeakProminence', 0.1); [peak_y, idx_y_raw] = findpeaks(corr_y_pos, 'MinPeakHeight', 0.1,'MinPeakDistance', 50, 'MinPeakProminence', 0.1); [peak_z, idx_z_raw] = findpeaks(corr_z_pos, 'MinPeakHeight', 0.1,'MinPeakDistance', 50, 'MinPeakProminence', 0.1); idx_x_raw = sort(idx_x_raw); idx_y_raw = sort(idx_y_raw); idx_z_raw = sort(idx_z_raw); idx_x = findFalseDetectedPeaks(idx_x_raw, lag_x, corr_x); idx_y = findFalseDetectedPeaks(idx_y_raw, lag_y, corr_y); idx_z = findFalseDetectedPeaks(idx_z_raw, lag_z, corr_z); Xwindow = m(i-window_size:i,3); Xwindow_mean_ts_diff = mean(diff(lag_x(idx_x) * sample_rate_ms)); %2.5 ms is the time between two samples at 400hz Xwindow_mean_bpm = (60000 / (Xwindow_mean_ts_diff)); figure(11); plot(lag_x, corr_x, lag_x(idx_x), corr_x(idx_x), 'r*', lag_x(idx_x_raw), corr_x(idx_x_raw), 'g*') %z hold ("on") m_label_ms = strcat(" mean ms: ", num2str(Xwindow_mean_ts_diff)); m_label_bpm = strcat(" mean bpm: ", num2str(Xwindow_mean_bpm)); title(strcat(" ", m_label_ms, " ", m_label_bpm)); hold ("off"); Ywindow = m(i-window_size:i,4); Ywindow_mean_ts_diff = mean(diff(lag_y(idx_y) * sample_rate_ms)); Ywindow_mean_bpm = (60000 / (Ywindow_mean_ts_diff)); figure(12); plot(lag_y, corr_y, lag_y(idx_y), corr_y(idx_y), 'r*', lag_y(idx_y_raw), corr_y(idx_y_raw), 'g*') %z hold ("on") m_label_ms = strcat(" mean ms: ", num2str(Ywindow_mean_ts_diff)); m_label_bpm = strcat(" mean bpm: ", num2str(Ywindow_mean_bpm)); title(strcat(" ", m_label_ms, " ", m_label_bpm)); hold ("off"); Zwindow = m(i-window_size:i,5); Zwindow_mean_ts_diff = mean(diff(lag_z(idx_z)* sample_rate_ms)); Zwindow_mean_bpm = (60000 / (Zwindow_mean_ts_diff)); figure(13); plot(lag_z, corr_z, lag_z(idx_z), corr_z(idx_z), 'r*', lag_z(idx_z_raw), corr_z(idx_z_raw), 'g*') %z hold ("on") m_label_ms = strcat(" mean ms: ", num2str(Zwindow_mean_ts_diff)); m_label_bpm = strcat(" mean bpm: ", num2str(Zwindow_mean_bpm)); title(strcat(" ", m_label_ms, " ", m_label_bpm)); hold ("off"); %breakpoints dummy for testing if(length(idx_x) > length(idx_x_raw)) a = 0; %breakpointdummy end if(length(idx_y) > length(idx_y_raw)) a = 0; %breakpointdummy end if(length(idx_z) > length(idx_z_raw)) a = 0; %breakpointdummy end %Find the most proper axis. We use 3 quantities: mean of corr. %value, sum of corr val. and number of peaks. Simple normalization %to get the axis that fullfills the quantities the most. idx_noZero_x = idx_x(lag_x(idx_x) ~= 0); idx_noZero_y = idx_y(lag_x(idx_y) ~= 0); idx_noZero_z = idx_z(lag_x(idx_z) ~= 0); corr_mean_x = geomean(corr_x(idx_noZero_x(corr_x(idx_noZero_x)>0))); corr_mean_y = geomean(corr_y(idx_noZero_y(corr_y(idx_noZero_y)>0))); corr_mean_z = geomean(corr_z(idx_noZero_z(corr_z(idx_noZero_z)>0))); corr_rms_x = rms(corr_x(idx_x(lag_x(idx_x) ~= 0))); corr_rms_y = rms(corr_y(idx_y(lag_y(idx_y) ~= 0))); corr_rms_z = rms(corr_z(idx_z(lag_z(idx_z) ~= 0))); num_peaks_x = 1;%length(idx_x); num_peaks_y = 1;%length(idx_y); num_peaks_z = 1;%length(idx_z); num_intersection_x = getNumberOfIntersections(corr_x, lag_x, 0.2); num_intersection_y = getNumberOfIntersections(corr_y, lag_y, 0.2); num_intersection_z = getNumberOfIntersections(corr_z, lag_z, 0.2); quantity_matrix = [corr_mean_x corr_mean_y corr_mean_z; corr_rms_x corr_rms_y corr_rms_z; num_intersection_x num_intersection_y num_intersection_z]; quantity_matrix_percent(1,:) = quantity_matrix(1,:) ./ sum(quantity_matrix(1,:)); quantity_matrix_percent(2,:) = quantity_matrix(2,:) ./ sum(quantity_matrix(2,:)); quantity_matrix_percent(3,:) = quantity_matrix(3,:) ./ sum(quantity_matrix(3,:)); quantity_factors = sum(quantity_matrix_percent) / 3; %TODO: Wenn ein quantity wert NaN ist, sind alle NaN... quantity_x = quantity_factors(1); quantity_y = quantity_factors(2); quantity_z = quantity_factors(3); %choose axis with sum(corr) nearest to 0 %{ corr_sum_xyz = [sum(corr_x) sum(corr_y) sum(corr_z)]; [~,idx_nearest_zero] = min(abs(corr_sum_xyz)); if(idx_nearest_zero == 1) window_mean_ts_diff = Xwindow_mean_ts_diff; window_mean_bpm = Xwindow_mean_bpm; elseif(idx_nearest_zero == 2) window_mean_ts_diff = Ywindow_mean_ts_diff; window_mean_bpm = Ywindow_mean_bpm; else window_mean_ts_diff = Zwindow_mean_ts_diff; window_mean_bpm = Zwindow_mean_bpm; end %} %quantity_x = num_intersection_x; %quantity_y = num_intersection_y; %quantity_z = num_intersection_z; if(quantity_x > quantity_y && quantity_x > quantity_z) window_mean_ts_diff = Xwindow_mean_ts_diff; window_mean_bpm = Xwindow_mean_bpm; elseif(quantity_y > quantity_z) window_mean_ts_diff = Ywindow_mean_ts_diff; window_mean_bpm = Ywindow_mean_bpm; else window_mean_ts_diff = Zwindow_mean_ts_diff; window_mean_bpm = Zwindow_mean_bpm; end if(isnan(window_mean_ts_diff) || isnan(window_mean_bpm)) %do nothing else bpm_per_window_ms = [bpm_per_window_ms, window_mean_ts_diff]; bpm_per_window = [bpm_per_window, window_mean_bpm]; end %TODO: if correlation value is lower then a treshhold, we are not conducting TODO: change to a real classification instead of a treshhold. end end %TODO: smooth the results using a moving avg or 1d kalman filter.(transition for kalman could be adding the last measured value) %remove the first 40% of the results, due to starting delays while recording. number_to_remove = round(abs(0.1 * length(bpm_per_window_ms))); num_all = length(bpm_per_window_ms); bpm_per_window_ms = bpm_per_window_ms(number_to_remove:num_all); bpm_per_window = bpm_per_window(number_to_remove:num_all); mean_final_ms = mean(bpm_per_window_ms); std_final_ms = std(bpm_per_window_ms); mean_final_bpm = mean(bpm_per_window); std_final_bpm = std(bpm_per_window); fprintf('%s: mean = %f bpm (%f ms) stddev = %f bpm (%f ms)\n', strrep(regexprep(filename,'^.*recording_',''),'.txt',''), mean_final_bpm, mean_final_ms, std_final_bpm, std_final_ms); end