function varargout = ltpda_arma_freq(varargin)
% LTPDA_ARMA_FREQ estimates the ARMA parameters of the transfer function
% relating an output to an input
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: LTPDA_ARMA_FREQ estimates the ARMA parameters of the
% transfer function relating an output to an input time series. The 
% algorithm uses the IV method to find a set of initial parameters 
% and then performs an iterative search based on the Optimization
% Toolbox.
%
% CALL:        b = ltpda_arma_freq(ax,ay,pl)
%
% INPUTS:      ax  - analysis object containig the input time series
%              ay  - analysis object containig the output time series
%              pl  - parameters list
%
% OUTPUTS:     b  - cdata type analysis object containing the ARMA
%                   parameters and the residual of the fit
%
% PARAMETERS:  
%  MaxIter     -  Maximum number of iterations to be performed by the  
%                 iterative search (default 4e2)
%  MaxFunEvals -  Maximum number of function evaluations (default 1e2)
%  TolX        -  Tolerance on the estimated parameter  (default 1e-6)
%  TolFun      -  Tolerance on the evaluated function   (default 1e-6)
%  UpBound     -  Array of parameters upper bounds for the iterative search
%                 (default is 1 for each parameter)
%  LowBound    -  Array of parameters lower bounds for the iterative search
%                 (default is 1 for each parameter)
%  ARMANum     -  MA order of the ARMA filter (default 2)
%  ARMADen     -  AR order of the ARMA filter (default 1)
%
% VERSION:     $Id: ltpda_arma_freq.m,v 1.1 2008/03/04 12:55:05 miquel Exp $
%
% HISTORY:     25-02-2008 M Nofrarias
%                 Creation
%
% 
% TODO: - Add parameters errors
%       - Pass from univariate to multivariate
% 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% Standard history variables 

ALGONAME = mfilename;
VERSION  = '$Id: ltpda_arma_freq.m,v 1.1 2008/03/04 12:55:05 miquel Exp $';
CATEGORY = 'SIGPROC';


%% Check if this is a call for parameters
if nargin == 1 && ischar(varargin{1})
  in = char(varargin{1});
  if strcmpi(in, 'Params')
    varargout{1} = getDefaultPL();
    return
  elseif strcmpi(in, 'Version')
    varargout{1} = VERSION;
    return
  elseif strcmpi(in, 'Category')
    varargout{1} = CATEGORY;
    return
  end
end

%% Capture input variables names
invars = {};
as     = [];
ps     = [];
for j=1:nargin
   invars = [invars cellstr(inputname(j))];
   if isa(varargin{j}, 'ao')
     as = [as varargin{j}];
   end
   if isa(varargin{j}, 'plist')
     ps = [ps varargin{j}];
   end
end


%% Check plist
if isempty(ps)
  pl = getDefaultPL();
else
  pl = combine(ps, getDefaultPL);
end


%% Fourier transform input varialbles

Fs=0.65;
L=length(xin)
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
NFFT=length(xin);

x_ = fft(xin,NFFT)/L;
y_ = fft(yout',NFFT)/L;

Px = 2*abs(x_(1:NFFT/2)); % single-sided amplitude spectrum.
Py = 2*abs(y_(1:NFFT/2)); % single-sided amplitude spectrum.

f = Fs/2*linspace(0,1,NFFT/2)';




%% Iterative search using the Optimization Toolbox


% Options for the lsqcurvefit function
opt = optimset(...
'MaxIter',find(pl, 'MaxIter'),...
'TolX',find(pl, 'TolXm'),...
'TolFun',find(pl, 'TolFun'),...
'MaxFunEvals',find(pl, 'MaxFunEvals'),...
'Display','off');

% Upper and Lower Bounds for the parameters 
ub=find(pl,'UpBound');
lb=find(pl,'LowBound');
if isempty(ub)
ub=ones(p+q,1);
pl = append(pl, param('UpBound', ub));
disp('! Parameters Upper Bound set to 1')
end
if isempty(lb)
lb=-ones(p+q,1);
pl = append(pl, param('LowBound', lb));
disp('! Parameters Lower Bound set to -1')
end    

% Call to the Optimization Toolbox function
[par,res]=lsqcurvefit(@(par,xdata)filter([par(1:p)],[1 par(p+1:p+q)],xdata),...
       par0,x,y,lb,ub,opt);
   
   
%% Build output ao
 
% New output history
h = history(ALGONAME, VERSION, pl,[as(1).hist as(2).hist]);
h = set(h,'invars', invars);

% Make output analysis object
b = ao([par; res]);
 
% Name, mfilename, description for this object
b = setnh(b,...
   'name', ['ARMA param. Input: ' sprintf('%s ', char(invars{1})) ' Output: ' char(invars{2})],...
   'mfilename', ALGONAME, ...
   'description', find(pl,'description'),...
   'hist', h);
 
varargout{1} = b;
 
 
%% --- FUNCTIONS ---


% Get default parameters
function plo = getDefaultPL();
disp('* creating default plist...');
plo = plist();
plo = append(plo, param('MaxIter',4e2));
plo = append(plo, param('MaxFunEvals',1e2));
plo = append(plo, param('TolX',1e-6));
plo = append(plo, param('TolFun',1e-6));
plo = append(plo, param('ARMANum',2));
plo = append(plo, param('ARMADen',1));
plo = append(plo, param('fs', 1));
disp('* done.');



% siz=length(Px);
% % siz=200;
% 
% xdataf=Px(1:siz)*1e5;
% ydataf=Py(1:siz)*1e5;
% % xdataf=Px.data.y;
% % ydataf=Py.data.y;
% global freq
% % freq=Px.data.x;
% freq =f(1:siz);


tic;
[x1f,resnorm1f]=...
    lsqcurvefit(@myfunf,p0,xdataf,ydataf,[-1 -1 -1],[1 1 1],options);
elaps=toc

% % 
% %%
% % Ff = ((x1f(1)+x1f(2)*exp(-i*2*pi*f/0.65))./(1+x1f(3)*exp(-i*2*pi*f/0.65)));
% % F = ((x1(1)+x1(2)*exp(-i*2*pi*f/0.65))./(1+x1(3)*exp(-i*2*pi*f/0.65)));
% % 
% % figure(1)
% % hold off
% % loglog(Px.data.x,ydataf./xdataf)
% % hold
% % loglog(Px.data.x,abs(Ff),'k')
% % loglog(Px.data.x,abs(F),'m')
% % loglog(f,abs((39.6e-3-39.5e-3*exp(-i*2*pi*f/0.65))./(1-0.996*exp(-i*2*pi*f/0.65))),'r')
% % 
% % 
% 
% %%
% 
% Fs=18;
% Lw=2;
% 
% Ff = ((x1f(1)+x1f(2)*exp(-i*2*pi*f/0.65))./(1+x1f(3)*exp(-i*2*pi*f/0.65)));
% F = ((x1(1)+x1(2)*exp(-i*2*pi*f/0.65))./(1+x1(3)*exp(-i*2*pi*f/0.65)));
% 
% figure(1)
% set(gca,'Fontsize',Fs)
% hold off
% loglog(f,Py./Px,'LineWidth',Lw)
% hold
% % % loglog(f,ydataf,'m')
% loglog(f,abs(Ff),'k','LineWidth',Lw)
% loglog(f,abs(F),'m','LineWidth',Lw)
% 
% % loglog(f,(x1(1)+x1(2)*exp(-i*2*pi*f/1.32))./(1+x1(3)*exp(-i*2*pi*f/1.32)),'r')
% loglog(f,abs((39.6e-3-39.5e-3*exp(-i*2*pi*f/0.65))./(1-0.996*exp(-i*2*pi*f/0.65))),'r')
% 
% legend('Empirical','Freq.domain','Time domain','Analytic')
%