Code for peak detection

问题:

I want to calculate if a real-time signal pass some thresholds in the first step. In the first step, I want to detect if the real signal pass under those thresholds (in order to detect a peak in the signal). My Matlab code:

k=1;  
t = 1;
l=1;

for i =1:length(sm) //sm my signal.

    if (sm(i) > 0.25)
        first(k) = i;
        k = k+1;
        if (sm(i) > 0.5)

            second(t) = i;
            t =t +1;
            if (sm(i) > 0.75)

                third(l) = i;
                l = l+1;
            end
        end
    end
end

Example:

![enter image description here][1]

I want to calculate the times that the signal pass over and under the three thresholds 0.25, 0.5, 0.75 and to return those windows. Basically the three main peaks that I have in my example.

Basically what am I trying to do is to use fastsmooth function and findpeaks.

signalSmoothed = fastsmooth(sm,50); plot(signalSmoothed)
[max_pk1 max_pk2] = findpeaks(signalSmoothed);
find(max_pk1 >0.5)
inversex = 1.01*max(signalSmoothed) - signalSmoothed;
[min_pk1 min_pk2]  = findpeaks(inversex);
find(min_pk1 >0.5) 

Which are the heuristics in order to take only the desired peaks? Moreover the depticted image is an offline example. Generally I want to perform the technique online.

EDIT: I wrongfully defined as peak my desired curve result which is the whole wave not just the max value.

回答1:

Here is a solution to get the points where the signal sm passes the thresholds 0.25, 0.50 and 0.75. The points can be converted into windows inside the data-range and get stored in W. Finally we can plot them easily in the same figure. Note that we need to do some checks in the local function getwindows to handle special cases, for example when the window starts outside the data-range. The detection of windows inside another window is done in the getwindowsspecial-function.

Here is the code:

function peakwindow
% generate sample data
rng(7);
sm = 2*rand(1,25)-0.5;
sm = interp1(1:length(sm),sm,1:0.01:100*length(sm));

% get points
firstup    = find(diff(sm > 0.25)==1);
secondup   = find(diff(sm > 0.50)==1);
thirdup    = find(diff(sm > 0.75)==1);
firstdown  = find(diff(sm < 0.25)==1);
seconddown = find(diff(sm < 0.50)==1);
thirddown  = find(diff(sm < 0.75)==1);

% plot the result
figure; hold on;
plot(sm,'k')
plot(firstup,sm(firstup),'*')
plot(firstdown,sm(firstdown),'*')
plot(secondup,sm(secondup),'*')
plot(seconddown,sm(seconddown),'*')
plot(thirdup,sm(thirdup),'*')
plot(thirddown,sm(thirddown),'*')

% get windows
W1 = getwindows(firstup,firstdown);
W2 = getwindows(secondup,seconddown);
W3 = getwindows(thirdup,thirddown);

% get special window
WS = getwindowsspecial(W1,W3);

% plot windows
plotwindow(W1,0.25,'r');
plotwindow(W2,0.50,'r');
plotwindow(W3,0.75,'r');
plotwindow(WS,0,'b-');


function W = getwindows(up,down)
if length(up)>1 && length(down)>1 && up(1)>down(1)
    down(1)=[]; % handle case when window begins out of bounds left
end
if length(up)<1 || length(down)<1;
    W = []; % handle if no complete window present
else
    % concatenate and handle case when a window ends out of bounds right
    W = [up(1:length(down));down]';
end


function plotwindow(W,y,lspec)
for i = 1:size(W,1)
    plot(W(i,:),[y,y],lspec)
end


% get windows of U where there is a window of H inside
function W = getwindowsspecial(U,H)
W = []; % empty matrix to begin with
for i = 1:size(U,1) % for all windows in U
    if any(H(:,1)>=U(i,1) & H(:,1)<=U(i,2))
        W = [W;U(i,:)]; % add window
    end
end

This is the result:

---

To see that the handling works properly, we can plot the result when initialized with rng(3):

Note that the window of 0.25 and 0.50 would start out of bounds left and therefore is not present in the plotted windows.