From 019dc635948d988f3055304ea9eadbcb6bf7f2d3 Mon Sep 17 00:00:00 2001 From: toni Date: Sun, 15 Oct 2017 21:04:33 +0200 Subject: [PATCH] work on #3 - viele Dinge versucht, leider fehlt irgendwie die richtige Idee... --- matlab/AutoCorrMethodNew_Watch.m | 313 ++++++++++++++++++------------- matlab/findFalseDetectedPeaks.m | 8 +- 2 files changed, 185 insertions(+), 136 deletions(-) diff --git a/matlab/AutoCorrMethodNew_Watch.m b/matlab/AutoCorrMethodNew_Watch.m index aaa4cd6..2ab290d 100644 --- a/matlab/AutoCorrMethodNew_Watch.m +++ b/matlab/AutoCorrMethodNew_Watch.m @@ -6,21 +6,21 @@ %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/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_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_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', ';'); %* ganz schlimm genau die hälfte +%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 @@ -41,163 +41,208 @@ measurements = dlmread('../../measurements/lgWear/recording_48bpm_4-4_176527527. %measurements = dlmread('../measurements/wearR/recording_180bpm_4-4_177064915.csv', ';'); * -%draw the raw acc data -m_idx = []; -m_idx = (measurements(:,2)==2); -m = measurements(m_idx, :); +files = dir(fullfile('../../measurements/lgWear/', '*.csv')); +%files = dir(fullfile('../../measurements/wearR/', '*.csv')); -%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); +for file = files' + + filename = [file.folder '/' file.name]; + measurements = dlmread(filename, ';'); -%put all together again -m = [t_interp', t_interp', m_interp]; + %draw the raw acc data + m_idx = []; + m_idx = (measurements(:,2)==2); + m = measurements(m_idx, :); -figure(1); -plot(m(:,1),m(:,3)) %x -legend("x", "location", "eastoutside"); + %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); -figure(2); -plot(m(:,1),m(:,4)) %y -legend("y", "location", "eastoutside"); + %put all together again + m = [t_interp', t_interp', m_interp]; -figure(3); -plot(m(:,1),m(:,5)) %z -legend("z", "location", "eastoutside"); + figure(1); + plot(m(:,1),m(:,3)) %x + legend("x", "location", "eastoutside"); -%waitforbuttonpress(); + figure(2); + plot(m(:,1),m(:,4)) %y + legend("y", "location", "eastoutside"); -%save timestamps -timestamps = m(:,1); -data = m(:,3); %only z + figure(3); + plot(m(:,1),m(:,5)) %z + legend("z", "location", "eastoutside"); -%TODO: Different window sizes for periods under 16.3 s -window_size = 4096; %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) + %save timestamps + timestamps = m(:,1); + data = m(:,3); %only z - %wait until window is filled with new data - if(mod(i,overlap) == 0) + %TODO: Different window sizes for periods under 16.3 s + window_size = 4096; %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) - %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"); + %wait until window is filled with new data + if(mod(i,overlap) == 0) - corr_x_pos = corr_x; - corr_y_pos = corr_y; - corr_z_pos = corr_z; + %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_pos<0)=0; - corr_y_pos(corr_y_pos<0)=0; - corr_z_pos(corr_z_pos<0)=0; + corr_x_pos = corr_x; + corr_y_pos = corr_y; + corr_z_pos = corr_z; - [peak_x, idx_x] = findpeaks(corr_x_pos, "MinPeakHeight", 0.1,"MinPeakDistance", 50); - [peak_y, idx_y] = findpeaks(corr_y_pos, "MinPeakHeight", 0.1,"MinPeakDistance", 50); - [peak_z, idx_z] = findpeaks(corr_z_pos, "MinPeakHeight", 0.1,"MinPeakDistance", 50); + corr_x_pos(corr_x_pos<0)=0; + corr_y_pos(corr_y_pos<0)=0; + corr_z_pos(corr_z_pos<0)=0; - mean_x = length(peak_x); - mean_y = length(peak_y); - mean_z = length(peak_z); + [peak_x, idx_x] = findpeaks(corr_x_pos, 'MinPeakHeight', 0.1,'MinPeakDistance', 100); + [peak_y, idx_y] = findpeaks(corr_y_pos, 'MinPeakHeight', 0.1,'MinPeakDistance', 100); + [peak_z, idx_z] = findpeaks(corr_z_pos, 'MinPeakHeight', 0.1,'MinPeakDistance', 100); - waitforbuttonpress(); + idx_x = sort(idx_x); + idx_y = sort(idx_y); + idx_z = sort(idx_z); - idx_x = sort(idx_x); - idx_y = sort(idx_y); - idx_z = sort(idx_z); + idx_x = findFalseDetectedPeaks(idx_x, lag_x, corr_x); + idx_y = findFalseDetectedPeaks(idx_y, lag_y, corr_y); + idx_z = findFalseDetectedPeaks(idx_z, 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)); - idx_x = findFalseDetectedPeaks(idx_x, lag_x, corr_x); - idx_y = findFalseDetectedPeaks(idx_y, lag_y, corr_y); - idx_z = findFalseDetectedPeaks(idx_z, lag_z, corr_z); + figure(11); + plot(lag_x, corr_x, lag_x(idx_x), corr_x(idx_x), 'r*') %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"); - 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)); + 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(11); - plot(lag_x, corr_x, lag_x(idx_x), corr_x(idx_x), 'r*') %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"); + figure(12); + plot(lag_y, corr_y, lag_y(idx_y), corr_y(idx_y), 'r*') %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"); - 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)); + 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(12); - plot(lag_y, corr_y, lag_y(idx_y), corr_y(idx_y), 'r*') %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"); + figure(13); + plot(lag_z, corr_z, lag_z(idx_z), corr_z(idx_z), 'r*') %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"); - 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)); + %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); - figure(13); - plot(lag_z, corr_z, lag_z(idx_z), corr_z(idx_z), 'r*') %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"); + 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))); - %window = data(i-window_size:i,:); - %window_timestamps = timestamps(i-window_size:i,:); + num_peaks_x = 1;%length(idx_x); + num_peaks_y = 1;%length(idx_y); + num_peaks_z = 1;%length(idx_z); + + quantity_matrix = [corr_mean_x corr_mean_y corr_mean_z; + corr_rms_x corr_rms_y corr_rms_z; + num_peaks_x num_peaks_y num_peaks_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; + + %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 + + %{ + 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. - %choose axis with most points - if(mean_x > mean_y && mean_x > mean_z) - window_mean_ts_diff = Xwindow_mean_ts_diff; - window_mean_bpm = Xwindow_mean_bpm; - elseif(mean_y > mean_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 - - bpm_per_window_ms = [bpm_per_window_ms, window_mean_ts_diff]; - bpm_per_window = [bpm_per_window, window_mean_bpm]; - - %TODO: choose axis with highest correlation values at peaks - - %TODO: if correlation value is lower then a treshhold, we are not conducting TODO: change to a real classification instead of a treshhold. - - - 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 -%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.4 * 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) - - - - - diff --git a/matlab/findFalseDetectedPeaks.m b/matlab/findFalseDetectedPeaks.m index 792dd4f..c1f633e 100644 --- a/matlab/findFalseDetectedPeaks.m +++ b/matlab/findFalseDetectedPeaks.m @@ -19,8 +19,12 @@ function idx = findFalseDetectedPeaks(idx_orig, lag_orig, corr_orig) %3 %factor_matrix = [1:length(factor_peaks); factor_peaks]; %we need this since octave doesn't have a repelem only repelemS. %diff_peaks_new = repelems(diff_peaks ./ factor_peaks, factor_matrix); - diff_peaks_new = repelem(diff_peaks ./ factor_peaks, factor_peaks); - + if(diff_peaks > 1) + diff_peaks_new = repelem(diff_peaks ./ factor_peaks, factor_peaks); + else + break; + end + %4 idx_new = round([idx_orig(1), cumsum(diff_peaks_new) + idx_orig(1)]'); sum_peaks_new = sum(corr_orig(idx_new));