[ideas wanted]: How to find "peaks" of sloping signal?

Shoutout to you all math wizards! I need to detect peaks of (possibly) sloping signals. i.e sometimes first "peak" is real peak (findpeaks), but other times "peak" is in upward sloping portion of signal. See example figure below: the blue triangle is what findpeaks finds as first peak. But I need to detect that "peak" before it (see red arrow).

All ideas are much appreciated!

data is in attachment (data.mat)

Answer 1

[Star Strider]

The online ability to run scripts seems to be broken.

I used MATLAB Online to devise this —

% LD = load('data.mat')
file = websave('data.mat','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1809778/data.mat');
LD = load(file);
data = LD.data;
x = linspace(0, 5.5, numel(data));                      % Supply Missing Coordinate
Fs = 1/x(2)

[FTs1,Fv] = FFT1(data,x);

figure
plot(Fv, abs(FTs1)*2)
xlim([0 3])
title('Fourier Transform (Informs Filter Design)')

data_filt = highpass(data, 0.0425, Fs, ImpulseResponse='iir');

[pksf,locsf] = findpeaks(data_filt);
[pks,locs] = findpeaks(data);

locs = [locsf(1) locs]
pks = data(locs)

figure
plot(x,  data)
hold on
plot(x(locs), data(locs), 'rv' )
hold off

figure
plot(x,  data_filt)
title('Filtered ‘data’')

function [FTs1,Fv] = FFT1(s,t)
% Arguments:
%    s: Signal Vector Or Matrix
%    t: Associated Time Vector 
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    Fv = Fv(:);
    FTs1 = FTs(Iv,:);
end

.

[Sven Larsen]

elegant solution! thank you so much :)

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

[Sven Larsen]

by the way, why have you hard coded value 0.0425 on :

data_filt = highpass(data, 0.0425, Fs, ImpulseResponse='iir');

how did you come up to that particular number? im applying your solution to my overall program, and need to understand the reasoning for that value.

[Star Strider]

That is the frequency that works best. (Signal processinig is frequently experimental.) I used the Fourier transform result to provide the first approximation, then went from there until I got a result that works.

Note that it is based on the ‘x’ vector. Basing it on the index vector is an optiion, and in that instance the sampling frequency is 1, and the cutoff frequency would have to be calculated differently, and likely experimentally, although scaling it as a function of the calculated sampling frequency here (0.0425/Fs) would likely work.

.

[Star Strider]

The option to run code now works again!

Anyway, when you mentioned that you want to use my code on other data sets, I decided to see if I could come up with something better that would not involve filtering, because I am not certain that the same approach will always work without re-designing thee filter each time.

This initially ‘truncates’ the ‘data’ vector to include only a small part near the top (defined by everything greater than or equal to 0.99 multiplied by the lowest ‘valley’ point returned by findpeaks, assuming all the valleys are positive). Treating the end points of that vector as ‘valleys’ for the time being, the code then creates points conneecting adjacent ‘valley’ points, connects them with a liine with the same number of ponts, and then subtracts that from the ‘data’ values in the same range, effectively ‘normalising’ them to a constant zero-valuee y- line reference. It then returns the locations of the peeaks in that range, refereenced to the original index vector. That will reeturn the best estimiate for the first ‘peak’, and imprecise estimates for the others, beeecause by design it ‘tilts’ the ‘data’ vector in each of those ranges. The code then finds the precise peaks of the original ‘data’ vector, and concatenates the first location to them to create the final complete ‘locs’ vector.

Now, the only part of this that resists a programmatic approach is the last step, because that requires manually creating the final ‘locs’ vector. (I am not certain that the ‘missing’ peak is always the first, as opposed to being either the first or the last. If it is always the first, then the last step can be automated as I did here.)

Try this —

file = 'data.mat';
LD = load(file);
data = LD.data;
x = linspace(0, 5.5, numel(data));                      % Supply Missing Coordinate
[datamin,datamax] = bounds(data);
datarng  = datamax-datamin;

[vlys,vlocs] = findpeaks(-data);

vlys_ofst = min(-vlys);
Lv = data >= (vlys_ofst * 0.99);                        % Logical Vector
% QQ = nnz(Lv)

x_tr = x(Lv);                                           % ‘Truncated’ ‘x’
data_tr = data(Lv);                                     % ‘Truncated’ ‘data’

[vlys_tr,vlocs_tr] = findpeaks(-data_tr);               % Find ‘Valleys6
vlocs_tr = [1 vlocs_tr numel(data_tr)];

for k = 1:numel(vlocs_tr)-1                             % Normalisation Loop
    idxrng = vlocs_tr(k) : (vlocs_tr(k+1));
    x_bl{k} = x_tr(idxrng);
    baseline = interp1(x_tr([min(idxrng) max(idxrng)]), data_tr([min(idxrng) max(idxrng)]), x_bl{k});
    data_bl{k} = data_tr(idxrng) - baseline;
    [pks{k},loc] = findpeaks(data_bl{k});
    locsc{k} = loc + idxrng(1);
end    

pks

pks = 1x5 cell array
    {[0.6507 1.1108]}    {[0.7926]}    {[0.7462]}    {[0.7361]}    {[0.8566]}

locs = [locsc{:}] + find(Lv, 1, 'first')                    % Normalised ‘locs’ With Respect To The First ‘Lv’ Index

locs = 1×6

46391 70337 109976 144683 180437 215416

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x_locs = x(locs)

x_locs = 1×6

0.9894 1.5001 2.3455 3.0857 3.8483 4.5943

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure                      % This Plot Demonstrates How The Normalisation Loop Works And Can Be Deleted
hold on
for k = 1:numel(x_bl)
    % data_bl{k}
    plot(x_bl{k}, data_bl{k})
end
hold off
xlim([0 6])
title('Individual Normalised Segments Of The ‘data’ Vector')

[fpks, plocs] = findpeaks(data)

fpks = 1×5

89.3733 89.5787 89.4956 89.2324 88.6728

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

plocs = 1×5

72491 110271 144190 179670 213368

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = [locs(1) plocs]

locs = 1×6

46391 72491 110271 144190 179670 213368

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
plot(x, data, '-b')
hold on
plot(x(locs), data(locs)+0.25, 'vr', MarkerFaceColor='r')
hold off
title('Final Plot Showing All Defined Peaks')

.