creating a 1D vector

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

Answer 1

[Meg Noah]

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),'+');

Answer 2

[Meg Noah]

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

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),'+')

[Meg Noah]

Solution as function and a program to call it is attached.

If there is something else that this code needs to do, please explain.