Finding correct MinPeakDistance using for loop

Question

Joel 2023-4-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1952163-finding-correct-minpeakdistance-using-for-loop

编辑： chicken vector 2023-4-24

I am trying to find the peaks when a signal (EFFEKT) exceeds a threshold (Abonnemang)

EFFEKT is an 8760x1 double, containing a value for every hour of the year.

Some peaks involves several values (see figure), thus I use this loop to make every peak only contain the biggest value (max_val). t0_pos1 is the starting point; when a peak crosses the threshold on the way up and t0_neg1 on the way down. The if-statement makes sure that the number of crossings is correlated with the number of peaks.

However the if statements is not activated for certain threshold values like in this example.

Any ideas on how I can overcome this?

load EFFEKT.mat
Abonnemang=21200;
dt = datetime(2021,1,1:1:365);
et = hours(0:23)';
t = repmat(dt,length(et),1) + repmat(et,1,length(dt));
t=t(:);
for x=1:numel(EFFEKT)-2 %REMOVE NEIGHBOURING PEAKS
    [max_val,locs] = findpeaks(EFFEKT,'MinPeakHeight',Abonnemang,'MinPeakDistance',x); % find the peak and its x axis location
    [t0_pos1,s0_pos1,t0_neg1,s0_neg1]= crossing_V7(EFFEKT,t,Abonnemang,'linear'); % positive (pos) and negative (neg) slope crossing points
    peak_width = t0_neg1- t0_pos1;
    if length(max_val)==length(t0_pos1)
        break
    end
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

chicken vector 2023-4-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1952163-finding-correct-minpeakdistance-using-for-loop#answer_1221998

编辑：chicken vector 2023-4-24

在 MATLAB Online 中打开

Your idea of increasing 'MinPeakDistance' to adjust the peaks found is not consistent.

For example, if you iterate manually, you can see that when length(max_val) is approaching length(t0_pos1), you are already starting to lose some valuable information as shown below in green.

There was a related question yesterday. You can use the method explained there, which is originally from this question to divide peaks in group.

clear;clc;
load EFFEKT.mat
Abonnemang = 21200;
t = datetime(2021,1,1,1:length(EFFEKT),0,0)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Credits @KSSV:
ad = EFFEKT' > Abonnemang;
ad0 = zeros(size(ad)) ;
ad0(logical(ad)) = find(ad) ;
ii = zeros(size(ad0));
jj = ad0 > 0;
ii(strfind([0,jj(:)'],[0 1])) = 1;
idx = cumsum(ii).*jj;
out = accumarray( idx(jj)',ad0(jj)',[],@(x){x'});
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nPeaks = length(out);
peaks = zeros(1,nPeaks);
peakTime = NaT(1,nPeaks);
firstCrossing = NaT(1,nPeaks);
lastCrossing = NaT(1,nPeaks);
for j = 1 : nPeaks
    [maxValue, maxIdx] = max(EFFEKT(out{j}));
    peaks(j) = maxValue;
    peakTime(j) = t(out{j}(maxIdx));
    firstCrossing(j) = t(out{j}(1));
    lastCrossing(j) = t(out{j}(end));
end
figure;
tiledlayout(1,2);
nexttile
hold on;
plot(t,EFFEKT);
yline(Abonnemang);
scatter(peakTime,peaks,'r')
title('Peaks')
set(gca,'FontSize',24)
peakIdx = 34;
plotIdx = out{peakIdx}(1)-2:out{peakIdx}(end)+2;
nexttile
hold on;
plot(t(plotIdx),EFFEKT(plotIdx));
yline(Abonnemang);
scatter(peakTime(peakIdx),peaks(peakIdx),'r')
xline(firstCrossing(peakIdx),'k')
xline(lastCrossing(peakIdx),'k')
title('Crossing')
set(gca,'FontSize',24)

Result:

On a sidenote, I suggest you to pay attention on what really depends on the looping variable, in your case x.

For example you are performing the 8760 times the command:

[t0_pos1,s0_pos1,t0_neg1,s0_neg1]= crossing_V7(EFFEKT,t,Abonnemang,'linear'); % positive (pos) and negative (neg) slope crossing points
peak_width = t0_neg1- t0_pos1;

But this does not depend on x, therefore it may be performed only once prior to the loop, largely increasing the speed of your code.