# creating a 1D vector

48 views (last 30 days)
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

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
nx = numel(ECG); x = 1:nx;
idx = find(ECG > 200);
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 on 11 Jan 2020
Here's a solution that has a step by step process.
clc
close all
clear all
%% Step 0
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
%% 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.
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
However, the original solution had an off by one error. But that doesn't even matter. The solution works.
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 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.