求助--两列信号的coherency的窗口数目问题? 附coherency程序(matlab)

shizhu018

请教各位达人，两列信号（时间序列）的coherency中的参数（取样间隔sampling interval、窗口数目number of windows）意义是什么，怎么选取？特别是窗口数目多大才算最好，为什么这样最好？ 默认的是8。
%function [s, c, ph, ci, phi] = cmtmPH(x,y,dt,NW,confn,qplot);
%
%simple multi-taper method coherence estimate and Monte Carlo
%confidence estimation procedure.
%
% Inputs:
%      x        - Input data vector 1.
%      y        - Input data vector 2.
%      dt       - sampling interval
%      NW       - number of windows to use (8 is often a good choice)
%      confn    - number of Monte Carlo iterations to use (default 50)
%                 in computing approximate 95\% confidence intervals. 0 or
%                 omit to forgo confidence intervals.  In some cases as many
%                 as 200 iterations are required for results to converge.
%      qplot    - plot the results, 1= Yes, 0 or omit = No.  The upper
%                 tickmarks indicate the bandwidth of the coherence and
%                 phase estimates.
%
% Outputs:
%      s       - frequency
%      c       - coherence
%      ph      - phase
%      ci      - 95% coherence confidence level
%      phi     - 95% phase confidence interval
%
%
%Peter Huybers
%MIT, 2003
[email=%phuyber@mit.edu]%phuyber@mit.edu[/email]
function [s, c, ph, ci, phi] = cmtmPH(x,y,dt,NW,confn,qplot);
x   = x(:)-mean(x); y=y(:)-mean(y);
%check input
if nargin<6, qplot=0;  end;
if nargin<5, confn=0; end;
if length(confn)==0, confn=0; end;
if nargin<4, NW=8;     end;
if length(NW)==0, NW=8;end;
if nargin<3, dt=1;     end;
if length(dt)==0, dt=1;end;
%define some parameters
N   = length(x);
k   = min(round(2*NW),N);
k   = max(k-1,1);
s   = (0:1/(N*dt):1/dt-1/(N*dt))';
pls=2:(N+1)/2+1;
if rem(length(y),2)==1; pls=pls(1:end-1); end;
%Compute the discrete prolate spheroidal sequences, requires the
%spectral analysis toolbox.
[E,V]=dpss(N,NW,k);
%Compute the windowed DFTs.
fkx=fft(E(:,1:k).*x(:,ones(1,k)),N);
fky=fft(E(:,1:k).*y(:,ones(1,k)),N);
if confn>0,
  Fx=mean(abs(fkx)')';
  Fy=mean(abs(fky)')';
end;
%Normalize the energy of each tapered transform
for ct=1:length(E(1,:));
  fkx(:,ct)=fkx(:,ct)/std(fkx(:,ct));
  fky(:,ct)=fky(:,ct)/std(fky(:,ct));
end;
%Compute coherence
Cxy= sum( (fkx.*conj(fky))' );
c  = abs(Cxy)./sqrt(sum([fkx.*conj(fkx)]').*sum([fky.*conj(fky)]'));
ph = angle(Cxy)*180/pi;
%Estimate confidence intervals
if confn>0,
  disp('Estimating confidence intervals');
  for iter=1:confn;
    if rem(iter,10)==0, disp(['iteration: ',num2str(iter)]); end;
    %Make filters for phase uncertainty estimates
    fxc=abs(fft(x));
    pl=find(fxc(pls)<2.5*std(fxc))+1;
    fxc2=exp(polyval(polyfit(log(s(pl)),log(fxc(pl)),2),log(s(pls))));
    if rem(length(x),2)==1;
      fxc2=[0; fxc2; flipud(fxc2)];
    else,
      fxc2=[0; fxc2; flipud(fxc2(2:end))];
    end;
      
    fyc=abs(fft(y));
    pl=find(fyc(pls)<2.5*std(fyc))+1;
    fyc2=exp(polyval(polyfit(log(s(pl)),log(fyc(pl)),2),log(s(pls))));
    if rem(length(y),2)==1;
      fyc2=[0; fyc2; flipud(fyc2)];
    else,
      fyc2=[0; fyc2; flipud(fyc2(2:end))];
    end;
    %coherence
    ys=randn(size(y));
    ys=ys-mean(ys); ys=ys/std(ys);
    ys=real(ifft(fft(ys).*fyc2));   
    xs=randn(size(x));
    xs=xs-mean(xs); xs=xs/std(xs);
    xs =real(ifft(fft(xs).*fxc2));
    [si, ci(iter,:), dum]=cmtmPH(xs,ys,dt,NW);
    %phase
    fy=fft(randn(size(y))+rand).*fft(y);
    fx=fft(randn(size(x))+rand).*fft(x);
    fy=fy/sum(abs(fy));
    fx=fx/sum(abs(fx));
    cb=c-(1-c).^3; pl=find(cb<0); cb(pl)=0;
    ys =real( ifft(fy.*sqrt(1-cb'.^2)));   
    ys =ys+real(ifft(fx.*cb'));    %NOTE: coherence is a biassed estimator.
                                   %Runs with known signals indicated a bias
                                   %of +.3 for incoherent processes, with
                                   %the biass tapering off toward higher
                                   %true coherence.  Thus an empiracally
                                   %derived  adjustment is made to c (i.e. cb).
    xs =real(ifft(fx));
    [si, ciph(iter,:), phi(iter,:)]=cmtmPH(xs,ys,dt,NW);
  end;
  %sorting and averaging to determine confidence levels.
  pl=round(.95*iter);  
  ci=sort(ci);   
  ci=mean(ci(pl,:));
  ci=mean(ci)*ones(size(ci));
  pl=round(.975*iter);  
  phi=sort(phi);  
  phi=[phi(pl,:); -phi(iter-pl+1,:)];
  phi=[mean(phi); -mean(phi)];
  phi(1,1:3)=mean(phi(1,1:3));
  phi(2,1:3)=mean(phi(2,1:3));
  phi(1,end-2:end)=mean(phi(1,end-2:end));
  phi(2,end-2:end)=mean(phi(2,end-2:end));
  temp=conv(phi(1,:),[.25 .5 .25]);
  phi(1,4:end-3)=temp(5:end-4);
  temp=conv(phi(2,:),[.25 .5 .25]);
  phi(2,4:end-3)=temp(5:end-4);
end;
%Cut to one-sided funtions
c = c(pls);
s = s(pls);
ph=ph(pls);
%plotting
if qplot==1,
  %coherence
  figure(gcf);
  %subplot(211);
  hold on;
  plot(s,c);
  h=ylabel('coherence'); font(h,15);
  if confn>0;
    plot(si,ci*ones(size(si)),'k:');
    pl=find(c>ci(1));
    %title([num2str(100*length(pl)/length(c),2),'% of estimates above 95% confidence level']);
  end;
  axis tight; h=axis; axis([h(1:2) 0 1.025]); h=axis;
  w  = NW/(dt*N);   %half-bandwidth of the dpss
  plot([s(1) h(2)],[1.02 1.02],'k');
  for ds=min(s):2*w:max(s);
    plot([ds ds],[.98 1.02],'k');
  end;
  %axis([h(1) 4 h(3) 1.025]);
  %set(gca,'xscale','log');
  %plot([1/1.5 1/1.5],h(3:4),'k--');
  %plot([1/10 1/10],h(3:4),'k--');
  font(gca,12);
  %h=text(1/70,.15,'(a)'); font(h,15);
end;