1. 程式人生 > 實用技巧 >peakdet: Peak detection using MATLAB 峰識別 峰面積計算 peak area 相關matlab基本詳解

peakdet: Peak detection using MATLAB 峰識別 峰面積計算 peak area 相關matlab基本詳解

轉自http://www.billauer.co.il/peakdet.html

peakdet: Peak detection using MATLAB

peakdet:使用MATLAB的峰值檢測

Here's a problem I encounter in several fields: Find the local maxima and minima in some noisy signal, which typically looks like the following graph:

在幾個領域,我遇到一個問題:在一些噪聲訊號中找到區間的最大值和最小值,它通常看起來像下面的圖:

The local maxima and minima are plotted as red and green stars on the graph. To the eye it's so obvious where they are, but making a computer find them can turn out tricky.

該區間的最大值和最小值被繪製成圖上的紅色和綠色的星星。用眼睛很明顯的發現他們在哪裡,但是用電腦程式設計發現他們是棘手的。

Let's start with whatnot to do: Using the well-known zero-derivate method. Due to the noise, which is always there in real-life signals, accidental zero-crossings of the first derivate occur, yielding false detections. The typical solution is to smooth the curve with some low-pass filter, usually killing the original signal at the same time. The result is usually that the algorithm goes horribly wrong where it's so obvious to the eye.

讓我們這樣開始:使用著名的零導數的方法。由於噪聲,它總是在現實生活中的訊號的過零點,偶然的一階導數時,產生錯誤檢測。典型的解決方案是光滑的曲線具有一定的低通濾波器,通常造成原始訊號的同時。結果通常是,演算法錯的離譜,眼睛看到的確實那麼明顯。

In many cases, we don't really care about maxima and minima in the mathematical sense. We can see the peaks and valleys, and we want the computer to find them. This is what "peakdet" does.

在許多情況下,我們真的不在乎的最大值和最小值的數學意義。我們可以看到山峰和山谷,我們想讓電腦找到他們。這就是“peakdet”。

The trick here is to realize, that a peak is the highest point betweem "valleys". What makes a peak is the fact that there are lower points around it. This strategy is adopted by "peakdet": Look for the highest point, around which there are points lower by X on both sides.

找到峰的技巧是這樣,一個峰是谷之間的最高點。峰是點,在該點的周圍有較它低的點。這種策略是由“peakdet”採用:尋找最高點,在x軸上,該點兩邊的點的y值都小於該點的y值。

Let's see an example: First, let's create the graph shown in the figure above:

例如:首先,讓我們建立上圖中所示的圖:

>>t=0:0.001:10;
>>x=0.3*sin(t) + sin(1.3*t) + 0.9*sin(4.2*t) + 0.02*randn(1, 10001);
>>figure; plot(x);

 t=0:0.001:10;
 x=0.3*sin(t) + sin(1.3*t) + 0.9*sin(4.2*t) + 0.02*randn(1, 10001);
 figure; plot(x);
>> [maxtab, mintab] = peakdet(x, 0.5);
>> hold on; plot(mintab(:,1), mintab(:,2), 'g*');
>> plot(maxtab(:,1), maxtab(:,2), 'r*');
[maxtab, mintab] = peakdet(x, 0.5);
hold on; plot(mintab(:,1), mintab(:,2), 'g*');
plot(maxtab(:,1), maxtab(:,2), 'r*');

Note the call to peakdet(): The first argument is the vector to examine, and the second is the peak threshold: We require a difference of at least 0.5 between a peak and its surrounding in order to declare it as a peak. Same goes with valleys.

注意呼叫peakdet():第一個引數是待檢測的向量,第二個是閾值:我們需要至少0.5的閥值,在峰和它其周圍以宣佈它作為一個峰值之間的差。谷同理。

The returned vectors "maxtab" and "mintab" contain the peak and valley points, as evident by their plots (note the colors).

返回向量”maxtab”和“Mintab”包含峰谷點,明顯的情節(注意顏色)。

The vector's X-axis values can be passed as a third argument (thanks to Sven Billiet for his contribution on this), in which case peakdet() returns these values instead of indices, as shown in the following example:

向量的x軸的座標值可以作為第三個引數傳遞(感謝Sven Billiet的貢獻在這),在這種情況下返回這些值代替peakdet()指標,如下面的示例所示:

>> figure; plot(t,x);
>> [maxtab, mintab] = peakdet(x, 0.5, t);
figure; plot(t,x);
[maxtab, mintab] = peakdet(x, 0.5, t);

And from here we continue like before, but note that the X axis represents "t" and not indices.

從這裡我們繼續像以前一樣,但請注意,X軸代表“t”而不是指數。

>> hold on; plot(mintab(:,1), mintab(:,2), 'g*');
>> plot(maxtab(:,1), maxtab(:,2), 'r*');
hold on; plot(mintab(:,1), mintab(:,2), 'g*');
plot(maxtab(:,1), maxtab(:,2), 'r*');
%合併後
t=0:0.001:10;
x=0.3*sin(t) + sin(1.3*t) + 0.9*sin(4.2*t) + 0.02*randn(1, 10001);
figure; 
plot(x); [maxtab, mintab] = peakdet(x, 0.5);
hold on; 
plot(mintab(:,1), mintab(:,2), 'g*');
plot(maxtab(:,1), maxtab(:,2), 'r*');
As for the implementation of this function: The work is done with a for-loop, which is considered lousy practice in MATLAB. Since I've never needed this function for anything else than pretty short vectors (< 100000 points), I also never bothered to try speeding it up. Compiling to MEX is a direct solution. I'm not sure if it's possible to vectorize this algorithm in MATLAB. I'll be glad to hear suggestions.
為實現此項功能:工作完成一個迴圈,這被認為是在MATLAB的糟糕的實踐。因為我從來都不需要這種功能的任何其他比很短向量(<100000分),我也從來沒有試著加快速度。編寫Mex是一個直接的解決方案。我不確定它是否可以量化,該演算法在MATLAB。我很樂意聽的建議。

A final note: If you happen to prefer Python,you could try this(someone has been kind enough to convert this function). There are also aversion in Cby Hong Xu and aversion in FORTRAN 90by Brian McNoldy. I haven't verified any of these.

最後一點:如果你碰巧喜歡Python,你可以試一下這個(有人已經足夠將這個功能)。也有C版本(Hong Xu)和FORTRAN 90(Brian McNoldy)。我還沒有證實任何這些。

And here is the function. Copy and save it as 'peakdet.m'. It's released to the public domain:

函式matlab程式碼如下。複製並儲存為“peakdet.m”。

注意:使用matlab中的file-new-script 建立m檔案,否則會出現亂碼,不能正常執行。

function [maxtab, mintab]=peakdet(v, delta, x)
%PEAKDET Detect peaks in a vector
%        [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
%        maxima and minima ("peaks") in the vector V.
%        MAXTAB and MINTAB consists of two columns. Column 1
%        contains indices in V, and column 2 the found values.
%      
%        With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
%        in MAXTAB and MINTAB are replaced with the corresponding
%        X-values.
%
%        A point is considered a maximum peak if it has the maximal
%        value, and was preceded (to the left) by a value lower by
%        DELTA.

% Eli Billauer, 3.4.05 (Explicitly not copyrighted).
% This function is released to the public domain; Any use is allowed.

maxtab = [];
mintab = [];

v = v(:); % Just in case this wasn't a proper vector

if nargin < 3
  x = (1:length(v))';
else 
  x = x(:);
  if length(v)~= length(x)
    error('Input vectors v and x must have same length');
  end
end
  
if (length(delta(:)))>1
  error('Input argument DELTA must be a scalar');
end

if delta <= 0
  error('Input argument DELTA must be positive');
end

mn = Inf; mx = -Inf;
mnpos = NaN; mxpos = NaN;

lookformax = 1;

for i=1:length(v)
  this = v(i);
  if this > mx, mx = this; mxpos = x(i); end
  if this < mn, mn = this; mnpos = x(i); end
  
  if lookformax
    if this < mx-delta
      maxtab = [maxtab ; mxpos mx];
      mn = this; mnpos = x(i);
      lookformax = 0;
    end  
  else
    if this > mn+delta
      mintab = [mintab ; mnpos mn];
      mx = this; mxpos = x(i);
      lookformax = 1;
    end
  end
end

function [maxtab, mintab]=peakdet(v, delta, x) 函式功能介紹

function [maxtab, mintab]=peakdet(v, delta, x)
%v,x軸順序取值(12、...10001)時,對應的y軸的值
%delta 差值 △X的意思
%x,x軸順序取值(12、...10001%maxtab 二維陣列,類似於座標圖上的點座標集合,每個點point(x,y)表示,它有n個point組成
%maxtab 第一列:順序存放點的x座標值,第二列:順序存放與x軸點對應的的y座標值,
%maxtab 每一行存放一個座標點point(x,y)
%PEAKDET Detect peaks in a vector
%        [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
%        maxima and minima ("peaks") in the vector V.
%        MAXTAB and MINTAB consists of two columns. Column 1
%        contains indices in V, and column 2 the found values.
%      
%        With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
%        in MAXTAB and MINTAB are replaced with the corresponding
%        X-values.
%
%        A point is considered a maximum peak if it has the maximal
%        value, and was preceded (to the left) by a value lower by
%        DELTA.

% Eli Billauer, 3.4.05 (Explicitly not copyrighted).
% This function is released to the public domain; Any use is allowed.

maxtab = [];%存放峰點座標
mintab = [];%存放谷點座標

v = v(:); % Just in case this wasn't a proper vector

if nargin < 3 %nargin:輸入變數的個數
  x = (1:length(v))'; %1到10001
else 
  x = x(:);
  if length(v)~= length(x) %~= 不等於 
    error('Input vectors v and x must have same length');
  end
end
  
if (length(delta(:)))>1 %輸入delta只能是一個值  並且大於等於0
  error('Input argument DELTA must be a scalar'); %scalar 標量
end

if delta <= 0
  error('Input argument DELTA must be positive');
end

mn = Inf; mx = -Inf;  %inf 當算出的結果大於某個數(這個數很大,
                      %比如10的很多次方),則MATLAB認為就是無窮大了,並返回 inf
mnpos = NaN; mxpos = NaN; % NaN:Not  a  Number 類似初始化賦值

lookformax = 1; %取峰和谷的標記  1:取峰,0取谷

for i=1:length(v) %length(v) v中值的個數 10001個
  this = v(i);%v(i)代表第i個點y軸值
  if this > mx, mx = this; mxpos = x(i); end %x(i)代表第i個點x軸值
  if this < mn, mn = this; mnpos = x(i); end
  
  if lookformax
    if this < mx-delta %y軸值的比較  滿足該條件的點在峰的右邊
      maxtab = [maxtab ; mxpos mx]; %給maxtab(二維向量)賦值,x軸值mxpos,y軸值mx
      mn = this; mnpos = x(i);
      lookformax = 0; %已取到峰,下一次取谷
    end  
  else
    if this > mn+delta %y軸值的比較  滿足該條件的點在谷的右邊
      mintab = [mintab ; mnpos mn];
      mx = this; mxpos = x(i);
      lookformax = 1; %已取到谷,下一次取峰
    end
  end
end

[轉]https://blog.csdn.net/wyx100/article/details/43083921

http://www.billauer.co.il/peakdet.html

https://github.com/xuphys/peakdetect/blob/master/peakdetect.c

http://terpconnect.umd.edu/~toh/spectrum/ipeak.html

https://terpconnect.umd.edu/~toh/spectrum/Introduction.html