findpeaks.m

此为频谱分析工具箱

function [peak,pmax,sy,ny,sdy,ndy,sddy,nddy,pmin,pmaxsav] = findpeaks(x,ys,thr, pos,width,last)
% findpeaks : Quick search of peaks in a signal.
%Syntax: [peakmat,pmax,sy,ny,sdy,ndy,sddy,nddy,pmin] = findpeaks(x,y,{threshold=1, pos,width,mode ='silent'}) or findpeaks(y)
%
% This function searches peaks into an (x,y) signal, using a specified
% 'threshold' in noise units, usually around 1 ; few peaks for a high value.
% 'peakmat' is an array of rows for each maximum identified
%          [ index left_width right_width  max_pos Intensity Width ].
% 'pmax' is a Dirac type peaks signal.
% Averaged signals (including derivatives) and noises are also returned.
% opt : pos, width = part of 'y' to be scanned for peaks around pos.
%       mode is 'silent', 'normal' or 'verbose' (do plot).

% Author:  EF <[e-mail]>
% Description: quick search of peaks in a signal.

% Part of 'Spectral tools'. E.Farhi. 07/96 rev 12/97
% uses : diff.m
%	 smooth.m 
%	 interp.m

if ~exist('last')
	last = 'normal';
end

if isempty(last)
	last = 'normal';
end

if strcmp(last,'silent')
	tmp = 0;
elseif strcmp(last,'normal')
	tmp = 1;
else
	tmp = 2;
end

if (nargin < 1)
	error('usage : [peaks,pmax,sy,ny,sdy,ndy,sddy,nddy] = findpeaks(x,y,{threshold=1}) or findpeaks(y)');
end

if (tmp>0)
	fprintf(1,'Peak Searching\n');
end

if (nargin == 1)
	ys=x;
	x=1:length(x);
end

if ~exist('thr') 
	thr = [];
end
if isempty(thr)
	thr = 1;
end

if ~exist('pos') | ~exist('width')
	pos = [];
	width = [];
end
if isempty(pos) |  isempty(width)
	nd = length(x);
	pos = x(ceil(nd/2));
	width = abs(x(nd)-x(1));
end

if thr == 0, thr = 0.0001; end

y=ys(:);

nd = length(y);
dx = x(2) - x(1);

%if any(diff(x) ~= dx) % step is not constant, need to interpolate
%	nx = linspace(x(1), x(nd), round((x(nd)-x(1))/dx));
%	y = interp1(x,y,nx);
%	x=nx;
%	plot(x,y)
%end

t=find( (x>=(pos-width)) & (x<=(pos+width)) );		% index zone around pos
t = t(:);
y = y(t); y=y(:); 
x = x(t); x=x(:); 
nd = length(y);	

sy = smooth (y); 
dy = diff (sy); dy(nd) = dy(nd-1); 
sdy = smooth (dy); 
ddy = diff (sdy); ddy(nd) = ddy(nd-1);
sddy = smooth (ddy);

t=(abs(y) <= median(abs(y)));		% noises computation
t = t(:);
nt = sum(t);
ny = std(t.*y - t.*sy);
t=(abs(dy) <= median(abs(dy)));
ndy = std(t.*dy - t.*sdy);
t=(abs(ddy) <= median(abs(ddy)));
nddy = std(t.*ddy - t.*sddy);

t = dy(1:(nd-1)).*dy(2:nd); 	% for 1st derivative sign changes.
t(nd) = sdy(nd)^2;
t = t(:);

pmax = (abs(dy) <= ndy/thr) | (t <= (ndy/thr)^2);	% extrema

pmin = ( pmax & (sddy > nddy*thr));
pmax = ( pmax & (sddy < -nddy*thr));

% now alterning maxima and minima ...

pmaxsav = pmax;
peaklist=find(pmax);

if (tmp > 0)
	fprintf(1,'* analysing extrema (2nd order interpolation), %i possible peaks.',length(peaklist));
	fprintf(1,'\nnumber        indexes     max_pos     width   intensity (exp int)');
end

index =1;		% index in original peak list (from pmax)
ip = 1;		% index in final peak table
peak = [ 0 0 0 0 0 ];
pmax = 0*y;

while (index <= length(peaklist))	 % scanning all peaks in pmax
	peakindex = peaklist(index);
	t = find(pmin(1:peakindex));	% search min before
	if (~isempty(t))
		minbefore = t(length(t));
	else
		minbefore = 1;
	end
	t = find(pmin(peakindex:nd));	% search min after
	if (~isempty(t))
		minafter = t(1) + peakindex -1;
	else
		minafter = nd;
	end
	t = minbefore:minafter;
	[dummy,t] = max(y(t));	% find real max beween 2 min.
	if (isempty(t))
		index = index+1;			% no max found
		% trying next in pmax
	else
	max_pos = t(1) + minbefore -1;			% max found
	t = max(y(minbefore),y(minafter));
	if ((y(max_pos) - t) <= ny)	% this is a noise peak.
		if (y(minbefore) > y(minafter)) % equivalent minima
			pmin(minbefore) = 0;
		else			% higher minimum removed
			pmin(minafter) = 0;
		end
		if (minafter == nd)
			break;			% finish
		else
			if (minbefore == 1)
				index = index +1;
			end  % trying again

		end
	else				% real peak
		if (~isempty(find((peak(:,1) - max_pos) == 0)))
			index = index + 1;		% already found
		else
		peak(ip,1) = max_pos;
		t = min(y(minbefore),y(minafter));
		t = 0.1*y(max_pos) + 0.9*t;	% 10 percent hight for peak
		halfafter = find(y(max_pos:minafter) < t);
		if (isempty(halfafter))
			halfafter = minafter;
		else
			halfafter = halfafter(1) + max_pos -1;
		end
		halfbefore = find(y(minbefore:max_pos) < t);
		if (isempty(halfbefore))
			halfbefore = minbefore;
		else
			halfbefore = halfbefore(length(halfbefore)) + minbefore -1 ;
		end
		peak(ip, 2) = max_pos - halfbefore;
		peak(ip, 3) = halfafter - max_pos;
		t = halfbefore:halfafter;
%disp([ max_pos length(t) ])
% might be done with polyfit but this is better...
		[ smoothed, p] = interpsp(x(t),y(t),x(max_pos),max(ceil(length(t)/2),3),2);
%		p = polyfit(x(t),y(t),2);

		S = x(max_pos);
		W = abs(x(halfafter) - x(halfbefore));
		In = y(max_pos);

		if (length(p) == 3)
			a=p(1); b=p(2); c=p(3);	% 2nd order interpolation
			delta = (b*b - 4*a*c);
			if ((a < 0) & (delta >0))
		  	  S=-b/2/a;
			  W=abs(sqrt(delta)/a);
			  In=a*S*S+b*S+c;
			else
				t = minbefore:minafter;
				[ smoothed, p] = interpsp(x(t),y(t),x(max_pos),max(ceil(length(t)/2),3),2);
%				p = polyfit(x(t),y(t),2);

				if (length(p) == 3)
					a=p(1); b=p(2); c=p(3);	% 2nd order interpolation
					delta = (b*b - 4*a*c);
delta = -1;
					if ((a < 0) & (delta >0))
		 		 	  S=-b/2/a;
					  W=abs(sqrt(delta)/a);
					  In=a*S*S+b*S+c;
					else
					  p = [];
					end
				end
			end
		else			% interpolation not performed
			p = [];
		end
		end
		if S > max(x) | S < min(x) | W > abs(max(x)-min(x))
			p = [];
		end
		if isempty(p)			% interpolation not performed
			S = x(max_pos);
			W = abs(x(halfafter) - x(halfbefore));
			In = y(max_pos);
		end

		t = [ minbefore minafter ];
		wx = find(x > S+W); 
		if ~isempty(wx), t = [ t wx(1)]; end
		wx = find(x >= S-W); 
		if ~isempty(wx), t = [ t wx(1) ]; end
		t = abs(t - max_pos) * 3;
		t = max(t);
		t = (max_pos -t):(max_pos+t);
		t = t(find(t>=1 & t<= nd));
		t = min(y(t)); % compared with nearer min peak
		peak(ip,4) = S;
		peak(ip,5) = In - t;
		peak(ip,6) = W;
		if (tmp > 0)
			fprintf(1,'\nPeak %3i : (%3i-%3i-%3i) S=%7.2f W=%7.2f In=%7.2f (%7.2f)', ip, peak(ip, 2), peak(ip, 1), peak(ip, 3), S, W, In, y(max_pos));
		end
		ip = ip + 1;
		index = index + 1;
		pmax(max_pos) = 1;
		end % from if (~isempty(find((peak(:,1) - max_pos) == 0)))
	end
	end
end
if tmp
	fprintf(1,'\n');
end
	
if (tmp > 1)
	fprintf(1,'\n%i peaks identified\n', ip-1);
	clf;
  	title('Peaks position in signal');
	plot(x,y,x,pmax.*y);
end

if (size(ys,2)==1)
	sy = sy';
	pmax = pmax';
	sdy = sdy';
	sddy = sddy';
end