diff --git a/AnnoResult.m b/AnnoResult.m index 9224c40..59fb1de 100644 --- a/AnnoResult.m +++ b/AnnoResult.m @@ -106,7 +106,7 @@ obj.missing_lead = num2cell(missing_lead); % VERSION MANUALLY UPDATED HERE - obj.version = {'1.2.2'}; + obj.version = {'1.3.0'}; end end diff --git a/Annoparams.m b/Annoparams.m index 3806371..d0a8353 100644 --- a/Annoparams.m +++ b/Annoparams.m @@ -40,6 +40,7 @@ % Load in ECG/peak detection maxBPM = 150; % Sets window for detecting R peaks pkthresh = 95; % Percentile of ECG signal above which a peak can be found + pkfilter = 0; % Additional filtering prior to peak detection on/off lowpass = 1; % Low pass wavelet filter on/off highpass = 1; % High pass wavelet filter on/off wavelet_level_lowpass = 1; % LPF freq is > samp freq/2^(wavelet_level_lowpass + 1) diff --git a/VCG.m b/VCG.m index 2a9f1f4..d059f1b 100644 --- a/VCG.m +++ b/VCG.m @@ -110,7 +110,7 @@ function l = length(obj); l = length(obj.VM); end - function QRS = peaks(obj,a); QRS = findpeaksecg(obj.VM, a.maxBPM, obj.hz, a.pkthresh); end + function QRS = peaks(obj,a); QRS = findpeaksecg(obj.VM, a.maxBPM, obj.hz, a.pkthresh, a.pkfilter); end function e = ecg(obj, transform_matrix) switch transform_matrix diff --git a/braveheart_batch_windows.prj b/braveheart_batch_windows.prj index 71e026a..a7dff0c 100644 --- a/braveheart_batch_windows.prj +++ b/braveheart_batch_windows.prj @@ -8,13 +8,13 @@ ${PROJECT_ROOT}\braveheart_batch_windows_resources\icon_24.png ${PROJECT_ROOT}\braveheart_batch_windows_resources\icon_16.png - 1.2.2 + 1.3.0 Hans F. Stabenau and Jonathan W. Waks braveheart.ecg@gmail.com BRAVHEART ECG/VCG Analysis Software -- Command Line Version -Version 1.2.2 +Version 1.3.0 ${PROJECT_ROOT}\braveheart_splash.png \braveheart_batch\ @@ -50,7 +50,7 @@ Version 1.2.2 varargin input arguments Description BRAVHEART ECG/VCG Analysis Software -- Command Line Version -Version 1.2.2 +Version 1.3.0 @@ -86,6 +86,7 @@ Version 1.2.2 ${PROJECT_ROOT}\ecg_axis.m ${PROJECT_ROOT}\elecpos.mat ${PROJECT_ROOT}\findHighpassLvl.m + ${PROJECT_ROOT}\findpeakswavelet.m ${PROJECT_ROOT}\getfiletype.m ${PROJECT_ROOT}\HRV_Calc.m ${PROJECT_ROOT}\leadidcodexyz.m diff --git a/braveheart_gui.fig b/braveheart_gui.fig index 47df29c..d365c77 100644 Binary files a/braveheart_gui.fig and b/braveheart_gui.fig differ diff --git a/braveheart_gui.m b/braveheart_gui.m index da8d73b..4afb7f4 100644 --- a/braveheart_gui.m +++ b/braveheart_gui.m @@ -44,7 +44,7 @@ % Edit the above text to modify the response to help braveheart_gui -% Last Modified by GUIDE v2.5 20-Jun-2024 10:57:50 +% Last Modified by GUIDE v2.5 11-Sep-2024 13:17:11 % Update the current L&F for mac button issues... % Windows will use the normal Windows theme @@ -411,7 +411,7 @@ function about_help_button_Callback(hObject, eventdata, handles) {'\fontsize{14}\it\bf \color[rgb]{0.09,0.078,0.377}';... 'BRAVE\color[rgb]{0.89,0.016,0.016}H\fontsize{11}EART';... '\fontsize{8}\rm\color{black}(Beth Israel Analysis of Vectors of the Heart)';... -'Version 1.2.2' ; +'Version 1.3.0' ; ' ' ;... 'Copyright 2016-2024 Hans F. Stabeneau and Jonathan W. Waks' ;... ' ' ;... @@ -439,7 +439,7 @@ function about_help_button_Callback(hObject, eventdata, handles) {'\fontsize{18}\it\bf \color[rgb]{0.09,0.078,0.377}';... 'BRAVE\color[rgb]{0.89,0.016,0.016}H\fontsize{14}EART\fontsize{11}';... '\rm\color{black}(Beth Israel Analysis of Vectors of the Heart)';... -'Version 1.2.2' ; +'Version 1.3.0' ; ' ' ;... 'Copyright 2016-2024 Hans F. Stabeneau and Jonathan W. Waks' ;... ' ' ;... @@ -4807,7 +4807,7 @@ function threshold_fig_button_Callback(hObject, eventdata, handles) % handles structure with handles and user data (see GUIDATA) aps = pull_guiparams(hObject, eventdata, handles); -threshold_figure_gui(handles.vcg.VM, aps.pkthresh, aps.maxBPM, handles.vcg.hz); +threshold_figure_gui(handles.vcg, aps); % --- Executes on button press in pushbutton150. @@ -5449,3 +5449,12 @@ function BRAVEHEART_GUI_WindowKeyPressFcn(hObject, eventdata, handles) end + + +% --- Executes on button press in pkfilter_checkbox. +function pkfilter_checkbox_Callback(hObject, eventdata, handles) +% hObject handle to pkfilter_checkbox (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of pkfilter_checkbox diff --git a/braveheart_gui_windows.prj b/braveheart_gui_windows.prj index 59a2996..f2d6a64 100644 --- a/braveheart_gui_windows.prj +++ b/braveheart_gui_windows.prj @@ -8,13 +8,13 @@ ${PROJECT_ROOT}\braveheart_gui_windows_resources\icon_24.png ${PROJECT_ROOT}\braveheart_gui_windows_resources\icon_16.png - 1.2.2 + 1.3.0 Hans F. Stabenau and Jonathan W. Waks bravheart.ecg@gmail.com BRAVEHEART ECG/VCG Analysis Software -- GUI Version -Version 1.2.2 +Version 1.3.0 ${PROJECT_ROOT}\braveheart_splash.png \braveheart_gui\ @@ -50,7 +50,7 @@ Version 1.2.2 varargin input arguments Description BRAVEHEART ECG/VCG Analysis Software -- GUI Version -Version 1.2.2 +Version 1.3.0 @@ -95,6 +95,7 @@ Version 1.2.2 ${PROJECT_ROOT}\ecg_formats.csv ${PROJECT_ROOT}\elecpos.mat ${PROJECT_ROOT}\findHighpassLvl.m + ${PROJECT_ROOT}\findpeakswavelet.m ${PROJECT_ROOT}\getfiletype.m ${PROJECT_ROOT}\HRV_Calc.m ${PROJECT_ROOT}\leadidcodexyz.m diff --git a/braveheart_userguide.pdf b/braveheart_userguide.pdf index 08ceb04..e939ba6 100644 Binary files a/braveheart_userguide.pdf and b/braveheart_userguide.pdf differ diff --git a/findpeaksecg.m b/findpeaksecg.m index 36f35d3..5d8aa40 100644 --- a/findpeaksecg.m +++ b/findpeaksecg.m @@ -20,10 +20,26 @@ % This software is for research purposes only and is not intended to diagnose or treat any disease. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -function QRS = findpeaksecg(VM, maxbpm, freq, pkthresh) +function QRS = findpeaksecg(VM, maxbpm, freq, pkthresh, filter) + +if filter == 1 + [QRS, ~, ~] = findpeakswavelet(VM, maxbpm, freq, pkthresh, 'sym4'); +else + q = quantile(abs(VM), 100); + [~, QRS] = findpeaks(abs(VM), ... + 'MinPeakHeight', q(pkthresh), 'MinPeakDistance', round(freq*60/maxbpm)); +end + +% QRS = QRS'; +% +% QRS +%QRS2 + +% figure +% hold on +% plot(VM) +% scatter(QRS,VM(QRS)) +% scatter(QRS2,VM(QRS2)) -q = quantile(abs(VM), 100); -[~, QRS] = findpeaks(abs(VM), ... - 'MinPeakHeight', q(pkthresh), 'MinPeakDistance', round(freq*60/maxbpm)); end \ No newline at end of file diff --git a/findpeakswavelet.m b/findpeakswavelet.m new file mode 100644 index 0000000..c4ed873 --- /dev/null +++ b/findpeakswavelet.m @@ -0,0 +1,128 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% BRAVEHEART - Open source software for electrocardiographic and vectorcardiographic analysis +% findpeakswavelet.m -- Find R peaks using wavelet filtered signals +% Copyright 2016-2024 Hans F. Stabenau and Jonathan W. Waks +% +% Source code/executables: https://github.com/BIVectors/BRAVEHEART +% Contact: braveheart.ecg@gmail.com +% +% BRAVEHEART is free software: you can redistribute it and/or modify it under the terms of the GNU +% General Public License as published by the Free Software Foundation, either version 3 of the License, +% or (at your option) any later version. +% +% BRAVEHEART is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; +% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +% See the GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License along with this program. +% If not, see . +% +% This software is for research purposes only and is not intended to diagnose or treat any disease. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [QRS, QRS_wavelet, R] = findpeakswavelet(VM, maxbpm, freq, pkthresh, wavelet) + +% Max level of wavelet decomposition based on signal length +max_lvl = floor(log2(length(VM))); + +% Get freq bands for decomposition at each level (col 1 lower bound, +% col2 upper bound) +fc = zeros(max_lvl,2); + +for i = 1:max_lvl + fc(i,2) = freq / (2^i); + fc(i,1) = freq / (2^(i+1)); +end + +fc(max_lvl,1) = 0; + +% Want to keep frequencies 10-40 Hz which should capture most of QRS energy +% without T wave (<10 Hz) +Fa = 40; +Fb = 10; + +% Find closest value to 40 in upper bound of freq cutoffs +U = fc(:,2) - Fa; + +% Take smallest positive value so dont undercut the freq band +U(U<0) = nan; +[~,idxU] = min(U); + +% Find closest value to 10 in lower bound of freq cutoffs +L = fc(:,1) - Fb; + +% Take smallest postive value so dont undercut the freq band +L(L<0) = nan; +[~,idxL] = min(L); + +% Decompose signal to max level using MODWT +% Signal is mirrored to reduce edge effects +M = modwt(mirror(VM), max_lvl, wavelet); + +% Reconstruct using only the coefficients in levels between idxL and idxU +% and then square the resultant signal (R) +Mrecon = zeros(size(M)); +Mrecon(idxU:idxL,:) = M(idxU:idxL,:); +R = imodwt(Mrecon,wavelet); +R = abs(R).^2; +R = middlethird(R); + +% Find peaks in the reconstructed signal which will be estimates of the +% actual R peaks. + +q = quantile(R, 100); + +[~, QRS_wavelet] = findpeaks(R, ... + 'MinPeakHeight', q(pkthresh), 'MinPeakDistance', round(freq*60/maxbpm)); + +% Now have to use the peaks in the wavelet signal as an estimate of where +% to search for true peaks in VM since the max in the wavelet +% reconstruction is not always going to correlate exactly with the max in +% each QRS complex in VM + +% Open a window of 100 ms on each side of each peak detected in the wavelet +% reconstruction signal +win_len = 100; +w = round(win_len * freq / 1000); +QRS = zeros(1,length(QRS_wavelet)); + +for i = 1:length(QRS_wavelet) + ll = QRS_wavelet - w; + ul = QRS_wavelet + w; + + % Prevent indices from being out of bounds for QRS complexes near + % start/end of signal + + ll(ll<0) = 1; + ul(ul>length(VM)) = length(VM); + + % Just take maximum as should only be 1 peak + [~,QRS(i)] = max(VM(ll(i):ul(i))); + + % Correct indexing + QRS(i) = QRS(i) + ll(i) - 1; + +end + +debug = 0; +if debug + q1 = quantile(abs(VM), 100); + q2 = quantile(abs(R), 100); + + s = max(VM)/max(R); + R=R*s; + + figure + plot(R) + hold on + plot(VM) + scatter(QRS_wavelet,R(QRS_wavelet)) + scatter(QRS, VM(QRS)) + + line([0 length(VM)],[q(pkthresh)*s q(pkthresh)*s]) + %line([0 5000],[q2(95) q2(95)]) + +end + + +end \ No newline at end of file diff --git a/pull_guiparams.m b/pull_guiparams.m index 3443d2a..2ece553 100644 --- a/pull_guiparams.m +++ b/pull_guiparams.m @@ -43,6 +43,13 @@ % Peak Threshold for QRS detection aps.pkthresh = str2double(get(handles.pkthresh, 'String')); + +% Peak Threshold filtering + if get(handles.pkfilter_checkbox, 'Value') == 1 + aps.pkfilter = 1; + else + aps.pkfilter = 0; + end % Lowpass Wavelet filtering if get(handles.wavelet_filter_box, 'Value') == 1 diff --git a/push_guiparams.m b/push_guiparams.m index cdf1419..d4105e6 100644 --- a/push_guiparams.m +++ b/push_guiparams.m @@ -38,6 +38,9 @@ function push_guiparams(aps, hObject, eventdata, handles) % Peak Threshold for QRS detection set(handles.pkthresh, 'String', aps.pkthresh); + +% Peak Threshold filtering + set(handles.pkfilter_checkbox, 'Value', aps.pkfilter); % Lowpass Wavelet filtering set(handles.wavelet_filter_box, 'Value', aps.lowpass); diff --git a/scpopen.m b/scpopen.m index 58a499c..da87d08 100644 --- a/scpopen.m +++ b/scpopen.m @@ -142,6 +142,10 @@ Sect1Len = Sect1Len - 3 - len; %% [tag,len,Sect1Len], %% DEBUGGING information field = fread(fid,[1,len],'uchar'); + + % Disable if field is not read correctly + if ~isempty(field) + if tag == 0, ListOfRecommendedTags(ListOfRecommendedTags==tag)=[]; HDR.Patient.Name = char(field); %% LastName @@ -195,8 +199,8 @@ elseif tag == 9, HDR.Patient.Race = field; elseif tag == 10, -field, - if (field(1)~=0) + field, + if (field(1)~=0) HDR.Patient.Medication = field; else HDR.Patient.Medication.Code = {field(2:3)}; @@ -274,10 +278,12 @@ else fprintf(HDR.FILE.stderr,'Warning SCOPEN: unknown tag %i (section 1)\n',tag); end; + + end end; if ~isempty(ListOfRequiredTags) - fprintf(HDR.FILE.stderr,'Warning SCPOPEN: the following tags are required but missing in file %s\n',HDR.FileName); - disp(ListOfRequiredTags); + %fprintf(HDR.FILE.stderr,'Warning SCPOPEN: the following tags are required but missing in file %s\n',HDR.FileName); + %disp(ListOfRequiredTags); end; if ~isempty(ListOfRecommendedTags) %fprintf(HDR.FILE.stderr,'Warning SCPOPEN: the following tags are recommended but missing in file %s\n',HDR.FileName); diff --git a/test_braveheart.m b/test_braveheart.m index be4116f..80a2fdb 100644 --- a/test_braveheart.m +++ b/test_braveheart.m @@ -50,6 +50,7 @@ ap = Annoparams(); % Blank Annoparams class ap.maxBPM = 150; % Sets window for detecting R peaks ap.pkthresh = 95; % Percentile of ECG signal above which a peak can be found + ap.pkfilter = 0; % Additional filtering prior to peak detection on/off ap.lowpass = 1; % Low pass wavelet filter on/off ap.highpass = 1; % High pass wavelet filter on/off ap.wavelet_level_lowpass = 2; % Lvl 2 at 500 Hz and Lvl 3 at 997 Hz is 62.5 Hz low-pass filter diff --git a/threshold_figure_gui.m b/threshold_figure_gui.m index 01242f6..30328fc 100644 --- a/threshold_figure_gui.m +++ b/threshold_figure_gui.m @@ -21,10 +21,40 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -function threshold_figure_gui(vm, th, maxbpm, freq) +function threshold_figure_gui(vcg, aps) + +% Pass in VCG object (vcg) and Annoparams object (aps) +vm = vcg.VM; +th = aps.pkthresh; +maxbpm = aps.maxBPM; +freq = vcg.hz; + +% Detect peaks in VM signal without filtering +aps.pkfilter = 0; +QRS = vcg.peaks(aps); + +if aps.spike_removal + [vcg_noppm, ~] = vcg.remove_pacer_spikes(QRS, aps); + locs = vcg_noppm.peaks(aps); +else + locs = vcg.peaks(aps); +end + +% Detect peaks in VM signal with filtering +aps.pkfilter = 1; +QRS_filt = vcg.peaks(aps); + +if aps.spike_removal + [vcg_noppm_filt, ~] = vcg.remove_pacer_spikes(QRS_filt, aps); + locs_filt = vcg_noppm_filt.peaks(aps); +else + locs_filt = vcg.peaks(aps); + vcg_noppm_filt = vcg; +end + +% Plot original VM signal data up = vm; -%down = vm; q = quantile(abs(vm), 100); @@ -34,27 +64,15 @@ function threshold_figure_gui(vm, th, maxbpm, freq) set(gcf, 'Position', [100, 800, 1600, 300]) % set figure size plot(vm,'linewidth',1','color','black'); hold on -p_green = plot(up,'linewidth',2','color','green'); +p_green = plot(up,'linewidth',1.5','color','green'); p_line = line([0 length(vm)],[q(th) q(th)],'linestyle','--','color','blue'); xlim([0 length(vm)]); color = [{'red'} {'blue'}]; -color = repmat(color,1,40); +color = repmat(color,1,100); d = round(freq*60/maxbpm); -[pks, locs] = findpeaks((vm),'MinPeakHeight', q(th), 'MinPeakDistance', d); -% -% P =zeros(1,100); -% vm_trim = vm(501:4500); -% for i = 1:100 -% [~, QRS] = findpeaks((vm_trim),'MinPeakHeight', q(i), 'MinPeakDistance', d); -% P(i) = length(QRS); -% end -% -% pcut = 50; -% P(1:pcut) = []; - YL = ylim; for i = 1:length(locs) @@ -71,22 +89,75 @@ function threshold_figure_gui(vm, th, maxbpm, freq) p_area_red = area(A,[nan nan],'facecolor','red','edgecolor','none','facealpha',0.15); p_area_overlap = area(A,[nan nan],'facecolor','[0.8510, 0.7255, 0.8745]','edgecolor','none','facealpha',1); -%est_thresh = max(find(P==max(P(P==mode(P))))); -%title(sprintf('Initial R Peak Detection (Threshold %d%%, %d bpm) - Estimated Optimal Threshold = %d%%',th, maxbpm,est_thresh+pcut),'fontweight','bold','fontsize',14); title(sprintf('Initial R Peak Detection (Threshold %d%%, %d bpm)',th, maxbpm),'fontweight','bold','fontsize',14); p_green_txt = sprintf('VM Signal > %d %%ile', th); p_line_txt = sprintf('%d %%ile', th); -legend([p_green p_line p_dot_red p_area_red p_dot_blue p_area_blue, p_area_overlap], ... +legend([p_green p_line p_dot_red p_area_red p_dot_blue p_area_blue p_area_overlap], ... {p_green_txt, p_line_txt,'Odd R Peak','Window Around Odd R Peaks','Even R Peak','Window Around Even R Peaks', 'Window Overlap'}, ... 'location','eastoutside','fontsize',12); +% Increase font size on mac due to pc/mac font differences + if ismac + fontsize(gcf,scale=1.25) + end +hold off + + + +% Figure using filtered QRST complex +%vm = vcg_noppm.VM; +[~, QRS_filt_raw, vm_filt] = findpeakswavelet(vcg_noppm_filt.VM, maxbpm, freq, th, 'sym4'); + +% Normalize to voltage of VM unfiltered signal +N = max(vm)/max(vm_filt); +vm_filt = vm_filt * N; + +q_filt = quantile(abs(vm_filt), 100); + +up_filt = vm_filt; +up_filt(up_filt %d %%ile', th); +p_dot_black = scatter(nan,nan,'black'); + +legend([ecg_grey p_vm_filt p_green_filt p_line_filt p_dot_black p_dot_red p_area_red p_dot_blue p_area_blue p_area_overlap], ... + {'VM Signal', 'Filtered VM Signal', p_green_txt_filt, p_line_txt, 'Peaks of Filtered VM Signal', 'Odd R Peak','Window Around Odd R Peaks','Even R Peak','Window Around Even R Peaks', 'Window Overlap'}, ... + 'location','eastoutside','fontsize',12); + % Increase font size on mac due to pc/mac font differences if ismac fontsize(gcf,scale=1.25) end +hold off end \ No newline at end of file