+
−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')
+
−%