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