ploteig.m

蒙托卡罗模拟奇异谱分析

function ploteig(E,V,c,opt)

% PLOTEIG - Plot an eigenspectrum with confidence limits
% Syntax: ploteig(E,V,c,opt)
%
% Inputs:  E - the basis associated with V.
%          V - an eigenspectrum (vector of length M).
%          c - an M by 2 matrix containing upper and
%              lower confidence limits.
%        opt - plotting option:
%              1 - plot rank on x-axis.
%              2 - plot dominant frequency on x-axis.
%              3 - plot rank on x-axis, and bolden the markers for 
%                  eigenvalues associated with nearly-sinusoidal eigenvectors. 
%              4 - plot dominant frequency on x-axis, and bolden the markers for 
%                  eigenvalues associated with nearly-sinusoidal eigenvectors. 
% 
%  Written by Eric Breitenberger.    Version date 5/24/96
%  Please send comments and suggestions to [e-mail]   

% PLOTEIG calls SHAPES, which is available in /pub/contrib/graphics/jnatools
% on the Mathworks FTP site (ftp.mathworks.com).

% Settings which can be varied according to taste:
sym_size=60;   % set relative size of symbols (smaller sym-size <=> bigger symbol)
aspect=1.3;    % set aspect ratio of plot
rftmin=.9;     % set minimum RFT value for bold markers.

N=length(V);

% Determine axis limits:
xmin=0;
xmax=N+1;
xrange=N+1;
ymin=min([c(:)' V]);
ymax=max([c(:)' V]);
ymin=10^(.1*floor(10*log10(ymin)));
ymax=10^(.1*ceil(10*log10(ymax)));
yrange=log10(ymax)-log10(ymin);

sym=eofsym(E, '1 or 0');  % Force EOFs to be 'symmetric' or 'anti-symmetric'
[mrft,f]=eoffreq(E);      % Find dominant frequencies of EOFs

if opt==2 | opt==4
  xmin=-.01;              % reset x-axis stuff
  xmax=.51;
  xrange=.52;
end

% build vector of symmetric EOF eigenvalues:
Vs=V(sym);
xs=find(sym);

% build vector of anti-symmetric EOF eigenvalues:
Vd=V(abs(sym-1));
xd=find(abs(sym-1));

% Set up for plotting - lots of monkeying around to get SHAPES to
% work with a semilog plot:
xscale=xrange/(sym_size*aspect);
yscale=log(10)*yrange/sym_size;

if opt==1 | opt==3  % Plot versus rank
  shapes(xd,log(Vd),'d',xscale,yscale);
  hold on
  shapes(xs,log(Vs),'s',xscale,yscale);
  hold off
  h=get(gca,'children');
  X=[];Y=[];
  for i=1:length(h), X(:,i)=get(h(i),'xdata')'; Y(:,i)=get(h(i),'ydata')'; end
  semilogy(X,exp(Y),'r'); 
  hold on
  cm=mean(c')';
  errorbar(1:N,cm,c(:,2)-cm,cm-c(:,1),'.y')
  h=get(gca,'children');
  set(h(1),'visible','off')
  axis([xmin xmax ymin ymax])

elseif opt==2 | opt==4 % Plot versus dom. frequency
  shapes(f(xd),log(Vd),'d',xscale,yscale);
  hold on
  shapes(f(xs),log(Vs),'s',xscale,yscale);
  hold off
  h=get(gca,'children');
  X=[];Y=[];
  for i=1:length(h), X(:,i)=get(h(i),'xdata')'; Y(:,i)=get(h(i),'ydata')'; end
  semilogy(X,exp(Y),'r'); 
  hold on
  cm=mean(c')';
  errorbar(f,cm,c(:,2)-cm,cm-c(:,1),'.y')
  h=get(gca,'children');
  set(h(1),'visible','off')
  axis([xmin xmax ymin ymax])
end

% Make bold markers if required:
% There are N+2 handles: the first 2 are for the errorbars; the next set
%   (ordinarily N/2 in number)) are for the anti-symmetric values
%   from largest to smallest; and the last set (also normally N/2
%   in number) are the symmetric eigenvalues from largest to smallest.
H=length(h);
Nsym=length(xs);
Nasym=length(xd);
if opt==3 | opt==4
  bold=find(mrft>rftmin);
  for b=1:length(bold)
    if sym(bold(b))==1
      index=find(xs==bold(b));
      set(h(H-Nsym+index),'linewidth',1.5,'color','c');
    else
      index=find(xd==bold(b));
      set(h(H-Nsym-Nasym+index),'linewidth',1.5,'color','c');
    end
  end
end

set(gca,'aspe',[aspect NaN])
hold off