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