# Help with loop script

2 views (last 30 days)
sam van Bohemen on 4 Dec 2019
Edited: Ridwan Alam on 8 Jan 2020
I have created a loop script to calculate peak amplitude for ECG data
x is the data file
TS is the total number of data points
L is the number of points sampled in each window
I have included an outlier code to remove peaks that are more than 3 std deviations from the mean.
I have written the code so I create a matrix that shows mean peak amplitude, peak std dev, number of peaks after outlier removal and number of peaks before outlier removal for each window.
In my code Bx represents the peaks detected for each window.
After I run my script Bx only shows the peaks detected for the last window
I want to alter my script so I can view Bx for all windows, Ideally in a matrix
My code is below:
Thanks!
xMatrix=reshape(x,L,TS/L);
%Filtering
[b,a]=butter(4,[0.0390625 0.46875],'bandpass'); %4th order bandpass butterworth 5-60Hz
xMatrix=filtfilt(b,a,xMatrix);
for N=(1:TS/L)
%Findpeaks
figure(1) % filtered ECG Plot x
plot(xMatrix(:,N));
[xpks,locs]=findpeaks(xMatrix(:,N),'MinPeakDistance',20,'MinPeakHeight',750);
Bx=rmoutliers(xpks,'mean');
[Bx,TFx]=rmoutliers(xpks);
xpeakmean=mean(Bx);
xpeakstd=std(Bx);
xspikes=numel(Bx);
xHR=numel(xpks);
findpeaks(xMatrix(:,N),'MinPeakDistance',20,'MinPeakHeight',750);
title('filtered ECG x')
xlabel('Samples')
ylabel('uV')
xlim([0 L])
ylim([-2000 2000])
%New Matrix
xpkMatrix(N,1)=xpeakmean;
xpkMatrix(N,2)=xpeakstd;
xpkMatrix(N,3)=xspikes;
xpkMatrix(N,4)=xHR;
end
##### 2 CommentsShowHide 1 older comment
sam van Bohemen on 8 Jan 2020
Hi Ridwan
I managed to figure it out
Thanks

Ridwan Alam on 5 Dec 2019
Edited: Ridwan Alam on 5 Dec 2019
If you can share x, I could test it properly. Otherwise, I think you just need 'hold on':
xMatrix=reshape(x,L,TS/L);
%Filtering
[b,a]=butter(4,[0.0390625 0.46875],'bandpass'); %4th order bandpass butterworth 5-60Hz
xMatrix=filtfilt(b,a,xMatrix);
figure(1) % filtered ECG Plot x
for N=(1:TS/L)
%Findpeaks
plot(xMatrix(:,N));hold on;
[xpks,locs]=findpeaks(xMatrix(:,N),'MinPeakDistance',20,'MinPeakHeight',750);
Bx=rmoutliers(xpks,'mean');
[Bx,TFx]=rmoutliers(xpks);
xpeakmean=mean(Bx);
xpeakstd=std(Bx);
xspikes=numel(Bx);
xHR=numel(xpks);
findpeaks(xMatrix(:,N),'MinPeakDistance',20,'MinPeakHeight',750);hold on;
title('filtered ECG x')
xlabel('Samples')
ylabel('uV')
xlim([0 L])
ylim([-2000 2000])
%New Matrix
xpkMatrix(N,1)=xpeakmean;
xpkMatrix(N,2)=xpeakstd;
xpkMatrix(N,3)=xspikes;
xpkMatrix(N,4)=xHR;
end
or new windows:
xMatrix=reshape(x,L,TS/L);
%Filtering
[b,a]=butter(4,[0.0390625 0.46875],'bandpass'); %4th order bandpass butterworth 5-60Hz
xMatrix=filtfilt(b,a,xMatrix);
% filtered ECG Plot x
for N=(1:TS/L)
%Findpeaks
figure;plot(xMatrix(:,N));hold on;
[xpks,locs]=findpeaks(xMatrix(:,N),'MinPeakDistance',20,'MinPeakHeight',750);
Bx=rmoutliers(xpks,'mean');
[Bx,TFx]=rmoutliers(xpks);
xpeakmean=mean(Bx);
xpeakstd=std(Bx);
xspikes=numel(Bx);
xHR=numel(xpks);
findpeaks(xMatrix(:,N),'MinPeakDistance',20,'MinPeakHeight',750);
title('filtered ECG x')
xlabel('Samples')
ylabel('uV')
xlim([0 L])
ylim([-2000 2000])
%New Matrix
xpkMatrix(N,1)=xpeakmean;
xpkMatrix(N,2)=xpeakstd;
xpkMatrix(N,3)=xspikes;
xpkMatrix(N,4)=xHR;
end
##### 2 CommentsShowHide 1 older comment
Ridwan Alam on 5 Dec 2019
Edited: Ridwan Alam on 8 Jan 2020
Hi Sam,
II tried to run the code without the Bx part, it did generate multiple plot windows showing both the signal and their peaks.
xnew = x(1:end-144);
L = 1000;
TS = length(xnew);
xMatrix=reshape(xnew,L,floor(TS/L));
%Filtering
[b,a]=butter(4,[0.0390625 0.46875],'bandpass'); %4th order bandpass butterworth 5-60Hz
xMatrix=filtfilt(b,a,xMatrix);
for i=1:TS/L
%Findpeaks
figure;plot(xMatrix(:,i));hold on;
findpeaks(xMatrix(:,i),'MinPeakDistance',20,'MinPeakHeight',750);
title('filtered ECG x')
xlabel('Samples')
ylabel('uV')
xlim([0 L])
ylim([-2000 2000])
% [xpks,locs]=findpeaks(xMatrix(:,i),'MinPeakDistance',20,'MinPeakHeight',750);
% Bx=rmoutliers(xpks,'mean');
% [Bx,TFx]=rmoutliers(xpks);
% xpeakmean=mean(Bx);
% xpeakstd=std(Bx);
% xspikes=numel(Bx);
% xHR=numel(xpks);
%New Matrix
% xpkMatrix(i,1)=xpeakmean;
% xpkMatrix(i,2)=xpeakstd;
% xpkMatrix(i,3)=xspikes;
% xpkMatrix(i,4)=xHR;
end
Please provide details about how far you got (make sure to show whole code) and what issues you are facing.