function [ret1,ret2]=harmonic(observs,per,dt,tstart,outtype)
%HARMONIC harmonic analysis of time series
% HARMONIC computes a least squares solution fitting the input time 
% series with waves of input periods. 
%
% Inputs:  observs - matrix of observations (ASSUMED that each series in 
%                    a column
%          per     - vector of periods  
%          dt      - observation time step (sampling interval);  the time
%                    units must be consistent with the specification of the
%                    periods specified in per above
%          tstart  - starting time of the observation vector (optional, def=0)
%          outtype - string indication type of output (optional, def='ap')
%                    'ca'   ==>> complex amplitudes (An+iBn)
%                    'ap'   ==>> real amplitude and phase (radians)
%                    'wave' ==>> return a reconstruction of the input 
%                                signal composed of the amplitudes and 
%                                phases; each component is returned in
%                                a separate column.  This option 
%                                is valid only for 1 vector of observations.
%
% Output:  The output consists of the harmonic compoments for 
%          the periods specified, plus the DC component (average,
%          residual) as the FIRST value returned. The output can 
%          be returned in the following forms:
%
%          1) [a,p]=harmonic(...)  returns the amp and phase separately in a,p.
%          2) ap=harmonic(...)  returns the amp and phase as  ap=[a p].
%          3) [re,im]=harmonic(...) returns the real and imag parts in re,im.
%          4) ri=harmonic(...) returns the real and imag parts as ri=[An+iBn].
%          5) wave=harmonic(...) (only if number of series==1)     
%
%          In #2,4, the amp/phase (or real/imag) coeffs are returned interlaced,
%          as in :           Ser #1  Ser #2  ... 
%                  0-freq    A0 B0   A0 B0   ...
%                  1st-freq  A1 B1   A1 B1   ...
%                    ...      ...     ...
% 
%          One or two output arguments must be specified.     
%
%          The function PERIODS returns a structure of tidal 
%          constituent names and periods that can be passed to 
%          HARMONIC for harmonic analysis of tidal signals.
%
% Call as: wave=harmonic(observs,per,dt,tstart,'wave');
%             The "wave" variable will contain one column per frequency,
%             plus one column at the end for the residual.  If the 'per'
%             vector is length 2, 'wave' will be 3 columns wide
%             and the total wave can be reconstructed as: 
%                 totwave=wave(:,1)+wave(:,2)+wave(:,3);
%      OR: [a,b]=harmonic(observs,per,dt,tstart,'ca');
%      OR: [amp,pha]=harmonic(observs,per,dt,tstart,'ap');
%             The phase will be in rad/time_unit.
%
% Written by: Brian O. Blanton

%             Fall 1996
%             Spring 1999, mod
%

if nargin==0
   disp('Call as: wave=harmonic(observs,per,dt,tstart);')
   disp('     or  [a,b]=harmonic(observs,per,dt,tstart,''ca'');')
   disp('     or  ab=harmonic(observs,per,dt,tstart,''ca'');')
   disp('     or  [amp,pha]=harmonic(observs,per,dt,tstart,''ap'');')
   disp('     or  ap=harmonic(observs,per,dt,tstart,''ap'');')
   return
end

if nargout==0
   error('Harmonic must be called with 1|2 output arguments.')
end

if nargin<3|nargin>5
   disp('Wrong number of input arguments to HARMONIC')
   return
elseif nargin==3
   outtype='ap';
   tstart=0;
elseif nargin==4
   outtype='ap';
else
   if ~isstr(outtype)
      disp('Fifth argument (outtype) must be a string')
      return;
   end 
   if ~strcmp(outtype,'ca')&~strcmp(outtype,'ap')&~strcmp(outtype,'wave')
      disp('Output type unknown')
      return;
   end
   [M,nser]=size(observs);
   if nser>1 & strcmp(outtype,'wave')
      error('Wave output cannot be specified with more than one timeseries.')
   end
end

% check size of inputs
[M,nser]=size(observs);   % get the number of observations and time series
if M < nser   % report warning that number of time series exceeds the number of 
           % observations;  user should probably transpose.
   disp('WARNING: Number of time series vectors exceeds number of observations.')
   disp('Maybe you should transpose the observation matrix.')
end

[m,n]=size(per);
if m~=1 & n~=1
   disp('Period vector must be 1-D')
   return;
end
per=per(:);


[m,n]=size(dt);
if m*n~=1
   disp('Dt must be scalar')
   return;
 end
 
[m,n]=size(tstart);
if m*n~=1
   disp('tstart must be scalar')
   return;
end

% Convert period vector to freq in radians/time_unit
sig=2*pi./per;

N=length(per);
t=dt*(0:M-1);

m=1:M;
t=24*tstart+(m-1)*dt;
CS=ones(M,2*N+1);
CS(:,2:2:2*N+1)=cos(sig*t)';
CS(:,3:2:2*N+1)=sin(sig*t)';

% Allocate space for the solution vector, one per time series
% and amp, phase, real/imag coeffs
AB=ones(2*N+1,nser);
An=ones(N+1,nser);
Bn=ones(N+1,nser);
A=ones(N+1,nser);
P=ones(N+1,nser);

% Loop over each series separately
for i=1:nser
   AB(:,i)=CS\observs(:,i);    % "Least Squares"
   % Extract An and Bn coeffs from the least squares solution.
   An(1,i)=AB(1,i);  % This is the residual
   Bn(1,i)=0;
   for j=2:N+1
      An(j,i)=AB(2*j-2,i);
      Bn(j,i)=AB(2*j-1,i);
   end
end
% Compute Amplitude and Phase
A=sqrt(An.*An+Bn.*Bn);
P=atan2(Bn,An);

% Resynthesize wave; only if nser==1 (fix this, do for all series)
% 
if nser==1
   WT=([0;sig]*t)';
   PHASE=((P*ones(size(t)))');
   AMP=((A*ones(size(t)))');
   wave=AMP.*cos(WT-PHASE);
end


switch outtype
case 'ap'
   if nargout==2
      ret1=A;
      ret2=P;
   else
      ret1=ones(size([A A]));
      [nrow,ncol]=size(ret1);
      ret1(:,1:2:ncol)=A;
      ret1(:,2:2:ncol)=P;        
      ret2=[];
   end
case 'ca'
   if nargout==2
      ret1=a;
      ret2=b;
   else
      ret1=An+sqrt(-1)*Bn;
      ret2=[];
   end
case 'wave'
   ret1=wave;
   ret2=[];
end


%
%        Brian O. Blanton
%        Department of Marine Sciences
%        15-1A Venable Hall
%        CB# 3300
%        Uni. of North Carolina
%        Chapel Hill, NC
%                 27599-3300
%
%        919-962-4466
%        blanton@marine.unc.edu
%