%Reproducing sine wave interpolation
wave1 = (pi/2) - ([1:201]-1)*pi/50;
nwave1 = sin(wave1);
plot(nwave1,'-b');
hold on
wave2 = (wave1 + pi/4);
nwave2 = sin(wave2);
plot(nwave2,'-r')

array = [nwave1;nwave2]';
x=[1:201];
nx=x;
t=[1:2];
nt = [1.5];

%Plot linear interpolation of sin wave
data = array';
[x,t]=meshgrid(x,t);
[nx,nt]=meshgrid(nx,nt);
inter = interp2(x,t,data,nx,nt,'linear');
plot(inter,'-g')

%What the result should be...
wave3 = sin(wave1 + pi/8);


%dimensions for "data" array are [space,time] or [rows,cols]
nspace=size(array,1);
ntime=size(array,2);


tto=[1.5:1.5];
ttn=[1:2];

sp=201;

%Above removes spatial mean, better to do it this way...
%[m,n] = size(array);
%A = mean(array,2);	%temporal mean
%xc = array - repmat(A,1,n); %temporal anomalies for each row

%this method didn't seem to work...
% Apply the Hilbert transform.
%data1=[];
%for ii=1:sp
%  data1=[data1;hilbert(array(ii,:))];
%end

%this hilbert did work...
data1=hilbert(array);

%the variance-covariance matrix...where nspace is number of spatial points
covm=data1'*data1/(nspace-1);

[tmp Val Vec]=svd(covm);

Vec2=Vec'*data1';

% Interpolate each time series and reconstruct the data;
%  aa is amplitude ag is phase
recon1=0;

% choose number of modes
nmodes=ntime;
%sum(per(1:nmodes))

for ii=1:nmodes
  tmpt=Vec(:,ii);
  ag=unwrap(angle(tmpt));
  aa=abs(tmpt);
  ag2=interp1(ttn',ag,tto');
  aa2=interp1(ttn',aa,tto');
  tmp=aa2.*exp(i*ag2);
  recon1=recon1+tmp*Vec2(ii,:);
end

recon1=real(recon1);
recon1=recon1';

%xc4 = recon1 + A; %temporal anomalies for each row

xc4=recon1;

plot(xc4,'.-');ylim([-1 1])

%interpolation difference (ceof - linear)
%diff = check - linear;
%contourf(diff.');
return;
diff = abs(wave3-xc4);
hold off
plot(diff)

