% SVD.M  2/01/2007
% This is supposed to be exacly like eof.m, except it uses the SVD of the data matrix instead of the 
% eigenanalysis
% This is an example of scaling EOFs
% Start by constructing a simple data set
% D.L. Hartmann ATM S 552
norm = inputwd('If you want to normalize input data, type 1, otherwise 0 (0): ',0)
xnoise = inputwd('Set the noise level: 0.5 is low 2.0 is high (1.0): ',1.0)
N=100
M=20
per=12.
yes='yes';
no='no '
ans='no '
if norm == 1
ans=yes;
end

% In MATLAB the second index is the row index, which goes across the columns and 
% gives the row data.
% Notice that the amplitude is made linearly dependent on the structure dimension
% index j=1:M
for i=1:N
 for j=1:M
  x(j,i)=5.0*(j/M)*(sin(2*pi*(j-1)/(M-1))*sin(2*pi*i/per)+xnoise*randn(1));
 end
end
xm=x;

% remove mean and possibly standard deviation from data time series
for j=1:M
  xbar=mean(x(j,:));
   if norm == 1
    sdev=std(x(j,:));
    else
    sdev=1;
   end
 for i=1:N
  xm(j,i)=(x(j,i)-xbar)/sdev;
 end
end

% calculate covariance matrix
%C=xm*xm'/N;
%C

% Do eigenanalysis
%[E,D]=eig(C);

% Do SVD of data matrix
[u s v] = svd(xm);
sizeu= size(u)
sizev= size(v)
sizes = size(s)

%calculate autocorrelation
for j=1:M
alpha(j,:)=auto(xm(j,:),1);
end
alpha
% the temporal autocorrelation is different for each spatial grid point,
% a problem.  Let's kluge this by using the mean of the mean and the maximum
% autocorrelation.  We'll take the absolute value before averaging, so that
% random negative correlations can't cancel out positive ones.  This is
% conservative.
ralph=(mean(abs(alpha(:,2)))+max(alpha(:,2)))/2.
tau=-1.0/log(ralph)
% Let's use Taylor/Leith for N*/N ---
Nstar=min(N/(2*tau),N)
N
% Let's use North et al for the 95% limits on the eigenvalues
factor=sqrt(2./Nstar)

% Plot Eigenvalues with North et al Whiskers
L=diag(s,0);
sizeL = size(L)
L=L.^2/N;
LE=L*factor;
for j=1:M
 ind(j)=j;
end
errorbar(ind,L,LE)
title(['Eigenvalue Spectrum:  Norm=  ' ans])
xlabel('index')
ylabel('Eigenvalue')

pause
% Plot eigenvectors
E=u;
plot(E(:,3),'g-.')
hold on
plot(E(:,2),'b--')
hold on
plot(E(:,1),'r')
hold off
MP=M+1
axis([0 MP -0.60 0.60])
%set(gca,'DataAspectRatio',[1 22 1])
title(['First Three Eigenvectors:  Norm=  ' ans])
xlabel('Structure dimension')
ylabel('Normalized EOF Amplitude')
pause

% Project eigenvectors onto original data to get PCs
Z=E'*x;
plot(Z(3,:)', 'g-.')
hold on
plot(Z(2,:)', 'b--')
hold on
plot(Z(1,:)','Color', 'r')
hold off
title(['First Three Principle Components:  Norm=  ' ans])
xlabel('Sampling dimension')
ylabel('PC Amplitude')
pause
%Normalize principle components time series to have zero mean and unit variance
ZM=Z;
for i=1:N
 for j=1:M
   zbar=mean(Z(j,:));
   sdev=std(Z(j,:));
  ZM(j,i)=(Z(j,i)-zbar)/sdev;
 end
end

% remove mean from original, unnormalized data time series
for i=1:N
 for j=1:M
  xbar=mean(x(j,:)); 
  xm(j,i)=x(j,i)-xbar;
 end
end

% Project data onto normalized PC time series

EM=xm*ZM'/N;

%plot first three dimensional eigenvector regressions
plot(EM(:,3),'g-.')
hold on
plot(EM(:,2),'b--')
hold on
plot(EM(:,1),'Color','r')
hold off
xmax=max(abs(EM(:,1)))*1.2;
axis([0 MP -xmax xmax])

title(['Regressions of First Three Normalized PCs on Unnormalized Data:  Norm=  ' ans])
xlabel('Structure dimension')
ylabel('Amplitude in Physical Units')