--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/diagnostics/ltpda_arma_freq.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,232 @@
+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')
+%