function peakf=PsdPeak2(X,interval,leadset)
% This function finds the frequency which has peak PSD in most
% leads. X is the input PSD matrix, and each column is for a different
% lead. The function sorts the input matrix and returns the index which
% most of the leads agree to have the max power.
peakf=-1;
if any(leadset)
tot=sum(leadset);
lint=length(interval);
[Xs,Xi]=sort(X);
Y=X(interval,:);
[m,n]=size(Y);
[Ys,Yi]=sort(Y);
Yi=Yi'+interval(1)-1;
hn=zeros(4,lint);
hn_shiftl=zeros(4,lint);
hn_shiftr=zeros(4,lint);
for i=1:4
if leadset(i)
hn(i,:)=hist(Yi(i,end-2:end),interval);
hn_shiftl(i,:)=[hn(i,2:end),0];
hn_shiftr(i,:)=[0 hn(i,1:end-1)];
end
end
skewed_hn=hn+hn_shiftl+hn_shiftr;
skewed_hn=skewed_hn>0;
test=ones(1,lint);
for i=1:4
if leadset(i)
test=test.*skewed_hn(i,:);
end
end
windowed_test=filter([1 1 1],1,test);
[tm,ti]=max(windowed_test);
if tm>0
if sum(windowed_test==tm)==1
peakf=interval(ti-1);
else
k=interval(find(windowed_test==tm)-1);
s=zeros(1,length(k));
for i=1:length(k);
s(i)=sum(any(Xi(end-2:end,leadset)==k(i))+any(Xi(end-2:end,leadset)==k(i)-1)+any(Xi(end-2:end,leadset)==k(i)+1));
end
[sm,si]=max(s);
if sum(s==sm)==1
peakf=k(si);
end
if all(diff(k)==1)
peakf=k(floor((length(k)+1)/2));
end
end
end
end