ltpda: m-toolbox/test/diagnostics/ltpda_arma

author	Daniele Nicolodi <nicolodi@science.unitn.it>
date	Wed, 23 Nov 2011 19:22:13 +0100
parents
children

rev	line source
0 f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	1 function varargout = ltpda_arma_freq(varargin)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	2 % LTPDA_ARMA_FREQ estimates the ARMA parameters of the transfer function
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	3 % relating an output to an input
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	4 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	6 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	7 % DESCRIPTION: LTPDA_ARMA_FREQ estimates the ARMA parameters of the
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	8 % transfer function relating an output to an input time series. The
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	9 % algorithm uses the IV method to find a set of initial parameters
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	10 % and then performs an iterative search based on the Optimization
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	11 % Toolbox.
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	12 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	13 % CALL: b = ltpda_arma_freq(ax,ay,pl)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	14 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	15 % INPUTS: ax - analysis object containig the input time series
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	16 % ay - analysis object containig the output time series
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	17 % pl - parameters list
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	18 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	19 % OUTPUTS: b - cdata type analysis object containing the ARMA
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	20 % parameters and the residual of the fit
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	21 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	22 % PARAMETERS:
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	23 % MaxIter - Maximum number of iterations to be performed by the
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	24 % iterative search (default 4e2)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	25 % MaxFunEvals - Maximum number of function evaluations (default 1e2)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	26 % TolX - Tolerance on the estimated parameter (default 1e-6)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	27 % TolFun - Tolerance on the evaluated function (default 1e-6)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	28 % UpBound - Array of parameters upper bounds for the iterative search
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	29 % (default is 1 for each parameter)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	30 % LowBound - Array of parameters lower bounds for the iterative search
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	31 % (default is 1 for each parameter)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	32 % ARMANum - MA order of the ARMA filter (default 2)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	33 % ARMADen - AR order of the ARMA filter (default 1)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	34 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	35 % VERSION: $Id: ltpda_arma_freq.m,v 1.1 2008/03/04 12:55:05 miquel Exp $
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	36 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	37 % HISTORY: 25-02-2008 M Nofrarias
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	38 % Creation
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	39 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	40 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	41 % TODO: - Add parameters errors
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	42 % - Pass from univariate to multivariate
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	43 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	44 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	46
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	47
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	48
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	49 %% Standard history variables
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	50
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	51 ALGONAME = mfilename;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	52 VERSION = '$Id: ltpda_arma_freq.m,v 1.1 2008/03/04 12:55:05 miquel Exp $';
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	53 CATEGORY = 'SIGPROC';
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	54
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	55
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	56 %% Check if this is a call for parameters
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	57 if nargin == 1 && ischar(varargin{1})
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	58 in = char(varargin{1});
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	59 if strcmpi(in, 'Params')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	60 varargout{1} = getDefaultPL();
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	61 return
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	62 elseif strcmpi(in, 'Version')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	63 varargout{1} = VERSION;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	64 return
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	65 elseif strcmpi(in, 'Category')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	66 varargout{1} = CATEGORY;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	67 return
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	68 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	69 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	70
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	71 %% Capture input variables names
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	72 invars = {};
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	73 as = [];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	74 ps = [];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	75 for j=1:nargin
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	76 invars = [invars cellstr(inputname(j))];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	77 if isa(varargin{j}, 'ao')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	78 as = [as varargin{j}];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	79 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	80 if isa(varargin{j}, 'plist')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	81 ps = [ps varargin{j}];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	82 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	83 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	84
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	85
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	86 %% Check plist
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	87 if isempty(ps)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	88 pl = getDefaultPL();
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	89 else
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	90 pl = combine(ps, getDefaultPL);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	91 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	92
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	93
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	94 %% Fourier transform input varialbles
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	95
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	96 Fs=0.65;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	97 L=length(xin)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	98 NFFT = 2^nextpow2(L); % Next power of 2 from length of y
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	99 NFFT=length(xin);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	100
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	101 x_ = fft(xin,NFFT)/L;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	102 y_ = fft(yout',NFFT)/L;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	103
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	104 Px = 2*abs(x_(1:NFFT/2)); % single-sided amplitude spectrum.
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	105 Py = 2*abs(y_(1:NFFT/2)); % single-sided amplitude spectrum.
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	106
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	107 f = Fs/2*linspace(0,1,NFFT/2)';
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	108
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	109
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	110
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	111
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	112 %% Iterative search using the Optimization Toolbox
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	113
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	114
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	115 % Options for the lsqcurvefit function
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	116 opt = optimset(...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	117 'MaxIter',find(pl, 'MaxIter'),...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	118 'TolX',find(pl, 'TolXm'),...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	119 'TolFun',find(pl, 'TolFun'),...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	120 'MaxFunEvals',find(pl, 'MaxFunEvals'),...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	121 'Display','off');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	122
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	123 % Upper and Lower Bounds for the parameters
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	124 ub=find(pl,'UpBound');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	125 lb=find(pl,'LowBound');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	126 if isempty(ub)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	127 ub=ones(p+q,1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	128 pl = append(pl, param('UpBound', ub));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	129 disp('! Parameters Upper Bound set to 1')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	130 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	131 if isempty(lb)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	132 lb=-ones(p+q,1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	133 pl = append(pl, param('LowBound', lb));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	134 disp('! Parameters Lower Bound set to -1')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	135 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	136
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	137 % Call to the Optimization Toolbox function
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	138 [par,res]=lsqcurvefit(@(par,xdata)filter([par(1:p)],[1 par(p+1:p+q)],xdata),...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	139 par0,x,y,lb,ub,opt);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	140
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	141
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	142 %% Build output ao
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	143
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	144 % New output history
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	145 h = history(ALGONAME, VERSION, pl,[as(1).hist as(2).hist]);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	146 h = set(h,'invars', invars);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	147
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	148 % Make output analysis object
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	149 b = ao([par; res]);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	150
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	151 % Name, mfilename, description for this object
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	152 b = setnh(b,...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	153 'name', ['ARMA param. Input: ' sprintf('%s ', char(invars{1})) ' Output: ' char(invars{2})],...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	154 'mfilename', ALGONAME, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	155 'description', find(pl,'description'),...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	156 'hist', h);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	157
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	158 varargout{1} = b;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	159
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	160
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	161 %% --- FUNCTIONS ---
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	162
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	163
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	164 % Get default parameters
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	165 function plo = getDefaultPL();
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	166 disp('* creating default plist...');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	167 plo = plist();
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	168 plo = append(plo, param('MaxIter',4e2));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	169 plo = append(plo, param('MaxFunEvals',1e2));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	170 plo = append(plo, param('TolX',1e-6));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	171 plo = append(plo, param('TolFun',1e-6));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	172 plo = append(plo, param('ARMANum',2));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	173 plo = append(plo, param('ARMADen',1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	174 plo = append(plo, param('fs', 1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	175 disp('* done.');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	176
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	177
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	178
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	179 % siz=length(Px);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	180 % % siz=200;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	181 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	182 % xdataf=Px(1:siz)*1e5;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	183 % ydataf=Py(1:siz)*1e5;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	184 % % xdataf=Px.data.y;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	185 % % ydataf=Py.data.y;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	186 % global freq
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	187 % % freq=Px.data.x;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	188 % freq =f(1:siz);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	189
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	190
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	191 tic;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	192 [x1f,resnorm1f]=...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	193 lsqcurvefit(@myfunf,p0,xdataf,ydataf,[-1 -1 -1],[1 1 1],options);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	194 elaps=toc
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	195
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	196 % %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	197 % %%
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	198 % % Ff = ((x1f(1)+x1f(2)exp(-i2pif/0.65))./(1+x1f(3)exp(-i2pif/0.65)));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	199 % % F = ((x1(1)+x1(2)exp(-i2pif/0.65))./(1+x1(3)exp(-i2pif/0.65)));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	200 % %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	201 % % figure(1)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	202 % % hold off
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	203 % % loglog(Px.data.x,ydataf./xdataf)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	204 % % hold
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	205 % % loglog(Px.data.x,abs(Ff),'k')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	206 % % loglog(Px.data.x,abs(F),'m')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	207 % % loglog(f,abs((39.6e-3-39.5e-3exp(-i2pif/0.65))./(1-0.996exp(-i2pif/0.65))),'r')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	208 % %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	209 % %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	210 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	211 % %%
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	212 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	213 % Fs=18;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	214 % Lw=2;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	215 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	216 % Ff = ((x1f(1)+x1f(2)exp(-i2pif/0.65))./(1+x1f(3)exp(-i2pif/0.65)));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	217 % F = ((x1(1)+x1(2)exp(-i2pif/0.65))./(1+x1(3)exp(-i2pif/0.65)));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	218 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	219 % figure(1)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	220 % set(gca,'Fontsize',Fs)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	221 % hold off
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	222 % loglog(f,Py./Px,'LineWidth',Lw)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	223 % hold
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	224 % % % loglog(f,ydataf,'m')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	225 % loglog(f,abs(Ff),'k','LineWidth',Lw)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	226 % loglog(f,abs(F),'m','LineWidth',Lw)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	227 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	228 % % loglog(f,(x1(1)+x1(2)exp(-i2pif/1.32))./(1+x1(3)exp(-i2pif/1.32)),'r')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	229 % loglog(f,abs((39.6e-3-39.5e-3exp(-i2pif/0.65))./(1-0.996exp(-i2pif/0.65))),'r')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	230 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	231 % legend('Empirical','Freq.domain','Time domain','Analytic')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	232 %

0

f0afece42f48 Import.