eoffreq.m

蒙托卡罗模拟奇异谱分析

function [mrft,f]=eoffreq(E,n)
% EOFFREQ - find the dominant frequencies of EOFs.
% Syntax: [mrft,f]=eoffreq(E); [mrft,f]=eoffreq(E,n);
%
% Given the eigenvector matrix E, eoffreq computes a 
% normalized reduced Fourier transform (RFT) for each
% EOF. The maximum values of the RFTs are returned in 
% 'mrft', and the frequencies at which the maxima 
% occur are returned in 'f'.
% 
% The second (optional) argument 'n' gives the number
% of frequencies which will be checked - the default
% is 500. The frequencies are determined by splitting
% the Nyquist interval into n sections.   
%
% Written by Eric Breitenberger.      Version 3/25/97
% Please send comments and suggestions to [e-mail]       
%

if nargin==1, n=500; end

s=eofsym(E); % eigenvector symmetries

[M,K]=size(E);

% center the eigenfunctions:
E=E-ones(M,1)*mean(E);

f=linspace(0,.5,n)';
F=zeros(n,K);
j=1:M;
j2=j-(M+1)/2;

Cc=f*j2;
Cs=sin(2*pi*Cc);
Cc=cos(2*pi*Cc);
r2c=sum((Cc.^2)')/M;
r2s=1-r2c;
r2s(1)=eps; % to avoid divide by zero.
if r2s(n)==0, r2s(n)=eps; end
for k=1:K
  if s(k)==1 
    F(:,k)=abs(Cc*E(:,k)).^2./r2c';
  elseif s(k)==0    
    F(:,k)=abs(Cs*E(:,k)).^2./r2s';
  elseif s(k)==-1
    F(:,k)=abs(Cc*E(:,k)).^2+abs(Cs*E(:,k)).^2;
  else
    error('Elements of s must be 1, 0, or -1.')
  end
end

[mrft,f]=max(F);
mrft=mrft/M;
f=.5*(f-1)/(n-1);