% This program computes the 1-D wavelet transform of a REAL time series x. 
%          format: [p,wx]=wavelet(x)
% p is the period scale in log2(p) (so as to properly label the "image" plot,
%    to be done by "wlplot")
% y is the cosine (symmetric, real) part of the wavelet spectrum.
%
% .......................................................................
% Reference: Zhou YH, Zheng DW, Liao XH, Wavelet analysis of 
%            interannual LOD, AAM and ENSO: 1997~98 El Niņo 
%            and 1998~99 La Niņa signals, Journal of Geodesy, Vol.75, 
%            Page 164-168, 2001.
% .......................................................................
%
function [p,wx]=wf(x)
n=length(x);
x=detrend(x)';
mu=pi*sqrt(2/log(2));
% NO NEED TO pad zeros to power of 2 to facilitate convolution later...
% nt=2^(nextpow2(nx))
% xx=[xx,zeros(1,nt-nx)];
% set up the (vertical) period scale...
dt=input('sampling interval = ');
pmin=input('minimum display period = ');
pmax=input('maximum dispaly period = ');
%np=ceil(log2(pmax/pmin))*4+1;
np=ceil(log(pmax/pmin)/log(2))*4+1;
%
fact=ceil(100/np);
if fact>1
   np=np*fact-fact+1;
end
%p=logspace(pmin,pmax,np);
%
pmin=log(pmin)/log(10);
pmax=log(pmax)/log(10);
p = (10).^ [pmin+(0:np-2)*(pmax-pmin)/(np-1), pmax];
%
% construct the wavelet and convolve...
a=mu*p/(2*pi);
y0=zeros(np,n);
for m=1:np
   nw2=ceil(5*p(m)/dt);
   nw=2*nw2+1;
%
%   tw=linspace(-nw2*dt,nw2*dt,nw);
%
   ddd1=-nw2*dt;
   ddd2=nw2*dt;
   tw=[ddd1+(0:nw-2)*(ddd2-ddd1)/(nw-1) ddd2];
%
   w=exp(-tw.^2/(2*a(m)^2)).*cos(mu*tw/a(m))/sqrt(2*pi*a(m));
%
   num=n+nw-1;
   nn=2^(nextpow2(num));
   yy=real(ifft(fft(x,nn).*fft(w,nn)));
%
%   yy=conv(x,w);
   y0(m,:)=yy(nw2+1:nw2+n)*dt/sqrt(a(m));
end
% interpolate the p scale to about 100 points...
%fact=ceil(100/np);
%if fact>1
%   np1=np*fact-fact+1;
%
%   zzad=interp(y0(:,1),fact);
%   [zm,zn]=size(zzad);
%
%   wx=zeros(np1,n);
%
%    wx=zeros(zm,n);
%
%   for i=1:n
%      wx(:,i)=interp(y0(:,i),fact);
%   end
%   else, wx=y0;
%end
wx=y0;