求助--两列信号的coherency的窗口数目问题? 附coherency程序(matlab)
请教各位达人，两列信号（时间序列）的coherency中的参数（取样间隔sampling interval、窗口数目number of windows）意义是什么，怎么选取？特别是窗口数目多大才算最好，为什么这样最好？ 默认的是8。%function = 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
phuyber@mit.edu]%phuyber@mit.edu
function = 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.
=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(').*sum('));
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=;
else,
fxc2=;
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=;
else,
fyc2=;
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));
=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
%derivedadjustment is made to c (i.e. cb).
xs =real(ifft(fx));
=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=;
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();
end;
axis tight; h=axis; axis(); h=axis;
w= NW/(dt*N); %half-bandwidth of the dpss
plot(,,'k');
for ds=min(s):2*w:max(s);
plot(,[.98 1.02],'k');
end;
%axis();
%set(gca,'xscale','log');
%plot(,h(3:4),'k--');
%plot(,h(3:4),'k--');
font(gca,12);
%h=text(1/70,.15,'(a)'); font(h,15);
end;
页:
[1]