function [Pav,fav,kav]=freqav(P,f,k,nfb)
%function [Pav,fav,kav]=freqav(P,f,k,nfb)
%does frequency band averaging for runspec
  
if nfb==1
  kav=repmat(k,length(f),1);
  P(:,4) = sqrt(P(:,5));P(:,5)=180/pi*angle(P(:,3));
  Pav=P(:,1:5);
  fav=f;
else
  lgf=log10(f)*2^(nfb-1);
  mn=floor(lgf(1));
  mx=ceil(lgf(end));
  nbnd=mx-mn;
  nf=1;i1=0;i2=1;lisg=0;
  for i=1:nbnd
    isg=find(lgf>mn+i1 & lgf<mn+i2);
    if ~isempty(isg)&(length(isg)>=lisg)
      Pav(nf,1:3)=mean(P(isg,1:3),1);
      Pav(nf,4) = abs(Pav(nf,3))./sqrt(Pav(nf,1).*Pav(nf,2));
      Pav(nf,5) = (180/pi)*angle(Pav(nf,3));
      kav(nf)=k*length(isg);
      fav(nf)=mean(f(isg));
      nf=nf+1;i1=i2;i2=i2+1;
      lisg=length(isg);%this ensures monotonic behavior
    else
      i2=i2+1;
    end
  end
  kav=kav(:);
end

  