Question about peak detection

I've written a script to find the percussive peaks of a provided icp signal, I'm very close but the markers I'm placing are off. I believe this is representative of my script marking them at the wrong index, just not sure how to fix this.

%% Plot ICP in Time Domain
clc; clear all; close all;

load icp;
t = (0:length(icp)-1)/fs;
figure('Color',[1 1 1]);
plot(t,icp);

%% Highpass, Remove DC Trend
fchp = 0.5;
Wp = fchp/(fs/2);
[B1,A1] = ellip(4,0.5,20,Wp,'high');
icphp = filtfilt(B1,A1,icp);
figure('Color',[1 1 1]);
h = plot(t,icphp);

figure('Color', [1 1 1]);
freqz(B1,A1,2^12,fs);
figure('Color',[1 1 1]);
zplane(B1,A1);


%% Lowpass, smooth signal
fclp = 3;
Wp = fclp/(fs/2);
[B2,A2] = ellip(2,0.5,20,Wp);
icplp = filtfilt(B2,A2,icphp);
figure('Color',[1 1 1]);
h = plot(t,icplp);

figure('Color', [1 1 1]);
freqz(B2,A2,2^12,fs);
figure('Color',[1 1 1]);
zplane(B2,A2);

%% Differentiator

icpd = diff([0;icphp]);
figure('Color',[1 1 1]);
h = plot(t,icpd);
set(h,'Color',[0.1 0.1 .1]);
set(h, 'Linewidth', 2);

%% Square Operator

icps = icpd.^2;
h = plot(t,icps);
set(h,'Color',[0.5 1 1]);
set(h,'Linewidth', 2);
xlabel('Time (s)');
ylabel('ICP');
box off; hold on;

%% Moving AVerage filter

B = 1/30*ones(30,1);
icpm = filtfilt(B,1,icps);
figure('Color', [1 1 1]);
h = plot(t,icpm);
set(h,'Color', [0.5 0.5 1]);
set(h,'Linewidth', 2);
title('Moving Average Filter');

%% Decision Logic 

x = icpm;
L = length(x);
id = find((x(1:L-2)<x(2:L-1))&(x(2:L-1)>x(3:L)))+1;
id2 = find(icp(id) > 1);
id2 = id(id2);
figure('Color',[1 1 1]);
h = plot(t, x, id/fs, x(id), 'r.');
set(h,'Markersize', 18);
title('Decision Logic');

%% Detection On Original Signal

figure('Color',[1 1 1]);

h = plot(t,icp,id2/fs,icp(id2),'r.');
set(h,'Markersize', 18);
set(h,'Linewidth', 2);

[Adam]

It's a bit difficult to understand what the outcome of all those operations will be on an unknown signal since you code starts by loading some data or other. Including the final plot so we can see what the problem actually is would help.

If you have the signal processing toolbos I would recommend using

doc findpeaks

though there is no reason, per se, why a simple 'home-made' method shouldn't also find the correct indices.

[Teapotimus]

This is just a very small portion of the whole signal. The script correctly places a marker on each wave, but the goal is to have it placed on the peak.

[darova]

[Adam]

Are the markers off by just one sample or more? And is it always the same amount? It looks like maybe more than 1 sample, but it can be hard to tell on steep sections of curve like that. If it is only 1 it may just be a matter of re-checking your code that calculates them to make sure the indices you are actually calculating are those of the peak point rather than the point just before the peak (e.g. using previous, current and next points to check for peaks, is your code returning the index of 'previous' rather than 'current', the centre point containing the actual peak?). There's often some +/- 1 logic that can easily get missed out or applied the wrong way with indexing things, especially when dividing by sample rates too.

[Teapotimus]

Should have attached this in the first place

Answer 1

[Cris LaPierre]

Do you have the Signal Processing toolbox? This can be accomplished with the findpeaks function.

[pk,loc]=findpeaks(icp,"MinPeakDistance",40);
t = (0:length(icp)-1)/fs;

figure
plot(t,icp)
hold on
plot(t(loc),pk,'r.')

Zoomed in: