0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 1 % AUTOCFIT performs a fitting loop to identify model order and parameters.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 2 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 4 % DESCRIPTION:
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 5 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 6 % Perform a fitting loop to automatically identify model order and
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 7 % parameters in s-domain. Model identification is performed by 'vcfit'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 8 % function. Output is a s-domain model expanded in partial fractions:
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 9 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 10 % r1 rN
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 11 % f(s) = ------- + ... + ------- + d
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 12 % s - p1 s - pN
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 13 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 14 % The function attempt to perform first the identification of a model
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 15 % with d = 0, then if the operation do not succeed, it try the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 16 % identification with d different from zero.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 17 % Identification loop stop when the stop condition is reached. Six
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 18 % stop criteria are available:
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 19 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 20 % Mean Squared Error
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 21 % Check if the normalized mean squared error is lower than the value
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 22 % specified in lrscond:
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 23 % mse < 10^(-1*lrsvarcond)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 24 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 25 % Mean Squared Error and variation
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 26 % Check if the normalized mean squared error is lower than the value specified in
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 27 % lrscond and that the relative variation of the mean squared error is lower
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 28 % than the value provided.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 29 % Checking algorithm is:
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 30 % mse < 10^(-1*lrsvarcond)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 31 % Dmse < 10^(-1*msevar)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 32 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 33 % Log Residuals difference
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 34 % Check if the minimum of the logarithmic difference between data and
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 35 % residuals is larger than a specified value. ie. if the conditioning
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 36 % value is 2, the function ensures that the difference between data and
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 37 % residuals is at lest 2 order of magnitude lower than data itsleves.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 38 % Checking algorithm is:
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 39 % lsr = log10(abs(y))-log10(abs(rdl));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 40 % min(lsr) > lrscond;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 41 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 42 % Log Residuals difference and Root Mean Squared Error
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 43 % Check if the log difference between data and residuals is in
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 44 % larger than the value indicated in lsrcond and that the variation of
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 45 % the root mean squared error is lower than 10^(-1*msevar).
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 46 % Checking algorithm is:
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 47 % lsr = log10(abs(y))-log10(abs(rdl));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 48 % (lsr > lrscond) && (mse < 10^(-1*lrsvarcond));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 49 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 50 % Residuals Spectral Flatness
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 51 % In case of a fit on noisy data, the residuals from a good fit are
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 52 % expected to be as much as possible similar to a white noise. This
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 53 % property can be used to test the accuracy of a fit procedure. In
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 54 % particular it can be tested that the spectral flatness coefficient of
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 55 % the residuals is larger than a certain qiantity sf such that 0<sf<1.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 56 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 57 % Residuals Spectral Flatness and root mean squared error
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 58 % Check that the spectral flatness coefficient of the residuals is
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 59 % larger than a certain qiantity sf such that 0<sf<1 and that the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 60 % variation of the root mean squared error is lower than
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 61 % 10^(-1*msevar).
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 62 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 63 % Once the loop iteration stops the parameters giving best Mean Square
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 64 % Error are output.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 65 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 66 % CALL:
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 67 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 68 % [res,poles,dterm,mresp,rdl,mse] = autocfit(y,f,params)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 69 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 70 % INPUT:
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 71 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 72 % - y are the data to be fitted. They represent the frequency response
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 73 % of a certain process.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 74 % - f is the frequency vector in Hz
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 75 % - params is a struct containing identification parameters
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 76 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 77 % params.spolesopt = 0 --> use external starting poles
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 78 % params.spolesopt = 1 --> use real starting poles
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 79 % params.spolesopt = 2 --> use logspaced complex starting poles.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 80 % Default option
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 81 % params.spolesopt = 3 --> use linspaced complex starting poles.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 82 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 83 % params.extpoles = [] --> a vector with the starting poles.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 84 % Providing a fixed set of starting poles fixes the function order so
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 85 % params.minorder and params.maxorder will be internally set to the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 86 % poles vector length.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 87 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 88 % params.fullauto = 0 --> Perform a fitting loop as far as the number
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 89 % of iteration reach Nmaxiter. The order of the fitting function will
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 90 % be that specified in params.minorder. If params.dterm is setted to
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 91 % 1 the function will fit only with direct term.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 92 % params.fullauto = 1 --> Parform a full automatic search for the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 93 % transfer function order. The fitting procedure will stop when the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 94 % stopping condition defined in params.ctp is satisfied. Default
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 95 % value.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 96 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 97 % params.Nmaxiter = # --> Number of maximum iteration per model order
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 98 % parformed. Default is 50.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 99 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 100 % params.minorder = # --> Minimum model trial order. Default is 2.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 101 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 102 % params.maxorder = # --> Maximum model trial order. Default is 25.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 103 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 104 % params.weightparam = 0 --> use external weights
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 105 % params.weightparam = 1 --> fit with equal weights (one) for each
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 106 % data point.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 107 % params.weightparam = 2 --> weight fit with the inverse of absolute
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 108 % value of data. Default value.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 109 % params.weightparam = 3 --> weight fit with the square root of the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 110 % inverse of absolute value of data.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 111 % params.weightparam = 4 --> weight fit with inverse of the square
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 112 % mean spread
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 113 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 114 % params.extweights = [] --> A vector of externally provided weights.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 115 % It has to be of the same size of input data.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 116 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 117 % params.plot = 0 --> no plot during fit iteration
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 118 % params.plot = 1 --> plot results at each fitting steps. default
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 119 % value.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 120 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 121 % params.ctp = 'lrs' --> check if the log difference between data and
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 122 % residuals is point by point larger than the value indicated in
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 123 % lsrcond. This mean that residuals are lsrcond order of magnitudes
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 124 % lower than data.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 125 % params.ctp = 'lrsmse' --> check if the log difference between data
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 126 % and residuals is larger than the value indicated in lsrcond and if
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 127 % the variation of root mean squared error is lower than
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 128 % 10^(-1*msevar).
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 129 % params.ctp = 'rft' --> check that the residuals spectral flatness
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 130 % coefficient is lerger than the value provided in lsrcond. In this
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 131 % case lsrcond must be such that 0<lsrcond<1.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 132 % params.ctp = 'rftmse' --> check that the residuals spectral flatness
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 133 % coefficient is lerger than the value provided in lsrcond and if
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 134 % the variation of root mean squared error is lower than
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 135 % 10^(-1*msevar). In this case lsrcond must be such that
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 136 % 0<lsrcond<1.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 137 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 138 % params.lrscond = # --> set conditioning value for point to point
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 139 % log residuals difference (params.ctp = 'lsr' and params.ctp =
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 140 % 'lrsmse') or set conditioning value for residuals spectral
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 141 % flateness (params.ctp = 'rft' and params.ctp = 'rftmse'). In the
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 142 % last case params.lrscond must be set to 0<lrscond<1.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 143 % Default is 2. See help for stopfit.m for further remarks.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 144 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 145 % params.msevar = # --> set conditioning value for root mean squared
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 146 % error variation. This allow to check that the variation of root
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 147 % mean squared error is lower than 10^(-1*msevar).Default is 7. See
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 148 % help for stopfit.m for further remarks.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 149 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 150 % params.stabfit = 0 --> Fit without forcing poles stability. Default
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 151 % value.
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 152 % params.stabfit = 1 --> Fit forcing poles stability
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 153 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 154 % params.dterm = 0 --> Try to fit without direct term
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 155 % params.dterm = 1 --> Try to fit with and without direct term
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 156 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 157 % params.spy = 0 --> Do not display the iteration progression
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 158 % params.spy = 1 --> Display the iteration progression
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 159 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 160 % OUTPUT:
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 161 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 162 % - res is the vector with model residues r
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 163 % - poles is the vector with model poles p
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 164 % - dterm is the model direct term d
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 165 % - mresp is the model frequency response calculated at the input
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 166 % frequencies
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 167 % - rdl are the residuals between data and model, at the input
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 168 % frequencies
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 169 % - mse magnitude squared error progression
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 170 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 172 % VERSION: $Id: autocfit.m,v 1.22 2010/01/27 17:56:11 luigi Exp $
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 173 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 174 % HISTORY: 08-10-2008 L Ferraioli
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 175 % Creation
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 176 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 177 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 178
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 179 function [res,poles,dterm,mresp,rdl,mse] = autocfit(y,f,params)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 180
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 181 utils.helper.msg(utils.const.msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 182
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 183 % Default input struct
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 184 defaultparams = struct('spolesopt',1, 'Nmaxiter',30, 'minorder',2,...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 185 'maxorder',25, 'weightparam',1, 'plot',1,...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 186 'ctp','chival','lrscond',2,'msevar',2,...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 187 'stabfit',0,'dterm',0,'spy',0,'fullauto',1,'extweights', [],...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 188 'extpoles', []);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 189
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 190 names = {'spolesopt','Nmaxiter','minorder',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 191 'maxorder','weightparam','plot',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 192 'ctp','lrscond','msevar',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 193 'stabfit','dterm','spy','fullauto','extweights','extpoles'};
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 194
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 195 % collecting input and default params
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 196 if ~isempty(params)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 197 for jj=1:length(names)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 198 if isfield(params, names(jj)) && ~isempty(params.(names{1,jj}))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 199 defaultparams.(names{1,jj}) = params.(names{1,jj});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 200 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 201 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 202 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 203
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 204 % collecting input params
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 205 spolesopt = defaultparams.spolesopt;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 206 Nmaxiter = defaultparams.Nmaxiter;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 207 minorder = defaultparams.minorder;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 208 maxorder = defaultparams.maxorder;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 209 weightparam = defaultparams.weightparam;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 210 check = defaultparams.plot;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 211 stabfit = defaultparams.stabfit;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 212 ctp = defaultparams.ctp;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 213 lrscond = defaultparams.lrscond;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 214 msevar = defaultparams.msevar;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 215 idt = defaultparams.dterm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 216 spy = defaultparams.spy;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 217 autosearch = defaultparams.fullauto;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 218 extweights = defaultparams.extweights;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 219 extpoles = defaultparams.extpoles;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 220
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 221 if check == 1
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 222 fitin.plot = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 223 fitin.ploth = figure; % opening new figure window
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 224 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 225 fitin.plot = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 226 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 227
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 228 if stabfit % fit with stable poles only
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 229 fitin.stable = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 230 else % fit without restrictions
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 231 fitin.stable = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 232 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 233
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 234 % Colum vector are preferred
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 235 [a,b] = size(y);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 236 if a < b % shifting to column
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 237 y = y.';
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 238 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 239 [Nx,Ny] = size(y);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 240
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 241 [a,b] = size(f);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 242 if a < b % shifting to column
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 243 f = f.';
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 244 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 245
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 246 % in case of externally provided poles
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 247 if ~isempty(extpoles)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 248 spolesopt = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 249 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 250 if spolesopt == 0 % in case of external poles
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 251 % Colum vector are preferred
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 252 [a,b] = size(extpoles);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 253 if a < b % shifting to column
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 254 extpoles = extpoles.';
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 255 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 256 [Npls,b] = size(extpoles);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 257 minorder = Npls;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 258 maxorder = Npls;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 259 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 260
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 261 if weightparam == 0 % in case of external weights
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 262 % Colum vector are preferred
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 263 [a,b] = size(extweights);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 264 if a < b % shifting to column
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 265 extweights = extweights.';
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 266 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 267 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 268
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 269 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 270 % Importing package
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 271 import utils.math.*
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 272
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 273 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 274 % Fitting
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 275
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 276 % decide to fit with or without direct term according to the input
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 277 % options
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 278 if autosearch
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 279 if idt % full auto identification
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 280 dterm_off = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 281 dterm_on = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 282 else % auto ident without dterm
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 283 dterm_off = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 284 dterm_on = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 285 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 286 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 287 if idt % fit only with dterm
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 288 dterm_off = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 289 dterm_on = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 290 else % fit without dterm
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 291 dterm_off = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 292 dterm_on = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 293 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 294 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 295
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 296 ext = zeros(Ny,1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 297
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 298 % starting banana mse
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 299 bmse = inf;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 300 cmse = inf;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 301
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 302 if dterm_off
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 303 utils.helper.msg(utils.const.msg.PROC1, ' Try fitting without direct term ')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 304 fitin.dterm = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 305
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 306 % Weighting coefficients
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 307 if weightparam == 0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 308 % using external weigths
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 309 utils.helper.msg(utils.const.msg.PROC1, ' Using external weights... ')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 310 weight = extweights;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 311 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 312 weight = utils.math.wfun(y,weightparam);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 313 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 314
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 315 % Do not perform the loop if autosearch is setted to false
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 316 if autosearch
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 317 order_vect = minorder:maxorder;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 318 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 319 order_vect = minorder:minorder;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 320 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 321
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 322 for N = order_vect
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 323
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 324 if spy
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 325 utils.helper.msg(utils.const.msg.PROC1, ['Actual_Order' num2str(N)])
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 326 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 327
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 328 % Starting poles
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 329 if spolesopt == 0 % in case of external poles
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 330 utils.helper.msg(utils.const.msg.PROC1, ' Using external poles... ')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 331 spoles = extpoles;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 332 else % internally calculated starting poles
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 333 pparams = struct('spolesopt',spolesopt, 'type','CONT', 'pamp', 0.01);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 334 spoles = utils.math.startpoles(N,f,pparams);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 335 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 336
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 337 % Fitting
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 338 M = 2*N;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 339 if M > Nmaxiter
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 340 M = Nmaxiter;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 341 elseif not(autosearch)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 342 M = Nmaxiter;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 343 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 344
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 345 clear mlr
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 346
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 347 for hh = 1:M
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 348 [res,spoles,dterm,mresp,rdl,mse] = utils.math.vcfit(y,f,spoles,weight,fitin); % Fitting
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 349
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 350 % decide to store the best result based on mse
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 351 %fprintf('iteration = %d, order = %d \n',hh,N)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 352 if norm(mse)<cmse
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 353 %fprintf('nice job \n')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 354 bres = res;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 355 bpoles = spoles;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 356 bdterm = dterm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 357 bmresp = mresp;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 358 brdl = rdl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 359 bmse = mse;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 360 cmse = norm(mse);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 361 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 362
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 363 if spy
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 364 utils.helper.msg(utils.const.msg.PROC1, ['Iter' num2str(hh)])
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 365 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 366
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 367 % ext = zeros(Ny,1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 368 if autosearch
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 369 for kk = 1:Ny
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 370 % Stop condition checking
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 371 mlr(hh,kk) = mse(:,kk);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 372 % decide between stop conditioning
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 373 if strcmpi(ctp,'lrs')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 374 yd = y(:,kk); % input data
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 375 elseif strcmpi(ctp,'lrsmse')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 376 yd = y(:,kk); % input data
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 377 elseif strcmpi(ctp,'rft')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 378 yd = mresp(:,kk); % model response
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 379 elseif strcmpi(ctp,'rftmse')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 380 yd = mresp(:,kk); % model response
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 381 elseif strcmpi(ctp,'chival')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 382 yd = y(:,kk); % model response
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 383 elseif strcmpi(ctp,'chivar')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 384 yd = y(:,kk); % model response
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 385 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 386 error('!!! Unable to identify appropiate stop condition. See function help for admitted values');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 387 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 388 [next,msg] = utils.math.stopfit(yd,rdl(:,kk),mlr(:,kk),ctp,lrscond,msevar);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 389 ext(kk,1) = next;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 390 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 391 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 392 for kk = 1:Ny
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 393 % storing mse progression
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 394 mlr(hh,kk) = mse(:,kk);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 395 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 396 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 397
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 398 if all(ext)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 399 utils.helper.msg(utils.const.msg.PROC1, msg)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 400 break
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 401 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 402
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 403 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 404 if all(ext)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 405 break
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 406 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 407
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 408 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 409 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 410
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 411 if dterm_on
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 412 if ~all(ext) % fit with direct term only if the fit without does not give acceptable results (in full auto mode)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 413 utils.helper.msg(utils.const.msg.PROC1, ' Try fitting with direct term ')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 414 fitin.dterm = 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 415
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 416 if autosearch
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 417 order_vect = minorder:maxorder;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 418 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 419 order_vect = minorder:minorder;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 420 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 421
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 422 for N = order_vect
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 423
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 424 if spy
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 425 utils.helper.msg(utils.const.msg.PROC1, ['Actual_Order' num2str(N)])
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 426 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 427
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 428 % Starting poles
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 429 if spolesopt == 0 % in case of external poles
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 430 utils.helper.msg(utils.const.msg.PROC1, ' Using external poles... ')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 431 spoles = extpoles;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 432 else % internally calculated starting poles
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 433 pparams = struct('spolesopt',spolesopt, 'type','CONT', 'pamp', 0.01);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 434 spoles = utils.math.startpoles(N,f,pparams);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 435 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 436
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 437 % Fitting
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 438 M = 2*N;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 439 if M > Nmaxiter
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 440 M = Nmaxiter;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 441 elseif not(autosearch)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 442 M = Nmaxiter;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 443 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 444
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 445 clear mlr
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 446
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 447 for hh = 1:M
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 448 [res,spoles,dterm,mresp,rdl,mse] = utils.math.vcfit(y,f,spoles,weight,fitin); % Fitting
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 449
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 450 % decide to store the best result based on mse
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 451 if norm(mse)<cmse
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 452 bres = res;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 453 bpoles = spoles;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 454 bdterm = dterm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 455 bmresp = mresp;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 456 brdl = rdl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 457 bmse = mse;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 458 cmse = norm(mse);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 459 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 460
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 461 if spy
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 462 utils.helper.msg(utils.const.msg.PROC1, ['Iter' num2str(hh)])
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 463 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 464
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 465 ext = zeros(Ny,1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 466 if autosearch
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 467 for kk = 1:Ny
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 468 % Stop condition checking
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 469 mlr(hh,kk) = mse(:,kk);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 470 % decide between stop conditioning
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 471 if strcmpi(ctp,'lrs')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 472 yd = y(:,kk); % input data
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 473 elseif strcmpi(ctp,'lrsmse')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 474 yd = y(:,kk); % input data
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 475 elseif strcmpi(ctp,'rft')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 476 yd = mresp(:,kk); % model response
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 477 elseif strcmpi(ctp,'rftmse')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 478 yd = mresp(:,kk); % model response
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 479 elseif strcmpi(ctp,'chival')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 480 yd = y(:,kk); % model response
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 481 elseif strcmpi(ctp,'chivar')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 482 yd = y(:,kk); % model response
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 483 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 484 error('!!! Unable to identify appropiate stop condition. See function help for admitted values');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 485 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 486 [next,msg] = utils.math.stopfit(yd,rdl(:,kk),mlr(:,kk),ctp,lrscond,msevar);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 487 ext(kk,1) = next;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 488 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 489 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 490 for kk = 1:Ny
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 491 % storing mse progression
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 492 mlr(hh,kk) = mse(:,kk);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 493 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 494 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 495
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 496 if all(ext)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 497 utils.helper.msg(utils.const.msg.PROC1, msg)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 498 break
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 499 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 500
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 501 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 502 if all(ext)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 503 break
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 504 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 505
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 506 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 507
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 508 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 509 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 510
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 511 poles = bpoles;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 512 clear mse
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 513 mse = mlr(:,:);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 514
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 515 res = bres;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 516 dterm = bdterm;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 517 mresp = bmresp;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 518 rdl = brdl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 519 mse = bmse;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 520
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 521 if all(ext) == 0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 522 utils.helper.msg(utils.const.msg.PROC1, ' Fitting iteration completed without reaching the prescribed accuracy. Try changing Nmaxiter or maxorder or accuracy requirements ')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 523 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 524 utils.helper.msg(utils.const.msg.PROC1, ' Fitting iteration completed successfully ')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 525 end