最近看到一个有意思的求数组局部极小值,极大值的代码,贴出来分享一下,源代码是matlab版的,我用我的较为暴力的诸多for循环将其修改为C++版的,不得不感叹matlab在矩阵运算上确实是很方便的!
局部极大值和极小值都能够求得,以代码中 Arr[NUM] = { 1.31,2.52, 2.52, 6.84, 5.48, 2.10, 6.77, 6.77, 1.22, 1.35,9.02 }为例,可以得到局部极大值三个,6.84, 6.77,9.02. 局部极小值三个:1.31,2.10,1.22.
如有错误,请指出,谢谢!
Matlab版,源版
%EXTREMA Gets the global extrema points from a time series.
% [XMAX,IMAX,XMIN,IMIN] = EXTREMA(X) returns the global minima and maxima
% points of the vector X ignoring NaN's, where
% XMAX - maxima points in descending order
% IMAX - indexes of the XMAX
% XMIN - minima points in descending order
% IMIN - indexes of the XMIN
%
% DEFINITION (from http://en.wikipedia.org/wiki/Maxima_and_minima):
% In mathematics, maxima and minima, also known as extrema, are points in
% the domain of a function at which the function takes a largest value
% (maximum) or smallest value (minimum), either within a given
% neighbourhood (local extrema) or on the function domain in its entirety
% (global extrema).
%
% Example:
% x = 2*pi*linspace(-1,1);
% y = cos(x) - 0.5 + 0.5*rand(size(x)); y(40:45) = 1.85; y(50:53)=NaN;
% [ymax,imax,ymin,imin] = extrema(y);
% plot(x,y,x(imax),ymax,'g.',x(imin),ymin,'r.')
%
% See also EXTREMA2, MAX, MIN
% Written by
% Lic. on Physics Carlos Adri醤 Vargas Aguilera
% Physical Oceanography MS candidate
% UNIVERSIDAD DE GUADALAJARA
% Mexico, 2004
%
% nubeobscura@hotmail.com
% From : http://www.mathworks.com/matlabcentral/fileexchange
% File ID : 12275
% Submited at: 2006-09-14
% 2006-11-11 : English translation from spanish.
% 2006-11-17 : Accept NaN's.
% 2007-04-09 : Change name to MAXIMA, and definition added.
%**********尊重原创,以上是代码具体来源*****************
x=[1.31, 2.52, 2.52, 6.84, 5.48, 2.10, 6.77, 6.77, 1.22, 1.35];
xmax = [];
imax = [];
xmin = [];
imin = [];
% Vector input?
Nt = numel(x);
if Nt ~= length(x)
error('Entry must be a vector.')
end
% NaN's:
inan = find(isnan(x));
indx = 1:Nt;
if ~isempty(inan)
indx(inan) = [];
x(inan) = [];
Nt = length(x);
end
% Difference between subsequent elements:
dx = diff(x);
% Is an horizontal line?
if ~any(dx)
return
end
% Flat peaks? Put the middle element:
a = find(dx~=0); % Indexes where x changes
lm = find(diff(a)~=1) + 1; % Indexes where a do not changes
d = a(lm) - a(lm-1); % Number of elements in the flat peak
a(lm) = a(lm) - floor(d/2); % Save middle elements
a(end+1) = Nt;
% Peaks?
xa = x(a); % Serie without flat peaks
b = (diff(xa) > 0); % 1 => positive slopes (minima begin)
% 0 => negative slopes (maxima begin)
xb = diff(b); % -1 => maxima indexes (but one)
% +1 => minima indexes (but one)
imax = find(xb == -1) + 1; % maxima indexes
imin = find(xb == +1) + 1; % minima indexes
imax = a(imax);
imin = a(imin);
nmaxi = length(imax);
nmini = length(imin);
% Maximum or minumim on a flat peak at the ends?
if (nmaxi==0) && (nmini==0)
if x(1) > x(Nt)
xmax = x(1);
imax = indx(1);
xmin = x(Nt);
imin = indx(Nt);
elseif x(1) < x(Nt)
xmax = x(Nt);
imax = indx(Nt);
xmin = x(1);
imin = indx(1);
end
return
end
% Maximum or minumim at the ends?
if (nmaxi==0)
imax(1:2) = [1 Nt];
elseif (nmini==0)
imin(1:2) = [1 Nt];
else
if imax(1) < imin(1)
imin(2:nmini+1) = imin;
imin(1) = 1;
else
imax(2:nmaxi+1) = imax;
imax(1) = 1;
end
if imax(end) > imin(end)
imin(end+1) = Nt;
else
imax(end+1) = Nt;
end
end
xmax = x(imax);
xmin = x(imin);
% NaN's:
if ~isempty(inan)
imax = indx(imax);
imin = indx(imin);
end
% Same size as x:
imax = reshape(imax,size(xmax));
imin = reshape(imin,size(xmin));
% Descending order:
[temp,inmax] = sort(-xmax);
clear temp
xmax = xmax(inmax);
imax = imax(inmax);
[xmin,inmin] = sort(xmin);
imin = imin(inmin);
C++版,较挫,for循环多,但是也能用
#include <iostream>
#include <vector>
using namespace std;
#define NUM 11
void main()
{
float Arr[NUM] = { 1.31,2.52, 2.52, 6.84, 5.48, 2.10, 6.77, 6.77, 1.22, 1.35,9.02 };
int num = NUM;
float diff[NUM-1];
vector <int> indexA, indexLm;
//Difference between subsequent elements
int n = 0;
for (int i = 0; i < NUM - 1; i++)
{
diff[i] = Arr[i + 1] - Arr[i];
if (diff[i] != 0) indexA.push_back(i); //元素发生变化
}
// Flat peaks? Put the middle element
vector <int> diffIndexA;
for (int i = 0; i < indexA.size()-1; i++)
{
int tmpdiff = indexA.at(i + 1) - indexA.at(i);
if (tmpdiff != 1) indexLm.push_back(i + 1);
}
vector <int> d;
for (int i = 0; i < indexLm.size(); i++)
{
int index = indexLm.at(i);
int tmp = indexA.at(index)-indexA.at(index-1);
d.push_back(tmp);
}
for (int i = 0; i < d.size(); i++)
{
int lmValue = indexLm.at(i);
int dvalue = d.at(i) / 2;
int tmp = indexA.at(lmValue) - dvalue;
indexA.at(lmValue) = tmp;
}
indexA.push_back(num-1);
//Peak?
vector <float> ArrIndexA; //Seris without flat peaks
for (int i = 0; i < indexA.size(); i++)
{
int tmpIndex = indexA.at(i);
float value = Arr[tmpIndex];
ArrIndexA.push_back(value);
}
vector <int> indexB;
for (int i = 0; i < ArrIndexA.size()-1; i++)
{
float diff = ArrIndexA.at(i + 1) - ArrIndexA.at(i);
if (diff>0) indexB.push_back(1); //1 positive slopes (minima begin)
else indexB.push_back(0); //0 negative slopes(maxima begin)
}
vector <int> xb;
for (int i = 0; i < indexB.size()-1; i++)
{
int diff = indexB.at(i + 1) - indexB.at(i);
xb.push_back(diff); //-1 maxima indexes; 0 minima indexes
}
vector <int> imax, imin, max, min; //maxima indexes; minima indexes
for (int i = 0; i < xb.size(); i++)
{
if (xb.at(i) == -1) imax.push_back(i + 1);
if (xb.at(i) == 1) imin.push_back(i + 1);
}
for (int i = 0; i < imax.size(); i++)
{
int index = imax.at(i);
int value = indexA.at(index);
max.push_back(value);
}
for (int i = 0; i < imin.size(); i++)
{
int index = imin.at(i);
int value = indexA.at(index);
min.push_back(value);
}
int nmax = max.size();
int nmin = min.size();
//Maximum or minumin on a flat peak at the ends?
vector <float> xmin, xmax;
if (nmax == 0 && nmin == 0)
{
if (Arr[0]>Arr[num - 1])
{
xmax.push_back(Arr[0]);
max.push_back(0);
xmin.push_back(Arr[num - 1]);
min.push_back(num - 1);
}
else if (Arr[0] < Arr[num - 1])
{
xmax.push_back(Arr[num - 1]);
max.push_back(num - 1);
xmin.push_back(Arr[0]);
min.push_back(0);
}
return ;
}
//Maximum or minumin at the ends?
if (0 == nmax)
{
max.push_back(0);
max.push_back(num - 1);
}
else if (0 == nmin)
{
min.push_back(0);
min.push_back(num - 1);
}
else
{
if (max.at(0) < min.at(0))
{
vector <int> tmp;
tmp.swap(min);
min.push_back(0);
for (int i = 0; i < nmin; i++)
{
min.push_back(tmp[i]);
}
}
else
{
vector <int>tmp;
tmp.swap(max);
max.push_back(0);
for (int i = 0; i < nmax; i++)
{
max.push_back(tmp[i]);
}
}
if (max.back()>min.back())
{
min.push_back(num - 1);
}
else
{
max.push_back(num - 1);
}
}
for (int i = 0; i < max.size(); i++)
{
int index = max[i];
float value = Arr[index];
xmax.push_back(value);
}
for (int j = 0; j < min.size(); j++)
{
int index = min[j];
float value = Arr[index];
xmin.push_back(value);
}
//Descending Order, bubble sort
for (int j = 0; j < xmax.size()-1; j++)
{
for (int i = 0; i < xmax.size() - j -1; i++)
{
float tmp;
if (xmax[i] < xmax[i + 1])
{
tmp = xmax[i];
xmax[i] = xmax[i + 1];
xmax[i + 1] = tmp;
}
}
}
system("pause");
}