Home > m > sigproc > frequency_domain > ltpda_arma_freq.m

ltpda_arma_freq

PURPOSE ^

LTPDA_ARMA_FREQ estimates the ARMA parameters of the transfer function

SYNOPSIS ^

function varargout = ltpda_arma_freq(varargin)

DESCRIPTION ^

 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.html,v 1.1 2008/03/01 12:29:33 hewitson Exp $

 HISTORY:     25-02-2008 M Nofrarias
                 Creation

 
 TODO: - Add parameters errors
       - Pass from univariate to multivariate
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_arma_freq(varargin)
0002 % LTPDA_ARMA_FREQ estimates the ARMA parameters of the transfer function
0003 % relating an output to an input
0004 %
0005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0006 %
0007 % DESCRIPTION: LTPDA_ARMA_FREQ estimates the ARMA parameters of the
0008 % transfer function relating an output to an input time series. The
0009 % algorithm uses the IV method to find a set of initial parameters
0010 % and then performs an iterative search based on the Optimization
0011 % Toolbox.
0012 %
0013 % CALL:        b = ltpda_arma_freq(ax,ay,pl)
0014 %
0015 % INPUTS:      ax  - analysis object containig the input time series
0016 %              ay  - analysis object containig the output time series
0017 %              pl  - parameters list
0018 %
0019 % OUTPUTS:     b  - cdata type analysis object containing the ARMA
0020 %                   parameters and the residual of the fit
0021 %
0022 % PARAMETERS:
0023 %  MaxIter     -  Maximum number of iterations to be performed by the
0024 %                 iterative search (default 4e2)
0025 %  MaxFunEvals -  Maximum number of function evaluations (default 1e2)
0026 %  TolX        -  Tolerance on the estimated parameter  (default 1e-6)
0027 %  TolFun      -  Tolerance on the evaluated function   (default 1e-6)
0028 %  UpBound     -  Array of parameters upper bounds for the iterative search
0029 %                 (default is 1 for each parameter)
0030 %  LowBound    -  Array of parameters lower bounds for the iterative search
0031 %                 (default is 1 for each parameter)
0032 %  ARMANum     -  MA order of the ARMA filter (default 2)
0033 %  ARMADen     -  AR order of the ARMA filter (default 1)
0034 %
0035 % VERSION:     $Id: ltpda_arma_freq.html,v 1.1 2008/03/01 12:29:33 hewitson Exp $
0036 %
0037 % HISTORY:     25-02-2008 M Nofrarias
0038 %                 Creation
0039 %
0040 %
0041 % TODO: - Add parameters errors
0042 %       - Pass from univariate to multivariate
0043 %
0044 %
0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0046 
0047 
0048 
0049 %% Standard history variables
0050 
0051 ALGONAME = mfilename;
0052 VERSION  = '$Id: ltpda_arma_freq.html,v 1.1 2008/03/01 12:29:33 hewitson Exp $';
0053 CATEGORY = 'SIGPROC';
0054 
0055 
0056 %% Check if this is a call for parameters
0057 if nargin == 1 && ischar(varargin{1})
0058   in = char(varargin{1});
0059   if strcmpi(in, 'Params')
0060     varargout{1} = getDefaultPL();
0061     return
0062   elseif strcmpi(in, 'Version')
0063     varargout{1} = VERSION;
0064     return
0065   elseif strcmpi(in, 'Category')
0066     varargout{1} = CATEGORY;
0067     return
0068   end
0069 end
0070 
0071 %% Capture input variables names
0072 invars = {};
0073 as     = [];
0074 ps     = [];
0075 for j=1:nargin
0076    invars = [invars cellstr(inputname(j))];
0077    if isa(varargin{j}, 'ao')
0078      as = [as varargin{j}];
0079    end
0080    if isa(varargin{j}, 'plist')
0081      ps = [ps varargin{j}];
0082    end
0083 end
0084 
0085 
0086 %% Check plist
0087 if isempty(ps)
0088   pl = getDefaultPL();
0089 else
0090   pl = combine(ps, getDefaultPL);
0091 end
0092 
0093 
0094 %% Fourier transform input varialbles
0095 
0096 Fs=0.65;
0097 L=length(xin)
0098 NFFT = 2^nextpow2(L); % Next power of 2 from length of y
0099 NFFT=length(xin);
0100 
0101 x_ = fft(xin,NFFT)/L;
0102 y_ = fft(yout',NFFT)/L;
0103 
0104 Px = 2*abs(x_(1:NFFT/2)); % single-sided amplitude spectrum.
0105 Py = 2*abs(y_(1:NFFT/2)); % single-sided amplitude spectrum.
0106 
0107 f = Fs/2*linspace(0,1,NFFT/2)';
0108 
0109 
0110 
0111 
0112 %% Iterative search using the Optimization Toolbox
0113 
0114 
0115 % Options for the lsqcurvefit function
0116 opt = optimset(...
0117 'MaxIter',find(pl, 'MaxIter'),...
0118 'TolX',find(pl, 'TolXm'),...
0119 'TolFun',find(pl, 'TolFun'),...
0120 'MaxFunEvals',find(pl, 'MaxFunEvals'),...
0121 'Display','off');
0122 
0123 % Upper and Lower Bounds for the parameters
0124 ub=find(pl,'UpBound');
0125 lb=find(pl,'LowBound');
0126 if isempty(ub)
0127 ub=ones(p+q,1);
0128 pl = append(pl, param('UpBound', ub));
0129 disp('! Parameters Upper Bound set to 1')
0130 end
0131 if isempty(lb)
0132 lb=-ones(p+q,1);
0133 pl = append(pl, param('LowBound', lb));
0134 disp('! Parameters Lower Bound set to -1')
0135 end    
0136 
0137 % Call to the Optimization Toolbox function
0138 [par,res]=lsqcurvefit(@(par,xdata)filter([par(1:p)],[1 par(p+1:p+q)],xdata),...
0139        par0,x,y,lb,ub,opt);
0140    
0141    
0142 %% Build output ao
0143  
0144 % New output history
0145 h = history(ALGONAME, VERSION, pl,[as(1).hist as(2).hist]);
0146 h = set(h,'invars', invars);
0147 
0148 % Make output analysis object
0149 b = ao([par; res]);
0150  
0151 % Name, mfilename, description for this object
0152 b = setnh(b,...
0153    'name', ['ARMA param. Input: ' sprintf('%s ', char(invars{1})) ' Output: ' char(invars{2})],...
0154    'mfilename', ALGONAME, ...
0155    'description', find(pl,'description'),...
0156    'hist', h);
0157  
0158 varargout{1} = b;
0159  
0160  
0161 %% --- FUNCTIONS ---
0162 
0163 
0164 % Get default parameters
0165 function plo = getDefaultPL();
0166 disp('* creating default plist...');
0167 plo = plist();
0168 plo = append(plo, param('MaxIter',4e2));
0169 plo = append(plo, param('MaxFunEvals',1e2));
0170 plo = append(plo, param('TolX',1e-6));
0171 plo = append(plo, param('TolFun',1e-6));
0172 plo = append(plo, param('ARMANum',2));
0173 plo = append(plo, param('ARMADen',1));
0174 plo = append(plo, param('fs', 1));
0175 disp('* done.');
0176 
0177 
0178 
0179 % siz=length(Px);
0180 % % siz=200;
0181 %
0182 % xdataf=Px(1:siz)*1e5;
0183 % ydataf=Py(1:siz)*1e5;
0184 % % xdataf=Px.data.y;
0185 % % ydataf=Py.data.y;
0186 % global freq
0187 % % freq=Px.data.x;
0188 % freq =f(1:siz);
0189 
0190 
0191 tic;
0192 [x1f,resnorm1f]=...
0193     lsqcurvefit(@myfunf,p0,xdataf,ydataf,[-1 -1 -1],[1 1 1],options);
0194 elaps=toc
0195 
0196 % %
0197 % %%
0198 % % Ff = ((x1f(1)+x1f(2)*exp(-i*2*pi*f/0.65))./(1+x1f(3)*exp(-i*2*pi*f/0.65)));
0199 % % F = ((x1(1)+x1(2)*exp(-i*2*pi*f/0.65))./(1+x1(3)*exp(-i*2*pi*f/0.65)));
0200 % %
0201 % % figure(1)
0202 % % hold off
0203 % % loglog(Px.data.x,ydataf./xdataf)
0204 % % hold
0205 % % loglog(Px.data.x,abs(Ff),'k')
0206 % % loglog(Px.data.x,abs(F),'m')
0207 % % 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')
0208 % %
0209 % %
0210 %
0211 % %%
0212 %
0213 % Fs=18;
0214 % Lw=2;
0215 %
0216 % Ff = ((x1f(1)+x1f(2)*exp(-i*2*pi*f/0.65))./(1+x1f(3)*exp(-i*2*pi*f/0.65)));
0217 % F = ((x1(1)+x1(2)*exp(-i*2*pi*f/0.65))./(1+x1(3)*exp(-i*2*pi*f/0.65)));
0218 %
0219 % figure(1)
0220 % set(gca,'Fontsize',Fs)
0221 % hold off
0222 % loglog(f,Py./Px,'LineWidth',Lw)
0223 % hold
0224 % % % loglog(f,ydataf,'m')
0225 % loglog(f,abs(Ff),'k','LineWidth',Lw)
0226 % loglog(f,abs(F),'m','LineWidth',Lw)
0227 %
0228 % % loglog(f,(x1(1)+x1(2)*exp(-i*2*pi*f/1.32))./(1+x1(3)*exp(-i*2*pi*f/1.32)),'r')
0229 % 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')
0230 %
0231 % legend('Empirical','Freq.domain','Time domain','Analytic')
0232 %

Generated on Tue 26-Feb-2008 10:52:52 by m2html © 2003