MATLAB Answers

creating a 1D vector

48 views (last 30 days)
Ken Chi
Ken Chi on 10 Jan 2020
Commented: Meg Noah on 12 Jan 2020
I have been given a task to create a function. However i dont believe i have completed the whole task.
. It only gives me a figure with the plotted with all the maximum points. how would i alter this to function? Have been working on this code for many hours and any hints very much appreciated

Answers (2)

Meg Noah
Meg Noah on 10 Jan 2020
Edited: Meg Noah on 11 Jan 2020
Here's one of many possible solutions (edited 2020/01/11 for off-by-one error in indexing):
clc
close all
clear all
ECG = dlmread('ECG.csv');
nx = numel(ECG); x = 1:nx;
idx = find(ECG > 200);
mask = zeros(nx,1);
mask(idx) = 1;
labeled = bwlabel(mask);
stats = regionprops(labeled,ECG,'ALL');
ipeak = zeros(length(stats),1);
for iblob = 1:length(stats)
ipeak(iblob) = stats(iblob).PixelIdxList(1) - 1 + ...
find( ECG(stats(iblob).PixelIdxList) == stats(iblob).MaxIntensity);
end
figure('color','white','position',[70 100 600 300]);
plot(x,ECG);
hold on
title({'Heartbeat data'})
xlabel('x')
ylabel('ECG(x)')
plot(x(ipeak),ECG(ipeak),'+');

Meg Noah
Meg Noah on 11 Jan 2020
Here's a solution that has a step by step process.
clc
close all
clear all
%% Step 0
% load the data
ECG = dlmread('ECG.csv');
nx = numel(ECG); x = 1:nx;
Fs = 350; % [samples/s] frequency
T = 1/Fs; % [s] sampling period
N = 3000; % [samples] Length of signal
t = (0:N-1)*T; % [s] Time vector
% visualization
figure('color','white','position',[70 100 900 500]);
plot(1e3*t,ECG,'k','DisplayName','ECG Amplitude')
title({'Heartbeat data'})
xlabel('t (milliseconds)')
ylabel('X(t)')
hold on
%% Step 1
% First, determine points of the ECG scan exceed a given threshold.
threshold = 100;
idx = find(ECG > threshold);
% visualization
mask = nan(size(ECG)); mask(idx) = 1;
plot(1e3*t,mask.*ECG,'b','DisplayName',['ECG > ' num2str(threshold)])
%% Step 2
% Next, determine which of the data-points found in the previous step
% correspond to the same heartbeat. (Hint: data-points corresponding
% to one heartbeat will have consecutive indices, whereas transitions
% between heartbeats correspond to discontinuities in the index. diff
% function might help.
idxStartNewHeartBeat = [idx(1) idx(find(diff(idx)>1)+1)']';
% visualization
plot(1e3*t(idxStartNewHeartBeat),ECG(idxStartNewHeartBeat),'+r', ...
'DisplayName',['Unique Heartbeat ECG > ' num2str(threshold)]);
%% Step 3
% For each of the groups of data-points (corresponding to different
% heartbeats), find the location of the element with the maximum amplitude,
% and store this location in a vector.
mask(isnan(mask)) = 0;
labeledHeartbeats = bwlabel(mask);
stats = regionprops(labeledHeartbeats,ECG,'MaxIntensity'); %#ok<MRPBW>
maxAmplitudeValues = [stats.MaxIntensity];
idxMaxIndex = zeros(size(maxAmplitudeValues));
for iheartbeat = 1:length(maxAmplitudeValues)
idxMaxIndex(iheartbeat) = find(labeledHeartbeats == iheartbeat & ...
ECG == maxAmplitudeValues(iheartbeat));
end
% visualization
plot(1e3*t(idxMaxIndex),ECG(idxMaxIndex),'or', ...
'DisplayName',['Maximum in Heartbeat']);
%% Step 4
% The output should be a 1D vector containing the indices (one for each
% heart beat) corresponding to the local maxima.
% visualization
fprintf(1,'Indices of the heartbeat local maxima\n');
for iheartbeat = 1:length(idxMaxIndex)
fprintf(1,'%d\n',idxMaxIndex(iheartbeat));
ht = text(1e3*t(idxMaxIndex(iheartbeat)), ...
ECG(idxMaxIndex(iheartbeat))+30, ...
['idx = ' num2str(idxMaxIndex(iheartbeat))]);
set(ht,'rotation',90);
end
ylim([-100 500]);
hl = legend('location','northoutside');
It produces this image
HeartBeatMaxima.png
However, the original solution had an off by one error. But that doesn't even matter. The solution works.
mask = zeros(nx,1);
mask(idx) = 1;
labeled = bwlabel(mask);
stats = regionprops(labeled,ECG,'ALL');
ipeak = zeros(length(stats),1);
for iblob = 1:length(stats)
ipeak(iblob) = stats(iblob).PixelIdxList(1) - 1 + ...
find( ECG(stats(iblob).PixelIdxList) == stats(iblob).MaxIntensity);
end
figure('color','white','position',[70 100 600 300]);
plot(x,ECG);
hold on
title({'Heartbeat data'})
xlabel('x')
ylabel('ECG(x)')
plot(x(ipeak),ECG(ipeak),'+')
  1 Comment
Meg Noah
Meg Noah on 12 Jan 2020
Solution as function and a program to call it is attached.
If there is something else that this code needs to do, please explain.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by