function truelocs=find_peaks(y,varnoise)
%%%% zt-- threshold for local z-score
% tslope-- threshold for slope around peaks. The maximum slope has to be
%          greater than threshold to be qualified for being a peak
% y-- smoothed spectrum

% The return values are:
%	PEAKS: a list of peak locations, in ticks
%	LEFTMIN:	a list of the closest minimum to the left of each peak
%	RIGHTMIN: a list of the closest minimum to the right of each peak


scalexy=length(y)/(max(y)-min(y)); %  y value is multiplied by scalexy when 
                                   %  calculating angle related quantities
scalexy=scalexy/1.3;  % 1.3 is an approximted ratio figure scope
y=y*scalexy;
varnoise=varnoise*scalexy^2;
tslope=sqrt(3)/2;
tslope=1;
tslope=tan(40*pi/180);

%% Step 1: Find all local maxima and associated local minima.
w=y;
x = diff(w);		% first difference
xl = x(1:end-1);	% left side
xr = x(2:end);		% right side
neartop = (xl > 0 & xr <= 0) | (xl == 0 & xr < 0);
nearbot = (xl < 0 & xr >= 0) | (xl == 0 & xr > 0);
xx = 1:length(x);
nb = [1 xx(nearbot)+1];	% tentative minimum
nt = xx(neartop)+1;		% tentative maximum


%% Setp 2: Find out where the closest minima
% are to both left and right. These might not agree for adjacent peaks
% because of the presence of plateaus.
leftMin = [1];	% list of resolvable minima to left of peak.
rightMin = [1];	% similar list to right of peak.
lennt=length(nt);
for i = 2:lennt
   temp = nt(i-1):nt(i);
   leftMin(i) = max(temp(w(temp)==min(w(temp))));
   rightMin(i) = min(temp(w(temp)==min(w(temp))));
end
tt=[nt(lennt):length(y)];
rmin=min(tt(w(tt)==min(w(tt))));
rightMin = [rightMin(2:end) rmin];

plotflag=0;
%% display in figure
if plotflag==1    
figure;plot(nt,y(nt),'o',leftMin,y(leftMin),'d',rightMin,y(rightMin),'*');
hold on;plot(y);
end

%% Step 3: Coalesce rising plateaus
len=length(nt);
tol_flat=tslope/10;
count=1;
flagpeak=zeros(1,len);
for i=1:len
    curpos=nt(i);
    lmin=leftMin(i);
    rmin=rightMin(i);
    dl=(y(curpos)-y(lmin))/(curpos-lmin);
    dr=(y(curpos)-y(rmin))/(rmin-curpos);
    if dl>tol_flat & dr>tol_flat
        Peak(count)=curpos;
        count=count+1;
    end
end

leftMin = [1];	% list of resolvable minima to left of peak.
rightMin = [1];	% similar list to right of peak.
tol=(max(y)-min(y))/200;
for i = 2:length(Peak),
   temp = Peak(i-1):Peak(i);
   midmin=min(w(temp));
   ltol=min(tol,(w(Peak(i))-midmin)/4);
   rtol=min(tol,(w(Peak(i-1))-midmin)/4);
   leftMin(i) = max(temp(abs(w(temp)-midmin)<ltol));
   rightMin(i) = min(temp(abs(w(temp)-midmin)<rtol));
end

tt=[Peak(length(Peak)):length(w)];
rmin=min(tt(w(tt)==min(w(tt))));
rightMin = [rightMin(2:end) rmin];
%rightMin = [rightMin(2:end) length(x)];

plotflag=0;
%% display in figure
if plotflag==1    
figure;plot(Peak,y(Peak),'o',leftMin,y(leftMin),'d',rightMin,y(rightMin),'*');
hold on;plot(y);
end


%% Step 3: keep peak candidates whose slope is above threshold tslope
lenp=length(Peak);
leny=length(y);
count=1;
winwidth=2;
%dx=abs(x(2)-x(1));
%tslope=tslope*dx;
peak_locs=[];
for i=1:lenp
   lwin=max(Peak(i)-winwidth,1);
   rwin=min(Peak(i)+winwidth,leny);
   win1=[leftMin(i):rightMin(i)];
   win2=[lwin:rwin];
   wind=intersect(win1,win2);
   ywin=y(wind); 
   curdy=diff(ywin);
   if max(abs(curdy))>tslope
      peak_locs(count)=Peak(i);
      count=count+1;
   end
end

%% Step 4: Keep peaks whose amp^2 is above noise variance and
%%         mean(peaks)/20
lenp=length(peak_locs);
if lenp==0
    truelocs=[];
    return;
end
count=1;
for i=1:lenp
%    if y(peak_locs(i)) > max(sqrt(varnoise),mean(y(peak_locs)))
    if y(peak_locs(i)) > mean(y(peak_locs))/10;
        truelocs(count)=peak_locs(i);
        count=count+1;
    end
end

plotflag=0;
%% display in figure
if plotflag==1    
figure;plot(peak_locs,y(peak_locs),'o');
grid on;
%,'o',leftMin,y(leftMin),'d',rightMin,y(rightMin),'*');
hold on;plot(y);
    loc1=peak_locs(1);
    plot([loc1 loc1+50],[y(loc1) y(loc1)+50*tslope],'LineWidth',3);
    hold off;
end
%peak_amp=y(peak_locs)';
%





