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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 %