From e5e316e1958e27696d7670e2492992d34ff38b68 Mon Sep 17 00:00:00 2001 From: ttt Date: Mon, 9 Jul 2018 16:54:44 +0530 Subject: added scilabs files --- Get.sci | 113 ++++++++++++++++++++++++++++++++ ar.sci | 102 +++++++++++++++++++++++++++++ armaX.sci | 51 +++++++++++++++ arx.sci | 110 +++++++++++++++++++++++++++++++ bj.sci | 122 +++++++++++++++++++++++++++++++++++ calModelPara.sci | 16 +++++ compare.sci | 145 +++++++++++++++++++++++++++++++++++++++++ compareBJ.sci | 56 ++++++++++++++++ dataSlice.sci | 67 +++++++++++++++++++ detrend.sci | 100 +++++++++++++++++++++++++++++ estpoly.sci | 136 +++++++++++++++++++++++++++++++++++++++ etfe.sci | 29 +++++++++ etfeeTest.sce | 13 ++++ fitch.sci | 42 ++++++++++++ frd.sci | 84 ++++++++++++++++++++++++ frdplot.sci | 11 ++++ iddata.sci | 68 ++++++++++++++++++++ iddataplot.sci | 69 ++++++++++++++++++++ identTime.sci | 32 ++++++++++ idpoly.sci | 142 ++++++++++++++++++++++++++++++++++++++++ impulseest.sci | 78 ++++++++++++++++++++++ iv.sci | 98 ++++++++++++++++++++++++++++ iv4.sci | 75 ++++++++++++++++++++++ misData.mat | Bin 0 -> 113297 bytes misdata.sci | 41 ++++++++++++ nInputData.sci | 15 +++++ oe.sci | 118 ++++++++++++++++++++++++++++++++++ pe.sci | 112 ++++++++++++++++++++++++++++++++ predict.sci | 192 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ rARX.mat | Bin 0 -> 20743 bytes rData.mat | Bin 0 -> 2319 bytes rarx.sci | 63 ++++++++++++++++++ sim.sci | 65 +++++++++++++++++++ spa.sci | 73 +++++++++++++++++++++ stem.sci | 19 ++++++ stepest.sci | 71 ++++++++++++++++++++ 36 files changed, 2528 insertions(+) create mode 100755 Get.sci create mode 100644 ar.sci create mode 100644 armaX.sci create mode 100644 arx.sci create mode 100644 bj.sci create mode 100644 calModelPara.sci create mode 100644 compare.sci create mode 100755 compareBJ.sci create mode 100644 dataSlice.sci create mode 100644 detrend.sci create mode 100644 estpoly.sci create mode 100644 etfe.sci create mode 100644 etfeeTest.sce create mode 100644 fitch.sci create mode 100644 frd.sci create mode 100644 frdplot.sci create mode 100755 iddata.sci create mode 100644 iddataplot.sci create mode 100644 identTime.sci create mode 100755 idpoly.sci create mode 100644 impulseest.sci create mode 100644 iv.sci create mode 100644 iv4.sci create mode 100644 misData.mat create mode 100644 misdata.sci create mode 100644 nInputData.sci create mode 100644 oe.sci create mode 100644 pe.sci create mode 100644 predict.sci create mode 100644 rARX.mat create mode 100644 rData.mat create mode 100644 rarx.sci create mode 100644 sim.sci create mode 100644 spa.sci create mode 100644 stem.sci create mode 100644 stepest.sci diff --git a/Get.sci b/Get.sci new file mode 100755 index 0000000..f14851a --- /dev/null +++ b/Get.sci @@ -0,0 +1,113 @@ +function []=Get(varargin) + //disp('Get') + tempData = varargin($) + f = fieldnames(tempData) +// disp(size(f,'*')) + if ~size(f,'*') then + if size(tempData,'*') == 0 then + printString = 'empty' + elseif isnan(tempData) then + printString = 'Nan' + elseif isinf(tempData) then + printString = 'inf' + else + printString = typeof(tempData) + end + mprintf('The type of ""Input"" is : %s \n',printString) + error(msprintf(gettext("%s: There is no objects in the ""Input"".\n"),"Get")) + else + maxLength = [] + for ii = 1:size(f,'*') + maxLength = [maxLength length(f(ii))] + end + maxLength = max(maxLength) + for ii = 1:size(f,'*') + blanckSpace = ' ' + for jj = 1:maxLength-length(f(ii)) + blanckSpace = blanckSpace + ' ' + end + mprintf('\t%s%s : ',blanckSpace,f(ii)) + objectData = tempData(f(ii)) + if typeof(objectData) == 'ce' then + //mprintf('cell') + sizeData = size(objectData) + prodData = prod(sizeData) + //pause + if prodData == 1 then + objectData = (cell2mat(objectData)) + objectData = getString(objectData) + tempString = '{' + objectData + '}' + mprintf('%s',tempString) + else + objectData = getString(sizeData,[]) + objectData = strsubst(objectData,']',' cell ]') + mprintf('%s',objectData) + end + elseif typeof(objectData) == 'string' then + if length(objectData) == 0 then + mprintf("{''''}") + else + mprintf('%s ',objectData) + end + + elseif typeof(objectData) == 'constant' then + //disp('in constant') + // + if ~size(objectData,'*') then + mprintf('%s%s','[',']') + //pause + //mprintf('%d',objectData) + elseif size(objectData,'*') == 1 then + if round(objectData)-objectData == 0 then + mprintf('%d',objectData) + else + mprintf('%.2f',objectData) + end + else + sizeData = size(objectData) + objectData = getString(sizeData,[]) + objectData = strsubst(objectData,']',' double ]') + mprintf('%s',objectData) + end + elseif typeof(objectData) == 'polynomial' then + if size(objectData,'*') == 1 then + polyData = pol2str(objectData) + mprintf('%s',polyData) + else + objectData = size(objectData) + objectData = getString(objectData,[]) + objectData = strsubst(objectData,']',' polynomial ]') + mprintf('%s',objectData) + end + elseif typeof(objectData) == 'hypermat' then + objectData = size(objectData) + objectData = getString(objectData,[]) + objectData = strsubst(objectData,']',' hypermat ]') + mprintf('%s',objectData) + end + mprintf('\n') + end + end +endfunction + + +// converting from matrix to string matrix +function Output = getString(varargin) + [lhs,rhs] = argn(0) + if rhs == 1 then + inBetween = ' ' + elseif rhs == 2 then + inBetween = 'x' + end + stringData = '[' + matData = varargin(1) + for ii = 1 : length(matData) + stringData = stringData + string(matData(ii)) //+ inBetween + if ii == length(matData) then + else + stringData = stringData + inBetween + end + end + stringData = stringData + ']' + Output = stringData +endfunction diff --git a/ar.sci b/ar.sci new file mode 100644 index 0000000..4239834 --- /dev/null +++ b/ar.sci @@ -0,0 +1,102 @@ +// Estimates Discrete time AR model +// A(q)y(t) = e(t) +// Current version uses random initial guess + +// Authors: Ashutosh,Harpreet,Inderpreet +// Updated(12-6-16) +function sys = ar(varargin) +// + [lhs , rhs] = argn(); + if ( rhs < 2 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "ar", rhs); + error(errmsg) + end + + z = varargin(1) + if typeof(z) == 'iddata' then + Ts = z.Ts;unit = z.TimeUnit + z = [z.OutputData z.InputData] + elseif typeof(z) == 'constant' then + Ts = 1;unit = 'seconds' + end + if ~iscolumn(z) then + errmsg = msprintf(gettext("%s: time series output data only"), "ar"); + error(errmsg); + end + + if (~isreal(z)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "ar"); + error(errmsg); + end + + n = varargin(2) + if (size(n,"*") ~=1 )then + errmsg = msprintf(gettext("%s: order should be nonnegative integer number "), "ar"); + error(errmsg); + end + + if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then + errmsg = msprintf(gettext("%s: values of order and delay matrix [na] should be nonnegative integer number "), "ar"); + error(errmsg); + end + + na = n; nb = 0; nk = 0; + // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + YDATA = z(:,1); + UDATA = zeros(size(z,1),1) + NDATA = size(UDATA,"*"); + function e = G(p,m) + e = YDATA - _objfun(UDATA,YDATA,p,na,nb,nk); + endfunction + tempSum = na+nb + p0 = linspace(0.1,0.9,tempSum)'; + [var,errl] = lsqrsolve(p0,G,size(UDATA,"*")); + err = (norm(errl)^2); + opt_err = err; + resid = G(var,[]); + a = 1-poly([var(nb+1:nb+na)]',"q","coeff"); + b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff"); + a = (poly([1,-coeff(a)],'q','coeff')) + t = idpoly(coeff(a),1,1,1,1,Ts) + + // estimating the other parameters + [temp1,temp2,temp3] = predict([YDATA UDATA],t) + [temp11,temp22,temp33] = pe([YDATA UDATA],t) + + estData = calModelPara(temp1,temp1,n(1)) + //pause + t.Report.Fit.MSE = estData.MSE + t.Report.Fit.FPE = estData.FPE + t.Report.Fit.FitPer = estData.FitPer + t.Report.Fit.AIC = estData.AIC + t.Report.Fit.AICc = estData.AICc + t.Report.Fit.nAIC = estData.nAIC + t.Report.Fit.BIC = estData.BIC + t.TimeUnit = unit + sys = t + //sys = idpoly(coeff(a),1,1,1,1,Ts) +// sys.TimeUnit = unit +endfunction + +function yhat = _objfun(UDATA,YDATA,x,na,nb,nk) + x=x(:) + q = poly(0,'q') + tempSum = nb+na + // making polynomials + b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff"); + a = 1 - poly([x(nb+1:nb+na)]',"q","coeff") + aSize = coeff(a);bSize = coeff(b) + maxDelay = max([length(aSize) length(bSize)]) + yhat = [YDATA(1:maxDelay)] + for k=maxDelay+1:size(UDATA,"*") + tempB = 0 + for ii = 1:size(bSize,'*') + tempB = tempB + bSize(ii)*UDATA(k-ii+1) + end + tempA = 0 + for ii = 1:size(aSize,"*") + tempA = tempA + aSize(ii)*YDATA(k-ii) + end + yhat = [yhat; [ tempA+tempB ]]; + end +endfunction diff --git a/armaX.sci b/armaX.sci new file mode 100644 index 0000000..93376e2 --- /dev/null +++ b/armaX.sci @@ -0,0 +1,51 @@ + +function sys = armaX(varargin) + [lhs , rhs] = argn(); + if ( rhs < 2 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "armaX", rhs); + error(errmsg) + end + + z = varargin(1) + if typeof(z) == 'iddata' then + Ts = z.Ts;unit = z.TimeUnit + z = [z.OutputData z.InputData] + elseif typeof(z) == 'constant' then + Ts = 1;unit = 'seconds' + end + if ((~size(z,2)==2) & (~size(z,1)==2)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "armaX"); + error(errmsg); + end + + if (~isreal(z)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "armaX"); + error(errmsg); + end + + n = varargin(2) + if (size(n,"*")<3| size(n,"*")>4) then + errmsg = msprintf(gettext("%s: The order and delay matrix [na nb nc nk] should be of size [3 or 4]"), "armaX"); + error(errmsg); + end + + if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then + errmsg = msprintf(gettext("%s: values of order and delay matrix [na nb nc nk] should be nonnegative integer number "), "armaX"); + error(errmsg); + end + + na = n(1); nb = n(2); nc = n(3); //nd = n(4);nf = n(5); + + if (size(n,"*") == 3) then + nk = 1 + else + nk = n(4); + end + + // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + YDATA = z(:,1); + UDATA = z(:,2); + + sys = estpoly(z,[na,nb,nc,0,0,nk]) + +endfunction diff --git a/arx.sci b/arx.sci new file mode 100644 index 0000000..0289cc8 --- /dev/null +++ b/arx.sci @@ -0,0 +1,110 @@ + +// Estimates Discrete time ARX model +// A(q)y(t) = B(q)u(t) + e(t) +// Current version uses random initial guess +// + +// Authors: Ashutosh,Harpreet,Inderpreet +// Updated(12-6-16) +function sys = arx(varargin) + [lhs , rhs] = argn(); + if ( rhs < 2 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "arx", rhs); + error(errmsg) + end + + z = varargin(1) + if typeof(z) == 'iddata' then + Ts = z.Ts;unit = z.TimeUnit + z = [z.OutputData z.InputData] + elseif typeof(z) == 'constant' then + Ts = 1;unit = 'seconds' + end + if ((~size(z,2)==2) & (~size(z,1)==2)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "arx"); + error(errmsg); + end + + if (~isreal(z)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "arx"); + error(errmsg); + end + + n = varargin(2) + if (size(n,"*")<2| size(n,"*")>3) then + errmsg = msprintf(gettext("%s: The order and delay matrix [na nb nk] should be of size [2 3]"), "arx"); + error(errmsg); + end + + if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then + errmsg = msprintf(gettext("%s: values of order and delay matrix [na nb nk] should be nonnegative integer number "), "arx"); + error(errmsg); + end + + na = n(1); nb = n(2); //nk = n(3); //nf = n(4); +// + if (size(n,"*") == 2) then + nk = 1 + else + nk = n(3); + end + + // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + YDATA = z(:,1); + UDATA = z(:,2); + NDATA = size(UDATA,"*"); + function e = G(p,m) + e = YDATA - _objfunarx(UDATA,YDATA,p,na,nb,nk); + endfunction + tempSum = na+nb + p0 = linspace(0.1,0.9,tempSum)'; + [var,errl] = lsqrsolve(p0,G,size(UDATA,"*")); + err = (norm(errl)^2); + opt_err = err; + resid = G(var,[]); + + a = 1-poly([var(nb+1:nb+na)]',"q","coeff"); + b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff"); + a = (poly([1,-coeff(a)],'q','coeff')) + t = idpoly(coeff(a),coeff(b),1,1,1,Ts) + + // estimating the other parameters + [temp1,temp2,temp3] = predict(z,t) + [temp11,temp22,temp33] = pe(z,t) + + estData = calModelPara(temp1,temp1,n(1)+n(2)) + //pause + t.Report.Fit.MSE = estData.MSE + t.Report.Fit.FPE = estData.FPE + t.Report.Fit.FitPer = estData.FitPer + t.Report.Fit.AIC = estData.AIC + t.Report.Fit.AICc = estData.AICc + t.Report.Fit.nAIC = estData.nAIC + t.Report.Fit.BIC = estData.BIC + t.TimeUnit = unit + sys = t + +endfunction + +function yhat = _objfunarx(UDATA,YDATA,x,na,nb,nk) + x=x(:) + q = poly(0,'q') + tempSum = nb+na + // making polynomials + b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff"); + a = 1 - poly([x(nb+1:nb+na)]',"q","coeff") + aSize = coeff(a);bSize = coeff(b) + maxDelay = max([length(aSize) length(bSize)]) + yhat = [YDATA(1:maxDelay)] + for k=maxDelay+1:size(UDATA,"*") + tempB = 0 + for ii = 1:size(bSize,'*') + tempB = tempB + bSize(ii)*UDATA(k-ii+1) + end + tempA = 0 + for ii = 1:size(aSize,"*") + tempA = tempA + aSize(ii)*YDATA(k-ii) + end + yhat = [yhat; [ tempA+tempB ]]; + end +endfunction diff --git a/bj.sci b/bj.sci new file mode 100644 index 0000000..f41ce65 --- /dev/null +++ b/bj.sci @@ -0,0 +1,122 @@ +// Estimates Discrete time BJ model +// y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t) +// Current version uses random initial guess +// Need to get appropriate guess from OE and noise models + +// Authors: Ashutosh,Harpreet,Inderpreet +// Updated(12-6-16) + +//function [theta_bj,opt_err,resid] = bj(varargin) +function sys = bj(varargin) + + [lhs , rhs] = argn(); + if ( rhs < 2 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "bj", rhs); + error(errmsg) + end + + z = varargin(1) + if typeof(z) == 'iddata' then + Ts = z.Ts;unit = z.TimeUnit + z = [z.OutputData z.InputData] + elseif typeof(z) == 'constant' then + Ts = 1;unit = 'seconds' + end + if ((~size(z,2)==2) & (~size(z,1)==2)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "bj"); + error(errmsg); + end + + if (~isreal(z)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "bj"); + error(errmsg); + end + + n = varargin(2) + if (size(n,"*")<4| size(n,"*")>5) then + errmsg = msprintf(gettext("%s: The order and delay matrix [nb nc nd nf nk] should be of size [4 5]"), "bj"); + error(errmsg); + end + + if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then + errmsg = msprintf(gettext("%s: values of order and delay matrix [nb nc nd nf nk] should be nonnegative integer number "), "bj"); + error(errmsg); + end + + nb = n(1); nc = n(2); nd = n(3); nf = n(4); + + if (size(n,"*") == 4) then + nk = 1 + else + nk = n(5); + end + + // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + YDATA = z(:,1); + UDATA = z(:,2); + NDATA = size(UDATA,"*"); + function e = G(p,m) + e = YDATA - _objfun(UDATA,p,nd,nc,nf,nb,nk); + endfunction + tempSum = nb+nc+nd+nf + p0 = linspace(0.5,0.9,tempSum)'; + [var,errl] = lsqrsolve(p0,G,size(UDATA,"*")); + + err = (norm(errl)^2); + opt_err = err; + resid = G(var,[]); + b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff"); + c = poly([1; var(nb+1:nb+nc)]',"q","coeff"); + d = poly([1; var(nb+nc+1:nb+nc+nd)]',"q","coeff"); + f = poly([1; var(nb+nd+nc+1:nd+nc+nf+nb)]',"q","coeff"); + t = idpoly(1,coeff(b),coeff(c),coeff(d),coeff(f),Ts) + + // estimating the other parameters + [temp1,temp2,temp3] = predict(z,t) + [temp11,temp22,temp33] = pe(z,t) + + estData = calModelPara(temp1,temp11,sum(n(1:4))) + //pause + t.Report.Fit.MSE = estData.MSE + t.Report.Fit.FPE = estData.FPE + t.Report.Fit.FitPer = estData.FitPer + t.Report.Fit.AIC = estData.AIC + t.Report.Fit.AICc = estData.AICc + t.Report.Fit.nAIC = estData.nAIC + t.Report.Fit.BIC = estData.BIC + t.TimeUnit = unit + sys = t +endfunction + +function yhat = _objfun(UDATA,x,nd,nc,nf,nb,nk) + x=x(:) + q = poly(0,'q') + tempSum = nb+nc+nd+nf + // making polynomials + b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff"); + c = poly([1; x(nb+1:nb+nc)]',"q","coeff"); + d = poly([1; x(nb+nc+1:nb+nc+nd)]',"q","coeff"); + f = poly([1; x(nb+nd+nc+1:nd+nc+nf+nb)]',"q","coeff"); + bd = coeff(b*d); cf = coeff(c*f); fc_d = coeff(f*(c-d)); + if size(bd,"*") == 1 then + bd = repmat(0,nb+nd+1,1) + end + maxDelay = max([length(bd) length(cf) length(fc_d)]) + yhat = [YDATA(1:maxDelay)] + for k=maxDelay+1:size(UDATA,"*") + bdadd = 0 + for i = 1:size(bd,"*") + bdadd = bdadd + bd(i)*UDATA(k-i+1) + end + + fc_dadd = 0 + for i = 1:size(fc_d,"*") + fc_dadd = fc_dadd + fc_d(i)*YDATA(k-i+1) + end + cfadd = 0 + for i = 2:size(cf,"*") + cfadd = cfadd + cf(i)*yhat(k-i+1) + end + yhat = [yhat; [ bdadd + fc_dadd - cfadd ]]; + end +endfunction diff --git a/calModelPara.sci b/calModelPara.sci new file mode 100644 index 0000000..e8fbc02 --- /dev/null +++ b/calModelPara.sci @@ -0,0 +1,16 @@ +function varargout = calModelPara(varargin) + y = varargin(1)+varargin(2) + N = size(y,'r') + ek = varargin(2) + np = varargin(3) + //pause + mse = sum(ek.^2)/N;//disp(mse) + fpe = mse * (1 + np/N)/(1 - np/N);//disp(fpe) + nrmse = 1 - sqrt(sum(ek^2))/sqrt(sum((y - mean(y))^2));//disp(nrmse) + AIC = N * log(mse) + 2 * np + N * size(y,'c') * (log(2 * %pi) + 1);//disp(AIC) + AICc = AIC * 2 * np * (np + 1)/(N - np - 1);//disp(AICc) + nAIC = log(mse) + 2 * np/N;//disp(nAIC) + BIC = N * log(mse) + N * size(y,'c') * (log(2 * %pi) + 1) + np * log(N);//disp(BIC) + //pause + varargout(1) = struct('MSE',mse,'FPE',fpe,'FitPer',nrmse*100,'AIC',AIC,'AICc',AICc,'nAIC',nAIC,'BIC',BIC) +endfunction diff --git a/compare.sci b/compare.sci new file mode 100644 index 0000000..4723130 --- /dev/null +++ b/compare.sci @@ -0,0 +1,145 @@ +function varargout = compare(varargin) + [lhs,rhs] = argn(0) +//------------------------------------------------------------------------------ +// checking the number of inputs + if rhs < 2 || rhs > 3 then + error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"compare")) + end +//------------------------------------------------------------------------------ + data = varargin(1) + model = varargin(2) + if rhs == 3 then + kStep = varargin(3) + elseif rhs == 2 then + kStep = %inf + end + +//------------------------------------------------------------------------------ + +// k step analysis + if typeof(kStep) <> 'constant' || isnan(kStep) then + error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"compare")) + end +// if given k step is infinity or [] +// if isinf(kStep) || ~size(kStep,'*') then +// kStep = 1 +// end +// checking the dimensions +// if size(kStep,'*') <> 1 || (ceil(kStep)-kStep) then +// error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"predict")) +// end +//------------------------------------------------------------------------------ +// checking the plant model + if typeof(model) ~= 'idpoly' then + error(msprintf(gettext("%s:Plant model must be ""idpoly"" type.\n"),"compare")) + end + modelSampleTime = model.Ts + modelTimeUnit = model.TimeUnit +//------------------------------------------------------------------------------ +//checking the data type + if typeof(data) <> 'iddata' && typeof(data) <> 'constant' then + error(msprintf(gettext("%s:Sample data must be ""iddata"" type or ""n x 2"" matrix type.\n"),"compare")) + end +// checking the plant data + if typeof(data) == 'iddata' then + if ~size(data.OutputData,'*') || ~size(data.InputData,'*') then + error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"compare")) + end + plantSampleTime = data.Ts + plantTimeUnit = data.TimeUnit + data = [data.OutputData data.InputData] + //disp('iddata') + elseif typeof(data) == 'constant' then + if size(data,'c') ~= 2 then + error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"compare")) + end + plantSampleTime = model.Ts + plantTimeUnit = model.TimeUnit + end +//------------------------------------------------------------------------------ +// comparing the sampling time + if modelSampleTime-plantSampleTime <> 0 then + error(msprintf(gettext("%s:The sample time of the model and plant data must be equal.\n"),"compare")) + end +// Comparing the time units + if ~strcmp(modelTimeUnit,plantTimeUnit) then + else + error(msprintf(gettext("%s:Time unit of the model and plant data must be equal.\n"),"compare")) + end +//------------------------------------------------------------------------------ +// ckecking the k step size. if it greater than number of sample size then the +// k step will become 1 +// if kStep >= size(data,'r') then +// kStep = 1 +// end +//------------------------------------------------------------------------------ + //storing the plant data + // B(z) C(z) + // y(n) = ---------- u(n) + ---------- e(n) + // A(z)*F(z) A(z)*D(z) + aPoly = poly(model.a,'q','coeff'); + bPoly = poly(model.b,'q','coeff'); + cPoly = poly(model.c,'q','coeff'); + dPoly = poly(model.d,'q','coeff'); + fPoly = poly(model.f,'q','coeff'); + Gq = bPoly/(aPoly*fPoly) + Hq = cPoly/(aPoly*dPoly) + //disp(kStep) + if kStep == %inf then + //disp('in inf') + outData = sim(data(:,2),model) + else + if kStep == 1 then + Wkq = Hq^-1 + elseif kStep > 1 then + adCoeff = coeff(aPoly*dPoly);adCoeff = adCoeff(length(adCoeff):-1:1); + adPoly = poly(adCoeff,'q','coeff') + cCoeff = model.c;cCoeff = cCoeff(length(cCoeff):-1:1); + cPoly = poly(cCoeff,'q','coeff') + hBar = clean((ldiv(cPoly,adPoly,kStep))',0.00001) + hBarPoly = poly(hBar,'q','coeff') + Wkq = hBarPoly*Hq^-1//*bPoly/(dPoly*fPoly) + end + // pause + WkqGq = Wkq * Gq + tempWkqGq = coeff(WkqGq.den) + if tempWkqGq(1) <> 1 then + WkqGq.num = WkqGq.num/tempWkqGq(1) + WkqGq.den = WkqGq.den/tempWkqGq(1) + end + Wkq1 = 1-Wkq + tempWkq1 = coeff(Wkq1.den) + if tempWkq1(1) == 1 then + Wkq1.num = Wkq1.num/tempWkq1(1) + Wkq1.den = Wkq1.den/tempWkq1(1) + end + z = poly(0,'z'); + WkqGqNum = horner(WkqGq.num,1/z);WkqGqDen = horner(WkqGq.den,1/z); + uPoly = WkqGqNum/WkqGqDen; + uPoly = syslin('d',uPoly);uData = flts(data(:,2)',uPoly); + + Wkq1Num = horner(Wkq1.num,1/z);Wkq1Den = horner(Wkq1.den,1/z); + yPoly = Wkq1Num/Wkq1Den; + yPoly = syslin('d',yPoly);yData = flts(data(:,1)',yPoly); + outData = (uData+yData)' + end + + tData = (1:size(data,'r'))'*plantSampleTime + fitData = fitValue(data(:,1),outData) + if lhs == 1 then + varargout(1) = [] + plot(tData,data(:,1),'m') + plot(tData,outData,'b') + legend('Plant Data','Model Data : '+string(fitData)+'%') + xtitle('Comparison of Time Response','Time ('+ plantTimeUnit+')','Amplitude') + xgrid() + elseif lhs == 2 then + varargout(1) = fitData + varargout(2) = outData + elseif lhs == 3 then + varargout(1) = outData + varargout(2) = tData + varargout(3) = fitData + end + +endfunction diff --git a/compareBJ.sci b/compareBJ.sci new file mode 100755 index 0000000..06d2fb8 --- /dev/null +++ b/compareBJ.sci @@ -0,0 +1,56 @@ +function varargout = compareBJ(varargin) + //varargin(1) -> idpoly data about oe + //varargin(2) -> [y u] matrix of "nx2" dimension + //disp('compareBj') + bjData = varargin(1) + //disp(typeof(bjData)) + plantData = varargin(2) + //disp(typeof(plantData)) + yData = plantData(:,1) + uData = plantData(:,2) + uData = [0;uData] + //storing the data in A,B,C,D,F matrix + // B(z) C(z) + // y(n)=---------- u(n)+ ----------- e(n) + // A(z)F(z) A(z)D(z) + polyA = poly(bjData.a,'x','coeff') + polyB = poly(bjData.b,'x','coeff') + polyC = poly(bjData.c,'x','coeff') + polyD = poly(bjData.d,'x','coeff') + polyF = poly(bjData.f,'x','coeff') + adf = polyA*polyD*polyF + bd = polyB*polyD + cf = polyC*polyF + delay = max(size(coeff(adf),'*'),size(coeff(bd),'*'),size(coeff(cf),'*')) + yHat = [0] + bdCoeff = coeff(bd) + adfCoeff = coeff(adf) + adfCoeff = -adfCoeff(2:length(adfCoeff)) + for ii = 1:length(uData) + uSum = 0;ySum = 0; + for jj = 1:length(bdCoeff) + if ii-jj <= 0 then + uSum = uSum + 0 + else + uSum = uSum + uData(ii-jj+1)*bdCoeff(jj) + end + end + for jj = 1:length(adfCoeff) + if ii-jj <= 0 then + ySum = ySum + 0 + else + ySum = ySum + yHat(ii-jj+1)*adfCoeff(jj) + end + end + yHat = [yHat; uSum+ySum] + end + tempStart = 1 + if size(yHat,'r')- size(yData,'r') > 0 then + tempStart = size(yHat,'r')-size(yData,'r')+1 + end + varargout(1) = yHat(tempStart:length(yHat)) +// plot(yHat(tempStart:length(yHat)),'m') +// plot(yData) +// xgrid() +// pause +endfunction diff --git a/dataSlice.sci b/dataSlice.sci new file mode 100644 index 0000000..6b3c160 --- /dev/null +++ b/dataSlice.sci @@ -0,0 +1,67 @@ +function sys = dataSlice(data,Start,End,Freq) + [lhs,rhs] = argn() + // storing the model data + modelData = data + // storing the statrting point + try + startData = Start + catch + startData = 1 + end + //storing the end point + try + endData = End + catch + endData = LastIndex(data) + end + //Storing the frequency + try + freqData = Freq + catch + freqData = 1 + end + // error message generate + if startData > endData then + error(msprintf(gettext("%s:Start index can not greater than End index.\n"),"dataSlice")) + end + if size(startData,'*') ~= 1 then + error(msprintf(gettext("%s:Start index must be non negative scalar integer number.\n"),"dataSlice")) + end + if size(endData,'*') ~= 1 then + error(msprintf(gettext("%s:End index must be non negative scalar integer number.\n"),"dataSlice")) + end + if ~freqData || size(freqData,'*') ~= 1 then + error(msprintf(gettext("%s:Frequency must be non negative scalar number.\n"),"dataSlice")) + end + //-------------------------------------------------------------------------- + if typeof(modelData) == 'constant' then + Ts = 1 + elseif typeof(modelData) == 'iddata' then + Ts = modelData.Ts + end + //-------------------------------------------------------------------------- + if freqData> Ts || modulo(Ts,freqData) then + warning(msprintf(gettext("%s:inconsistent frequency.\n"),"dataSlice")) + freqData = Ts + end + if typeof(modelData) == 'constant' then + temp = modelData(startData:Ts/freqData:endData,:) + elseif typeof(modelData) == 'iddata' then + tempY = modelData.OutputData;tempU = modelData.InputData + tempY = tempY(startData:Ts/freqData:endData,:);tempU = tempU(startData:Ts/freqData:endData,:) + temp = iddata(tempY,tempU,Ts/freqData) + temp.TimeUnit = modelData.TimeUnit + end + sys = temp +endfunction + +function varargout = LastIndex(modelData) + //finding the sample size + if typeof(modelData) == "constant" then + varargout(1) = length(modelData(:,1)) + + elseif typeof(modelData) == "iddata" then + temp = modelData.OutputData + varargout(1) = length(temp(:,1)) + end +endfunction diff --git a/detrend.sci b/detrend.sci new file mode 100644 index 0000000..993857e --- /dev/null +++ b/detrend.sci @@ -0,0 +1,100 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) Bruno Pincon +// +// Copyright (C) 2012 - 2016 - Scilab Enterprises +// +// This file is hereby licensed under the terms of the GNU GPL v2.0, +// pursuant to article 5.3.4 of the CeCILL v.2.1. +// This file was originally licensed under the terms of the CeCILL v2.1, +// and continues to be available under such terms. +// For more information, see the COPYING file which you should have received +// along with this program. +// +function [y] = detrend(x, flag, bp) + // + // this function removes the constant or linear or + // piecewise linear trend from a vector x. If x is + // matrix this function removes the trend of each + // column of x. + // + // flag = "constant" or "c" to removes the constant trend + // (simply the mean of the signal) + // flag = "linear" or "l" to remove the linear or piecewise + // linear trend. + // + // To define the piecewise linear trend the function needs the + // breakpoints and these must be given as the third argument bp. + // + // The "instants" of the signal x are in fact from 0 to m-1 + // (m = length(x) if x is a vector and m = size(x,1) in case + // x is a matrix). So bp must be reals in [0 m-1]. + // + + rhs = argn(2) + if rhs < 1 | rhs > 3 then + error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"detrend",1,3)); + elseif rhs == 1 + flag = "linear"; bp = [] + elseif rhs == 2 + bp = [] + end + + if type(x)~=1 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",1)); + end + if type(flag)~=10 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",2)); + end + if ~(type(bp)==1 & isreal(bp)) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",3)); + end + + [mx,nx] = size(x) + if mx==1 | nx==1 then + x_is_vector = %t; x = x(:); m = mx*nx; n = 1 + elseif mx*nx == 0 then + y = []; return + else + x_is_vector = %f; m = mx; n = nx + end + + + if flag == "constant" | flag == "c" then + y = x - ones(m,1)*mean(x,"r") + elseif flag == "linear" | flag == "l" + bp = unique([0 ; bp(:) ; m-1]) + // delete all the breakpoints outside [0,m-1] + while bp(1) < 0, bp(1)=[], end + while bp($) > m-1, bp($)=[], end + // breakpoints are 0-based so add one to + // compare them with signal vector indices (1-based) + bp = bp + 1; + nbp = length(bp); + // build the least square matrix with hat functions + // (as a basis for continuous piecewise linear functions) + A = zeros(m, nbp) + k1 = 1 + delta_bp = diff(bp) + for j = 2:nbp-1 + k2 = ceil(bp(j))-1 + ind = (k1:k2)' + A(ind,j-1) = (bp(j) - ind)/delta_bp(j-1) + A(ind,j) = (ind - bp(j-1))/delta_bp(j-1) + k1 = k2+1 + end + ind = (k1:m)' + A(ind,nbp-1) = (m - ind)/delta_bp(nbp-1) + A(ind,nbp) = (ind - bp(nbp-1))/delta_bp(nbp-1) + // solve the least square pb and retrieve the fitted + // piecewise linear func off the signal + y = x - A*(A\x) + else + error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n") ,.. + "detrend",2,"''constant'',''c'',''linear'',''l''")); + end + + if x_is_vector then + y = matrix(y,mx,nx) + end + +endfunction diff --git a/estpoly.sci b/estpoly.sci new file mode 100644 index 0000000..f1ea48d --- /dev/null +++ b/estpoly.sci @@ -0,0 +1,136 @@ +// Estimates Discrete time estpoly model +// y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t) +// Current version uses random initial guess +// Need to get appropriate guess from OE and noise models + +// Authors: Ashutosh,Harpreet,Inderpreet +// Updated(12-6-16) + +//function [theta_estpoly,opt_err,resid] = estpoly(varargin) +function sys = estpoly(varargin) + + [lhs , rhs] = argn(); + if ( rhs < 2 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "estpoly", rhs); + error(errmsg) + end + + z = varargin(1) + if typeof(z) == 'iddata' then + Ts = z.Ts;unit = z.TimeUnit + z = [z.OutputData z.InputData] + elseif typeof(z) == 'constant' then + Ts = 1;unit = 'seconds' + end + if ((~size(z,2)==2) & (~size(z,1)==2)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "estpoly"); + error(errmsg); + end + + if (~isreal(z)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "estpoly"); + error(errmsg); + end + + n = varargin(2) + if (size(n,"*")<5| size(n,"*")>6) then + errmsg = msprintf(gettext("%s: The order and delay matrix [na nb nc nd nf nk] should be of size [5 6]"), "estpoly"); + error(errmsg); + end + + if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then + errmsg = msprintf(gettext("%s: values of order and delay matrix [nb nc nd nf nk] should be nonnegative integer number "), "estpoly"); + error(errmsg); + end + + na = n(1); nb = n(2); nc = n(3); nd = n(4);nf = n(5); + + if (size(n,"*") == 5) then + nk = 1 + else + nk = n(6); + end + + // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + YDATA = z(:,1); + UDATA = z(:,2); + NDATA = size(UDATA,"*"); + function e = G(p,m) + e = YDATA - _oestpolyfun(UDATA,p,na,nb,nc,nd,nf,nk);//_oestpolyfun(UDATA,p,nd,nc,nf,nb,nk); + endfunction + tempSum = na+nb+nc+nd+nf + p0 = linspace(0.0001,0.001,tempSum)'; + [var,errl] = lsqrsolve(p0,G,size(UDATA,"*")); + //disp(errl) + err = (norm(errl)^2); + //disp(err) + opt_err = err; + resid = G(var,[]); + x = var +// b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff"); +// c = poly([1; var(nb+1:nb+nc)]',"q","coeff"); +// d = poly([1; var(nb+nc+1:nb+nc+nd)]',"q","coeff"); +// f = poly([1; var(nb+nd+nc+1:nd+nc+nf+nb)]',"q","coeff"); + a = poly([1; x(1:na)]',"q","coeff"); + b = poly([repmat(0,nk,1);x(na+1:na+nb)]',"q","coeff"); + c = poly([1; x(na+nb+1:na+nb+nc)]',"q","coeff"); + d = poly([1; x(na+nb+nc+1:na+nb+nc+nd)]',"q","coeff"); + f = poly([1; x(na+nb+nd+nc+1:na+nd+nc+nf+nb)]',"q","coeff"); + t = idpoly(coeff(a),coeff(b),coeff(c),coeff(d),coeff(f),Ts) + + //t = sys;//idpoly(1,coeff(b),coeff(c),coeff(d),coeff(f),Ts) + + // estimating the other parameters + [temp1,temp2,temp3] = predict(z,t) + [temp11,temp22,temp33] = pe(z,t) + //pause + estData = calModelPara(temp1,temp11,na+nb+nc+nd+nf) + //pause + t.Report.Fit.MSE = estData.MSE + t.Report.Fit.FPE = estData.FPE + t.Report.Fit.FitPer = estData.FitPer + t.Report.Fit.AIC = estData.AIC + t.Report.Fit.AICc = estData.AICc + t.Report.Fit.nAIC = estData.nAIC + t.Report.Fit.BIC = estData.BIC + t.TimeUnit = unit + sys = t + + //sys.TimeUnit = unit +endfunction + +function yhat = _oestpolyfun(UDATA,x,na,nb,nc,nd,nf,nk)//(UDATA,x,nd,nc,nf,nb,nk) + x=x(:) + q = poly(0,'q') + tempSum = na+nb+nc+nd+nf + // making polynomials + a = poly([1; x(1:na)]',"q","coeff"); + b = poly([repmat(0,nk,1);x(na+1:na+nb)]',"q","coeff"); + c = poly([1; x(na+nb+1:na+nb+nc)]',"q","coeff"); + d = poly([1; x(na+nb+nc+1:na+nb+nc+nd)]',"q","coeff"); + f = poly([1; x(na+nb+nd+nc+1:na+nd+nc+nf+nb)]',"q","coeff"); + bd = coeff(b*d); cf = coeff(c*f); fc_d = coeff(f*(c-a*d)); + if size(bd,"*") == 1 then + bd = repmat(0,nb+nd+1,1) + end + //pause + maxDelay = max([length(bd) length(cf) length(fc_d)]) + yhat = [YDATA(1:maxDelay)] + for k=maxDelay+1:size(UDATA,"*") + bdadd = 0 + for i = 1:size(bd,"*") + bdadd = bdadd + bd(i)*UDATA(k-i+1) + end + + fc_dadd = 0 + for i = 1:size(fc_d,"*") + fc_dadd = fc_dadd + fc_d(i)*YDATA(k-i+1) + end + cfadd = 0 + for i = 2:size(cf,"*") + cfadd = cfadd + cf(i)*yhat(k-i+1) + end + //pause + yhat = [yhat; [ bdadd + fc_dadd - cfadd ]]; + end +endfunction diff --git a/etfe.sci b/etfe.sci new file mode 100644 index 0000000..b09b7c4 --- /dev/null +++ b/etfe.sci @@ -0,0 +1,29 @@ +function varargout = etfe(varargin) + [lhs,rhs] = argn() + data = varargin(1) + if rhs == 1 then + n = 128 + elseif rhs == 2 then + n = varargin(2) + end + y = data.OutputData; + u = data.InputData + N = size(y,'r') + v = linspace(1,N,n) + y1 = y(v);u1 = u(v) + //y1 = y((1:(N-1)/(n-1):N));u1 = u((1:(N-1)/(n-1):N)) + //y1($) = y(N);u1($) = u(N) + data12 = [y1,u1] + z = [fft(y1),fft(u1)] + z = z/size(z,'r') + magData1 = abs(z(:,1));magData2 = abs(z(:,2)) + argData1 = phasemag(z(:,1),'m');argData2 = phasemag(z(:,2),'m') + magData = magData1./magData2;argData = argData1-argData2 + argData = [cosd(argData) sind(argData)] + data = [magData.*argData(:,1) magData.*argData(:,2)] + output = data(:,1)+%i*data(:,2) + resp = output(1:ceil(length(output)/2)) + frq = (1: ceil(n/2)) * %pi/floor(n/2) + output = frd(frq,resp,1) + varargout(1)= output +endfunction diff --git a/etfeeTest.sce b/etfeeTest.sce new file mode 100644 index 0000000..46f250c --- /dev/null +++ b/etfeeTest.sce @@ -0,0 +1,13 @@ +loadmatfile('rData.mat') +z = [fft(data1(:,1)),fft(data1(:,2))] +z = z/size(z,'r') +magData1 = abs(z(:,1));magData2 = abs(z(:,2)) +argData1 = phasemag(z(:,1),'m');argData2 = phasemag(z(:,2),'m') +magData = magData1./magData2;argData = argData1-argData2 +argData = [cosd(argData) sind(argData)] +data = [magData.*argData(:,1) magData.*argData(:,2)] +output = data(:,1)+%i*data(:,2) +resp = output(1:ceil(length(output)/2)) +frq = (1: ceil(256/2)) * %pi/floor(256/2) +output = frd(frq,resp,1) +disp(output) diff --git a/fitch.sci b/fitch.sci new file mode 100644 index 0000000..b72b472 --- /dev/null +++ b/fitch.sci @@ -0,0 +1,42 @@ +function varargout = fitch(varargin) + [lhs , rhs] = argn() + if ( rhs <> 1 ) then + errmsg = msprintf(gettext("%s: Wrong number of input arguments"), "fitch"); + error(errmsg) + elseif typeof(varargin(1)) <> "idpoly" then + error(msprintf(gettext("%s:Input model must be ""idpoly"" type.\n"),"fitch")) + end + model = varargin(1) + MSE = model.Report.Fit.MSE + FPE = model.Report.Fit.FPE + FitPer = model.Report.Fit.FitPer + AIC = model.Report.Fit.AIC + AICc = model.Report.Fit.AICc + nAIC = model.Report.Fit.nAIC + BIC = model.Report.Fit.BIC + t = tlist(['fitch','MSE','FPE','FitPer','AIC','AICc','nAIC','BIC'],MSE,FPE,FitPer,AIC,AICc,nAIC,BIC) + varargout(1) = t +endfunction + +function %fitch_p(mytlist) + f = fieldnames(mytlist) + maxLength = [] + for ii = 1:size(f,'*') + maxLength = [maxLength length(f(ii))] + end + maxLength = max(maxLength) + for ii = 1:size(f,'*') + blanckSpace = ' ' + for jj = 1:maxLength-length(f(ii)) + blanckSpace = blanckSpace + ' ' + end + mprintf('\t%s%s : ',blanckSpace,f(ii)) + objectData = mytlist(f(ii)) + if ceil(objectData)-objectData then + mprintf("%.4f",objectData) + else + mprintf("%d",objectData) + end + mprintf("\n") + end +endfunction diff --git a/frd.sci b/frd.sci new file mode 100644 index 0000000..c3838b4 --- /dev/null +++ b/frd.sci @@ -0,0 +1,84 @@ +function varargout = frd(varargin) + [lhs,rhs] = argn(0) + if rhs < 2 || rhs > 4 then + errmsg = msprintf(gettext("%s: Wrong numbers of input arguments."), "frd"); + error(errmsg) + end + frequency = varargin(2) + freqUnit = 'rad/TimeUnit' + if ~iscolumn(frequency) then + errmsg = msprintf(gettext("%s: frequency must be a finite row vector."), "frd"); + error(errmsg) + end + respData = varargin(1) +// pause + if size(frequency,'r') <> size(respData,'r') then + errmsg = msprintf(gettext("%s: input output matrix dimension must be equal."), "frd"); + error(errmsg) + end + if rhs == 2 then + Ts = 0 + elseif rhs >2 then + Ts = varargin(3) + end + if Ts < 0 || size(Ts,'*') <> 1 || typeof(Ts) <> 'constant' then + errmsg = msprintf(gettext("%s: Sampling time must be a scalar non negative real number."), "frd"); + error(errmsg) + end + // saving the spectrum value + if rhs == 4 then + spect = varargin(4) + else + spect = [] + end + /// matching its dimensions + if ~size(spect) then + elseif size(frequency,'r') <> size(spect,'r') then + errmsg = msprintf(gettext("%s: Numbers of power spectra must be equal to the numbers of frequency."), "frd"); + error(errmsg) + end + TUnit = 'seconds' + t = tlist(['frd','Frequency','FrequencyUnit','ResponseData','Ts','TimeUnit','Spect'],frequency,freqUnit,respData,Ts,TUnit,spect) + varargout(1) = t +endfunction + +//overloading +function %frd_p(varargin) + myTlist = varargin(1) + f = fieldnames(myTlist) + freqData = myTlist.Frequency + tempRespData= myTlist.ResponseData + for jj = 1:size(tempRespData,'c') + respData = tempRespData(:,jj) + mprintf("\t -------------------------") + mprintf("\n") + mprintf("\t Frequency \t Response") + mprintf("\n") + mprintf("\t -------------------------") + mprintf("\n") + for ii = 1:length(myTlist.Frequency) + temp = '' + if real(respData(ii))>=0 then + temp = temp + ' ' + end + temp = temp + string(real(respData(ii))) + // temp = string(real(respData(ii))) + if imag(respData(ii)) > 0 then + temp = temp +"+" + end + if ~imag(respData(ii)) then + else + temp = temp + string(imag(respData(ii))) +"i" + end + // temp = temp + string(imag(respData(ii))) + " i" + mprintf("\n\t %f \t %s",freqData(ii),temp)//real(respData(ii)),imag(respData(ii))) + end + mprintf("\n\n") + end + if ~myTlist.Ts then + mprintf("\n Continuous Domain frequency response.") + else + mprintf("\n Sampling Time = "+string(myTlist.Ts)+" "+myTlist.TimeUnit) + mprintf("\n Discrete Domain frequency response.") + end +endfunction diff --git a/frdplot.sci b/frdplot.sci new file mode 100644 index 0000000..c923429 --- /dev/null +++ b/frdplot.sci @@ -0,0 +1,11 @@ +function frdplot(varargin) + [lhs rhs] = argn(0) + if rhs <> 1 then + error(msprintf(gettext("%s: Wrong number of input arguments."),"frdplot")) + end + frdData = varargin(1) + if typeof(frdData) <> 'frd' then + error(msprintf(gettext("%s:Wrong type for input argument %d: ""frd"" expected."),"frdplot",1)) + end + bode((frdData.Frequency)',(frdData.ResponseData)') +endfunction diff --git a/iddata.sci b/iddata.sci new file mode 100755 index 0000000..e748029 --- /dev/null +++ b/iddata.sci @@ -0,0 +1,68 @@ +function varargout = iddata(varargin) + [lhs,rhs] = argn(0) + if rhs == 2 | rhs == 1 then + Ts = 1 + elseif rhs == 3 then + Ts = varargin(3) + else + error(msprintf(gettext("%s:Incorrect number of input arguments.\n"),"iddata")) + end + if Ts <= 0 || typeof(Ts) <> 'constant' || ~size(Ts,'*') then + error(msprintf(gettext("%s:Inconsist sampling time ""Ts"", ""non negative real number"" expected.\n"),"iddata")) + end + if rhs == 1 then + OutputData = varargin(1) + InputData = [] + elseif rhs == 2 | rhs == 3 then + OutputData = varargin(1) + InputData = varargin(2) + if isrow(InputData) then + InputData = InputData' + end + end + if size(OutputData,'*') & size(InputData,'*') then + if size(OutputData,'r') ~= size(InputData,'r') then + error(msprintf(gettext("%s:The numbers of the plant out datas must be equal to the numbers of the plant input datas.\n"),"iddata")) + end + end + t = tlist(['iddata','OutputData','InputData','Ts','TimeUnit'],OutputData,InputData,Ts,'seconds') + varargout(1) = t +endfunction + +function %iddata_p(mytlist) + f = fieldnames(mytlist) + if ~size(mytlist(f(1)),'*') & ~size(mytlist(f(2)),'*') then + mprintf(' Empty sample data.\n') + else + outputSize = size(mytlist(f(1))) + inputSize = size(mytlist(f(2))) + if prod(outputSize) then + sampleSize = max(outputSize) + elseif prod(inputSize) then + sampleSize = max(inputSize) + end + mprintf(' Time domain sample data having %d samples.',sampleSize) + if round(mytlist(f(3)))-mytlist(f(3)) == 0 then + mprintf('\n Sampling Time = %d',mytlist(f(3))) + else + mprintf('\n Sampling Time = %f',mytlist(f(3))) + end + mprintf(' %s',mytlist(f(4))) + mprintf('\n') + if prod(outputSize) then + mprintf('\n Output channel \n') + for ii = 1:min(outputSize) + yString = 'y' + string(ii) + mprintf(' %s\n',yString) + end + end + if prod(inputSize) then + mprintf('\n Input channel \n') + for ii = 1:min(inputSize) + uString = 'u' + string(ii) + mprintf(' %s\n',uString) + end + end + end + mprintf('\n') +endfunction diff --git a/iddataplot.sci b/iddataplot.sci new file mode 100644 index 0000000..4b1818d --- /dev/null +++ b/iddataplot.sci @@ -0,0 +1,69 @@ +function iddataplot(varargin) + [lhs rhs] = argn(0) + if rhs <> 1 then + error(msprintf(gettext("%s: Wrong number of input arguments."),"iddataplot")) + end + iddataData = varargin(1) + if typeof(iddataData) <> 'iddata' then + error(msprintf(gettext("%s:Wrong type for input argument %d: ""iddata"" expected."),"iddataplot",1)) + end +// figure() +// xtitle('Plant i/o Data','Time('+iddataData.TimeUnit+')','Amplitude') + uData = iddataData.InputData + yData = iddataData.OutputData + timeLength = max(size(uData,'r'),size(yData,'r')) + timeData = ((1:1:timeLength)*iddataData.Ts)' + //timeData = timeData(1:length(iddataData)-1) + if size(uData,'*') && size(yData,'*') then + firstIndex = 2 + else + firstIndex = 1 + end + // ploting y data + h = gcf() + if h.figure_name == "Plant Input Output Data"; then + clf() + end + if size(yData,'*') then + secondIndex = size(yData,'c') + if secondIndex == 1 then + outputString = 'y' + else + outputString = [] + for ii = 1:secondIndex + outputString = [outputString 'y'+string(ii)] + end + end + //disp(outputString) + for ii = 1:secondIndex + subplot(firstIndex,secondIndex,ii);plot(timeData,yData(:,ii)) + xtitle(outputString(ii)) + end + end + if size(uData,'*') then + secondIndex = size(uData,'c') + // disp(secondIndex) + if secondIndex == 1 then + outputString = 'u' + else + outputString = [] + for ii = 1:secondIndex + outputString = [outputString 'u'+string(ii)] + end + end + //disp(outputString) + for ii = 1:secondIndex + if size(yData,'*') then + temp = 1 + else + temp = 0 + end + subplot(firstIndex,secondIndex,ii+secondIndex*temp);plot2d(timeData,uData(:,ii),2) + //pause + xtitle(outputString(ii)) + end + end + h = gcf() + //disp(h) + h.figure_name= "Plant Input Output Data"; +endfunction diff --git a/identTime.sci b/identTime.sci new file mode 100644 index 0000000..c4c47c8 --- /dev/null +++ b/identTime.sci @@ -0,0 +1,32 @@ +function varargout = identTime(varargin) + [lhs , rhs] = argn() + if ( rhs <> 1 ) then + errmsg = msprintf(gettext("%s: Wrong number of input arguments"), "identTime"); + error(errmsg) + elseif typeof(varargin(1)) <> "iddata" then + error(msprintf(gettext("%s:Input model must be ""iddata"" type.\n"),"identTime")) + end + plantData = varargin(1) + //disp('yolo') + inputData = plantData.InputData;inputData = size(inputData,'r') + outputData = plantData.OutputData;outputData = size(outputData,'r') + sampleNumb = max(inputData,outputData) + timeData = (0:sampleNumb-1)*plantData.Ts + t = tlist(['identTime','samples','start','end','Frequency','TimeSeries'],sampleNumb,0,timeData($),1/plantData.Ts,timeData) + varargout(1) = t +endfunction + +function %identTime_p(mytlist) + f = fieldnames(mytlist) + mprintf("\t samples : %d\n",mytlist.samples) + mprintf("\t start : %d\n",mytlist.start) + mprintf("\t end : %d\n",mytlist.end) + if ceil(mytlist.Frequency)-mytlist.Frequency then + mprintf("\t Frequency : %.4f\n",mytlist.Frequency) + else + mprintf("\t Frequency : %d\n",mytlist.Frequency) + end + timeData = mytlist.TimeSeries + mprintf("\t TimeSeries : %.2f, %.2f, . ,%.2f",timeData(1),timeData(2),timeData($)) + +endfunction diff --git a/idpoly.sci b/idpoly.sci new file mode 100755 index 0000000..0be3754 --- /dev/null +++ b/idpoly.sci @@ -0,0 +1,142 @@ +function sys = idpoly(varargin) + [lhs,rhs] = argn(0) + tempCell = cell(6,1) + tempTs = 0 + for ii = 1:rhs + if typeof(varargin(ii)) == 'string' & varargin(ii) == 'Ts' then + tempCell{6,1} = varargin(ii+1) + break + elseif ~size(varargin(ii),'*') then + tempCell{ii,1} = 1 + elseif size(varargin(ii),'*') then + tempCell{ii,1} = varargin(ii) + end + end + //disp(tempCell) + for ii = 1:6 + if ~size(cell2mat(tempCell(ii,1)),'*') & ii ~= 6 then + tempCell{ii,1} = 1 + elseif ~size(cell2mat(tempCell(ii,1)),'*') & ii == 6 then + tempCell{ii,1} = -1 + end + end + //storing the data in A,B,C,D,F matrix + // B(z) D(z) + // y(n)=---------- u(n)+ ----------- e(n) + // A(z)F(z) A(z)D(z) + A = cell2mat(tempCell(1,1)); B = cell2mat(tempCell(2,1)); + C = cell2mat(tempCell(3,1)); D = cell2mat(tempCell(4,1)); + F = cell2mat(tempCell(5,1)); Ts = cell2mat(tempCell(6,1)); + if A(1,1) ~= 1 then + error(msprintf(gettext("%s: The first coefficient of ""A(z)"" polynomial must be 1.\n"),"idpoly")) + elseif C(1,1) ~= 1 then + error(msprintf(gettext("%s: The first coefficient of ""C(z)"" polynomial must be 1.\n"),"idpoly")) + elseif D(1,1) ~= 1 then + error(msprintf(gettext("%s: The first coefficient of ""D(z)"" polynomial must be 1.\n"),"idpoly")) + elseif F(1,1) ~= 1 then + error(msprintf(gettext("%s: The first coefficient of ""F(z)"" polynomial must be 1.\n"),"idpoly")) + end + report = struct('MSE',0,'FPE',0,'FitPer',0,'AIC',0,'AICc',0,'nAIC',0,'BIC',0) + errors = [0 0 0] + report = struct('Fit',report,'Uncertainty',errors) + t = tlist(['idpoly','a','b','c','d','f','Variable','TimeUnit','Ts','Report'],A,B,C,D,F,'z^-1','seconds',Ts,report) + //t = tlist(['idpoly','a','b','c','d','f','Variable','TimeUnit','Ts'],A,B,C,D,F,'z^-1','seconds',Ts) + + sys = t +endfunction + + +function %idpoly_p(mytlist) + f = fieldnames(mytlist) + //A polynomial + if mytlist(f(1)) == 1 && size(mytlist(f(1)),'*') == 1 then + else + mprintf('\n A(z) =') + temp = poly2str(mytlist(f(1))) + mprintf('%s\n',temp) + end + + //B polynomial + if mytlist(f(2)) == 1 then + else + mprintf('\n B(z) =') + temp = poly2str(mytlist(f(2))) + mprintf('%s\n',temp) + end + + //C polynomial + if mytlist(f(3)) == 1 && size(mytlist(f(3)),'*') == 1 then + else + mprintf('\n C(z) =') + temp = poly2str(mytlist(f(3))) + mprintf('%s\n',temp) + end + //D polynomial + if mytlist(f(4)) == 1 && size(mytlist(f(4)),'*') == 1 then + elseif size(mytlist(f(4)),'*') > 1 then + mprintf('\n D(z) =') + temp = poly2str(mytlist(f(4))) + mprintf('%s\n',temp) + end + + //F polynomial + if mytlist(f(5)) == 1 && size(mytlist(f(5)),'*') == 1 then + else + mprintf('\n F(z) =') + temp = poly2str(mytlist(f(5))) + mprintf('%s\n',temp) + end + + mprintf('\n Sampling Time = ') + + if mytlist.Ts == -1 then + mprintf('undefined') + else + if (ceil(mytlist.Ts)-mytlist.Ts) == 0 then + mprintf('%d %s',mytlist.Ts,mytlist.TimeUnit) + else + mprintf('%0.4f %s',mytlist.Ts,mytlist.TimeUnit) + end + end + //disp(mytlist.Ts) + mprintf('\n') + if mytlist.Report.Fit.MSE then + temp = ['MSE','FPE','FitPer','AIC','AICc','nAIC','BIC'] + spaces = ' ' + for ii = 1:size(temp,'c') + digiLength = length(string(round(mytlist.Report.Fit(temp(ii))))) + digiLength = digiLength + 5-length(temp(ii)) + blank = '' + for jj = 1:digiLength+1 + blank = blank + " " + end + spaces = spaces+blank+temp(ii)+' ' + end + mprintf('\n') + mprintf(spaces) + mprintf("\n %.4f %.4f %.4f %.4f %.4f %.4f %.4f",mytlist.Report.Fit.MSE,mytlist.Report.Fit.FPE,mytlist.Report.Fit.FitPer,mytlist.Report.Fit.AIC,mytlist.Report.Fit.AICc,mytlist.Report.Fit.nAIC,mytlist.Report.Fit.BIC) + end + + +endfunction + + +function strout = poly2str(h) + temp = poly(h,'x','coeff') + temp = pol2str(temp) + temp = strsubst(temp,'-',' - ') + temp = strsubst(temp,'x^',' z^-') + temp = strsubst(temp,'x',' z^-1') + temp = strsubst(temp,'*','') + temp = strsubst(temp,'+',' + ') + [ind which]=strindex(temp,'-') + //disp(ind) +// disp(which) + if ind(1,1) ~= 2 then + temp = ' ' + temp + elseif ind(1,1) == 2 then + temp = part(temp,4:length(temp)) + temp = ' -' + temp + end + strout = temp +endfunction diff --git a/impulseest.sci b/impulseest.sci new file mode 100644 index 0000000..8817bcc --- /dev/null +++ b/impulseest.sci @@ -0,0 +1,78 @@ +function varargout = impulseest(varargin) + [lhs,rhs] = argn(0) + //checking the number of inputs + if rhs > 2 then + error(msprintf(gettext("%s: Unexpected number of input arguments "),"impulseest")) + end + modelData = varargin(1) + if typeof(modelData) <> "idpoly" then + error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"impulseest")) + end + //adding noise + if rhs == 2 then + noiseFlag = varargin(2) + if typeof(noiseFlag) <> 'boolean' then + error(msprintf(gettext("%s: Last input data must be ""boolean"".type "),"impulseest")) + end + else + noiseFlag = %F + end + z = poly(0,'z') + aPoly = poly(modelData.a(length(modelData.a):-1:1),'z','coeff') + bPoly = poly(modelData.b,'z','coeff') + fPoly = poly(modelData.f(length(modelData.f):-1:1),'z','coeff') + afPoly = aPoly*fPoly + + bCoeff = modelData.b + extra = 1 + if ~bCoeff(1,1) then + afCoeff = coeff(afPoly) + bLength = length(bCoeff);afLength = length(afCoeff) + if bLength == afLength then + extra = 1 + else + extra = z^-(bLength-afLength) + end + end + + bPoly = poly(modelData.b(length(modelData.b):-1:1),'z','coeff') + sys = syslin('d',bPoly,afPoly)*extra + if size(find(gsort(abs(roots(sys.den)))>1),'*') then + tempRoot = clean(roots(sys.den)) + [sorted index] = gsort(abs(real(tempRoot))) + tempRoot = tempRoot(index) + n = 26/log10(abs(real(tempRoot(1)))) + + else + finalValue = horner(sys,1)*1 + n = 10 + tempValue = sum(ldiv(sys.num,sys.den,n)) + if finalValue >= 0 then + while finalValue*0.99 >= tempValue || n >=1000 + n = n+5 + tempValue = sum(ldiv(sys.num,sys.den,n)) + end + elseif finalValue < 0 then + while finalValue*0.99 =1000 + n = n+5 + tempValue = sum(ldiv(sys.num,sys.den,n)) + end + if n > 100 then + n = 100 + end + end + end + uData = [1 zeros(1,n)] + yData = flts(uData,sys) + timeData = (0:(n))*modelData.Ts +// pause + if noiseFlag then + varargout(1) = yData' + else + stem(timeData,yData) + h = gcf() + + h.figure_name= "Plant Impulse Response" + varargout(1) = [] + end +endfunction diff --git a/iv.sci b/iv.sci new file mode 100644 index 0000000..a819a4e --- /dev/null +++ b/iv.sci @@ -0,0 +1,98 @@ +function varargout = iv(varargin) + + [lhs , rhs] = argn(0); + if ( rhs < 2 || rhs > 3) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2 or 3"), "iv", rhs); + error(errmsg) + end + plantData = varargin(1) + if typeof(plantData) == 'iddata' then + Ts = plantData.Ts;unit = plantData.TimeUnit + plantData = [plantData.OutputData plantData.InputData] + elseif typeof(plantData) == 'constant' then + Ts = 1;unit = 'seconds' + end + if ((~size(plantData,2)==2) & (~size(plantData,1)==2)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "iv"); + error(errmsg); + end + + if (~isreal(plantData)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "arx"); + error(errmsg); + end + + n = varargin(2) + if (size(n,"*")<2| size(n,"*")>3) then + errmsg = msprintf(gettext("%s: The order and delay matrix [na nb nk] should be of size [2 or 3]"), "iv"); + error(errmsg); + end + + if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then + errmsg = msprintf(gettext("%s: values of order and delay matrix [na nb nk] should be nonnegative integer number "), "iv"); + error(errmsg); + end + na = n(1);nb = n(2) + if size(n,'*') == 2 then + nk = 1 + elseif size(n,'*') == 3 then + nk = n(3) + end + yData = plantData(:,1) + uData = plantData(:,2) + noOfSample = size(plantData,'r') + nb1 = nb + nk - 1 + n = max(na, nb1) + if rhs == 3 then + if typeof(varargin(3)) <> 'constant' then + errmsg = msprintf(gettext("%s: Incompatible last input argument. "), "iv"); + error(errmsg) + elseif size(varargin(3),'r') <> size(uData,'r') then + errmsg = msprintf(gettext("%s: number of samples of output must be equal to the dimensions of plant data "), "iv"); + error(errmsg); + end + x = varargin(3) + elseif rhs == 2 + arxModel = arx(plantData,[na nb nk]) + x = sim(uData,arxModel) + end + phif = zeros(noOfSample,na+nb) + psif = zeros(noOfSample,na+nb) + // arranging samples of y matrix + for ii = 1:na + phif(ii+1:ii+noOfSample,ii) = yData + psif(ii+1:ii+noOfSample,ii) = x + end + // arranging samples of u matrix + for ii = 1:nb + phif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData + psif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData + end + lhs = psif'*phif + lhsinv = pinv(lhs) + //pause + theta = lhsinv * (psif)' * [yData;zeros(n,1)] + ypred = (phif * theta) + ypred = ypred(1:size(yData,'r')) + e = yData - ypred + sigma2 = norm(e)^2/(size(yData,'r') - na - nb) + vcov = sigma2 * pinv((phif)' * phif) + t = idpoly([1; -theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,1) + + // estimating the other parameters + [temp1,temp2,temp3] = predict(z,t) + [temp11,temp22,temp33] = pe(z,t) + + estData = calModelPara(temp1,temp11,na+nb) + //pause + t.Report.Fit.MSE = estData.MSE + t.Report.Fit.FPE = estData.FPE + t.Report.Fit.FitPer = estData.FitPer + t.Report.Fit.AIC = estData.AIC + t.Report.Fit.AICc = estData.AICc + t.Report.Fit.nAIC = estData.nAIC + t.Report.Fit.BIC = estData.BIC + t.TimeUnit = unit + //sys = t + varargout(1) = t +endfunction diff --git a/iv4.sci b/iv4.sci new file mode 100644 index 0000000..32fd778 --- /dev/null +++ b/iv4.sci @@ -0,0 +1,75 @@ +function varargout = iv4(varargin) + [lhs, rhs] = argn(0) + + plantData = varargin(1) + orderData = varargin(2) + na = orderData(1);nb = orderData(2) + // arranging na ,nb,nk + if size(orderData,"*") == 2 then + nk = 1 + elseif size(orderData,'*') == 3 then + nk = orderData(3) + end + nb1 = nb + nk - 1 + n = max(na, nb1) + // arranging the plant data + if typeof(plantData) == 'constant' then + Ts = 1;unitData = 'second' + elseif typeof(plantData) == 'iddata' then + Ts = plantData.Ts;unitData = plantData.TimeUnit + plantData = [plantData.OutputData plantData.InputData] + end + noOfSample = size(plantData,'r') + // finding the iv model + ivTest = iv(plantData,[na nb nk]); + // residual + [aTemp,bTemp,cTemp] = pe(plantData,ivTest); + Lhat = ar(aTemp,na+nb); + x = sim(plantData(:,2),ivTest); + yData = plantData(:,1);uData = plantData(:,2) + Yf = filter(Lhat.a,Lhat.b,[plantData(:,1);zeros(n,1)]); + phif = zeros(noOfSample,na+nb) + psif = zeros(noOfSample,na+nb) + // arranging samples of y matrix + for ii = 1:na + phif(ii+1:ii+noOfSample,ii) = -yData + psif(ii+1:ii+noOfSample,ii) = -x + end + // arranging samples of u matrix + for ii = 1:nb + phif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData + psif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData + end + // passing it through the filters + for ii = 1:na+nb + phif(:,ii) = filter(Lhat.a,Lhat.b,phif(:,ii)); + psif(:,ii) = filter(Lhat.a,Lhat.b,psif(:,ii)); + end + lhs = psif'*phif + lhsinv = pinv(lhs) + theta = lhsinv * (psif)' * Yf + ypred = (phif * theta) + ypred = ypred(1:size(yData,'r')) + e = yData - ypred + sigma2 = norm(e)^2/(size(yData,'r') - na - nb) + vcov = sigma2 * pinv((phif)' * phif) + t = idpoly([1; theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,Ts) + + // estimating the other parameters + [temp1,temp2,temp3] = predict(z,t) + [temp11,temp22,temp33] = pe(z,t) + + estData = calModelPara(temp1,temp11,na+nb) + //pause + t.Report.Fit.MSE = estData.MSE + t.Report.Fit.FPE = estData.FPE + t.Report.Fit.FitPer = estData.FitPer + t.Report.Fit.AIC = estData.AIC + t.Report.Fit.AICc = estData.AICc + t.Report.Fit.nAIC = estData.nAIC + t.Report.Fit.BIC = estData.BIC + t.TimeUnit = unitData + //sys = t + varargout(1) = t + //varargout(1) = idpoly([1; -theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,Ts) +endfunction diff --git a/misData.mat b/misData.mat new file mode 100644 index 0000000..6d26b91 Binary files /dev/null and b/misData.mat differ diff --git a/misdata.sci b/misdata.sci new file mode 100644 index 0000000..c1ced50 --- /dev/null +++ b/misdata.sci @@ -0,0 +1,41 @@ +function varargout = misdata(varargin) + [lhs,rhs] = argn(0) +//------------------------------------------------------------------------------ +// checking the number of inputs + if rhs <> 1 then + error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"misdata")) + end +//------------------------------------------------------------------------------ + ioData = varargin(1) + if typeof(ioData) <> "iddata" then + error(msprintf(gettext("%s: Plant input data must be ""iddata"" type. "),"misdata")) + end + inputMat = ioData.InputData;inputMat = linearINTRP(inputMat,abs(ioData.Ts));ioData.InputData = inputMat; + outputMat = ioData.OutputData;outputMat = linearINTRP(outputMat,abs(ioData.Ts));ioData.OutputData = outputMat; + varargout(1) = ioData +endfunction + +function varargout = linearINTRP(matData,Ts) + // looking for overall nan values + nanData = isnan(matData);nanIndex = find(nanData == %T) + if ~size(nanIndex,'*') then + varargout(1) = matData + else + tempMat = [] + matSize = size(matData,'r') + // looking for nan in each column + for ii = 1:size(matData,'c') + nanData = isnan(matData(:,ii));nanIndex = find(nanData == %T); + if ~size(nanData,'*') then + tempMat = [tempMat matData(,ii)] + else + timeData = (linspace(1*Ts,matSize*Ts , matSize))'; + nanMat = isnan(matData(:,ii)); + data = matData(:,ii) + data(nanMat) = interp1(timeData(~nanMat), data(~nanMat), timeData(nanMat)); + tempMat = [tempMat data] + end + end + varargout(1) = tempMat + end +endfunction diff --git a/nInputData.sci b/nInputData.sci new file mode 100644 index 0000000..59a3b64 --- /dev/null +++ b/nInputData.sci @@ -0,0 +1,15 @@ +function varargout = nInputSeries(varargin) + [lhs rhs] = argn(0) + if rhs <> 1 then + error(msprintf(gettext("%s: Wrong number of input arguments."),"iddataplot")) + end + iddataData = varargin(1) + if typeof(iddataData) <> 'iddata' then + error(msprintf(gettext("%s:Wrong type for input argument %d: ""iddata"" expected."),"nInputSeries",1)) + end + if ~size(iddataData.InputData,'*') then + varargout(1) = 1 + else + varargout(1) = size(iddataData.InputData,'c') + end +endfunction diff --git a/oe.sci b/oe.sci new file mode 100644 index 0000000..8e249ec --- /dev/null +++ b/oe.sci @@ -0,0 +1,118 @@ +// Copyright (C) 2015 - IIT Bombay - FOSSEE +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Authors: Harpreet, Ashutosh +// Organization: FOSSEE, IIT Bombay + +function sys = oe(varargin) + +// Estimates Discrete time BJ model +// y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t) +// Current version uses random initial guess +// Need to get appropriate guess from OE and noise models + + [lhs , rhs] = argn(); + if ( rhs < 2 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "oe_2", rhs); + error(errmsg) + end + + + z = varargin(1) + if typeof(z) == 'iddata' then + Ts = z.Ts;unit = z.TimeUnit + z = [z.OutputData z.InputData] + elseif typeof(z) == 'constant' then + Ts = 1;unit = 'seconds' + end + if ((~size(z,2)==2) & (~size(z,1)==2)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "oe_2"); + error(errmsg); + end + + if (~isreal(z)) then + errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "oe_2"); + error(errmsg); + end +// + n = varargin(2) + if (size(n,"*")<2| size(n,"*")>3) then + errmsg = msprintf(gettext("%s: The order and delay matrix [nb nf nk] should be of size [2 4]"), "oe_2"); + error(errmsg); + end + + if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then + errmsg = msprintf(gettext("%s: values of order and delay matrix [nb nf nk] should be nonnegative integer number "), "oe_2"); + error(errmsg); + end +// + nb= n(1); nf = n(2); +// + if (size(n,"*") == 2) then + nk = 1 + else + nk = n(3); + end + + // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively + YDATA = z(:,1); + UDATA = z(:,2); + NDATA = size(UDATA,"*"); + function e = G(p,m) + e = YDATA - _objoefun(UDATA,YDATA,p,nf,nb,nk); + endfunction + tempSum = nf+nb + p0 = linspace(0.04,0.041,tempSum)'; + [var,errl] = lsqrsolve(p0,G,size(UDATA,"*")); + err = (norm(errl)^2); + opt_err = err; + resid = G(var,[]); + f = poly([1; var(nb+1:nb+nf)],"q","coeff"); + b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff"); + t = idpoly(1,coeff(b),1,1,coeff(f),Ts) + + // estimating the other parameters + [temp1,temp2,temp3] = predict(z,t) + [temp11,temp22,temp33] = pe(z,t) + + estData = calModelPara(temp1,temp11,n(1)+n(2)) + //pause + t.Report.Fit.MSE = estData.MSE + t.Report.Fit.FPE = estData.FPE + t.Report.Fit.FitPer = estData.FitPer + t.Report.Fit.AIC = estData.AIC + t.Report.Fit.AICc = estData.AICc + t.Report.Fit.nAIC = estData.nAIC + t.Report.Fit.BIC = estData.BIC + t.TimeUnit = unit + sys = t + //sys = t + //sys.TimeUnit = unit +endfunction + +function yhat = _objoefun(UDATA,YDATA,x,nf,nb,nk) + x=x(:) + q = poly(0,'q') + tempSum = nb+nf + // making polynomials + b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff"); + f = -1*poly([x(nb+1:nb+nf)]',"q","coeff") + fSize = coeff(f);bSize = coeff(b) + maxDelay = max([length(fSize) length(bSize)]) + yhat = [YDATA(1:maxDelay)] + for k=maxDelay+1:size(UDATA,"*") + tempB = 0 + for ii = 1:size(bSize,'*') + tempB = tempB + bSize(ii)*UDATA(k-ii+1) + end + tempF = 0 + for ii = 1:size(fSize,"*") + tempF = tempF + fSize(ii)*yhat(k-ii) + end + yhat = [yhat; [ tempB + tempF ]]; + end +endfunction diff --git a/pe.sci b/pe.sci new file mode 100644 index 0000000..3a097a8 --- /dev/null +++ b/pe.sci @@ -0,0 +1,112 @@ +function varargout = pe(varargin) + [lhs,rhs] = argn(0) + + +//------------------------------------------------------------------------------ + data = varargin(1) + model = varargin(2) + if rhs == 3 then + kStep = varargin(3) + elseif rhs == 2 then + kStep = 1 + end +//------------------------------------------------------------------------------ +// k step analysis + if typeof(kStep) <> 'constant' || isnan(kStep) then + error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"pe")) + end +// if given k step is infinity or [] + if isinf(kStep) || ~size(kStep,'*') then + kStep = 1 + end +// checking the dimensions + if size(kStep,'*') <> 1 || (ceil(kStep)-kStep) then + error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"pe")) + end +//------------------------------------------------------------------------------ +// checking the plant model + if typeof(model) ~= 'idpoly' then + error(msprintf(gettext("%s:Plant model must be ""idpoly"" type.\n"),"pe")) + end + modelSampleTime = model.Ts + modelTimeUnit = model.TimeUnit +//------------------------------------------------------------------------------ +//checking the data type + if typeof(data) <> 'iddata' && typeof(data) <> 'constant' then + error(msprintf(gettext("%s:Sample data must be ""iddata"" type or ""n x 2"" matrix type.\n"),"pe")) + end +// checking the plant data + if typeof(data) == 'iddata' then + if ~size(data.OutputData,'*') || ~size(data.InputData,'*') then + error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"pe")) + end + plantSampleTime = data.Ts + plantTimeUnit = data.TimeUnit + data = [data.OutputData data.InputData] + //disp('iddata') + elseif typeof(data) == 'constant' then + if size(data,'c') ~= 2 then + error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"pe")) + end + plantSampleTime = model.Ts + plantTimeUnit = model.TimeUnit + end +//------------------------------------------------------------------------------ +// comparing the sampling time + if modelSampleTime-plantSampleTime <> 0 then + error(msprintf(gettext("%s:The sample time of the model and plant data must be equal.\n"),"pe")) + end +// Comparing the time units + if ~strcmp(modelTimeUnit,plantTimeUnit) then + else + error(msprintf(gettext("%s:Time unit of the model and plant data must be equal.\n"),"pe")) + end +//------------------------------------------------------------------------------ + sampleLength = size(data,'r') + // one step ahead prediction + [y1 ,x0] = predict(data,model) + errorData1 = data(:,1)-y1 + if kStep == 1 then + eCapData = errorData1 + elseif kStep > 1 then + aPoly = poly(model.a,'q','coeff'); + dPoly = poly(model.d,'q','coeff'); + adCoeff = coeff(aPoly*dPoly);adCoeff = adCoeff(length(adCoeff):-1:1); + adPoly = poly(adCoeff,'q','coeff') + cCoeff = model.c;cCoeff = cCoeff(length(cCoeff):-1:1); + cPoly = poly(cCoeff,'q','coeff') + hBar = clean((ldiv(cPoly,adPoly,kStep))',0.00001) + hBar = hBar(length(hBar):-1:1) + hBarLength = length(hBar) + errorData1 = [zeros(hBarLength,1);errorData1] + eCapData = [] + + for ii = 1:sampleLength+1 + eCapData = [eCapData; hBar*errorData1(ii:ii+hBarLength-1)] + end + end + + timeData = (modelSampleTime:modelSampleTime:(sampleLength)*modelSampleTime)' + pseudoData = size(eCapData,'r') + eCapData = eCapData(abs(pseudoData-sampleLength)+1:$) + //pause + if lhs == 1 then + clf() + plot(timeData,eCapData) + axisData = gca() + + tempTimeUnit = 'Time('+modelTimeUnit+')' + xtitle('Predicted Response',tempTimeUnit,'y') + xgrid + //pause + varargout(1) = 0 + + elseif lhs == 2 then + varargout(1) = eCapData + varargout(2) = 0 + elseif lhs == 3 then + varargout(1) = eCapData + varargout(2) = timeData + varargout(3) = 0 + end +endfunction diff --git a/predict.sci b/predict.sci new file mode 100644 index 0000000..fcd4b23 --- /dev/null +++ b/predict.sci @@ -0,0 +1,192 @@ +//[y,x0] = predict(data,idpoly,k) +//References +//Digital Control(11.1.2) by Kanna M.Moudgalya +//System Identification Theory for User Second Edition (3.2) by Lennart Ljung +// Code Aurthor - Ashutosh Kumar Bhargava + +function varargout = predict(varargin) + [lhs,rhs] = argn(0) +//------------------------------------------------------------------------------ +// checking the number of inputs + if rhs < 2 || rhs > 3 then + error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"predict")) + end +//------------------------------------------------------------------------------ + data = varargin(1) + model = varargin(2) + if rhs == 3 then + kStep = varargin(3) + elseif rhs == 2 then + kStep = 1 + end + +//------------------------------------------------------------------------------ + +// k step analysis + if typeof(kStep) <> 'constant' || isnan(kStep) then + error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"predict")) + end +// if given k step is infinity or [] + if isinf(kStep) || ~size(kStep,'*') then + kStep = 1 + end +// checking the dimensions + if size(kStep,'*') <> 1 || (ceil(kStep)-kStep) then + error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"predict")) + end +//------------------------------------------------------------------------------ +// checking the plant model + if typeof(model) ~= 'idpoly' then + error(msprintf(gettext("%s:Plant model must be ""idpoly"" type.\n"),"predict")) + end + modelSampleTime = model.Ts + modelTimeUnit = model.TimeUnit +//------------------------------------------------------------------------------ +//checking the data type + if typeof(data) <> 'iddata' && typeof(data) <> 'constant' then + error(msprintf(gettext("%s:Sample data must be ""iddata"" type or ""n x 2"" matrix type.\n"),"predict")) + end +// checking the plant data + if typeof(data) == 'iddata' then + if ~size(data.OutputData,'*') || ~size(data.InputData,'*') then + error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"predict")) + end + plantSampleTime = data.Ts + plantTimeUnit = data.TimeUnit + data = [data.OutputData data.InputData] + //disp('iddata') + elseif typeof(data) == 'constant' then + if size(data,'c') ~= 2 then + error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"predict")) + end + plantSampleTime = model.Ts + plantTimeUnit = model.TimeUnit + end +//------------------------------------------------------------------------------ +// comparing the sampling time + if modelSampleTime-plantSampleTime <> 0 then + error(msprintf(gettext("%s:The sample time of the model and plant data must be equal.\n"),"predict")) + end +// Comparing the time units + if ~strcmp(modelTimeUnit,plantTimeUnit) then + else + error(msprintf(gettext("%s:Time unit of the model and plant data must be equal.\n"),"predict")) + end +//------------------------------------------------------------------------------ +// ckecking the k step size. if it greater than number of sample size then the +// k step will become 1 + if kStep >= size(data,'r') then + kStep = 1 + end +//------------------------------------------------------------------------------ + //storing the plant data + // B(z) C(z) + // y(n) = ---------- u(n) + ---------- e(n) + // A(z)*F(z) A(z)*D(z) + aPoly = poly(model.a,'q','coeff'); + bPoly = poly(model.b,'q','coeff'); + cPoly = poly(model.c,'q','coeff'); + dPoly = poly(model.d,'q','coeff'); + fPoly = poly(model.f,'q','coeff'); + Gq = bPoly/(aPoly*fPoly) + Hq = cPoly/(aPoly*dPoly) + + if kStep == 1 then + Wkq = Hq^-1 + elseif kStep > 1 then + adCoeff = coeff(aPoly*dPoly);adCoeff = adCoeff(length(adCoeff):-1:1); + adPoly = poly(adCoeff,'q','coeff') + cCoeff = model.c;cCoeff = cCoeff(length(cCoeff):-1:1); + cPoly = poly(cCoeff,'q','coeff') + hBar = clean((ldiv(cPoly,adPoly,kStep))',0.00001) + hBarPoly = poly(hBar,'q','coeff') + Wkq = hBarPoly*Hq^-1 + end + WkqGq = Wkq * Gq + tempWkqGq = coeff(WkqGq.den) + if tempWkqGq(1) <> 1 then + WkqGq.num = WkqGq.num/tempWkqGq(1) + WkqGq.den = WkqGq.den/tempWkqGq(1) + end + Wkq1 = 1-Wkq + tempWkq1 = coeff(Wkq1.den) + if tempWkq1(1) == 1 then + Wkq1.num = Wkq1.num/tempWkq1(1) + Wkq1.den = Wkq1.den/tempWkq1(1) + end + //pause +//------------------------------------------------------------------------------ + // storing the plant data + uCoeff = coeff(WkqGq.num*Wkq1.den) + yCoeff = coeff(WkqGq.den*Wkq1.num) + yCapCoeff = coeff(WkqGq.den*Wkq1.den) + //pause + lengthuCoeff = length(uCoeff) + lengthyCoeff = length(yCoeff) + lengthyCapCoeff = length(yCapCoeff) +//------------------------------------------------------------------------------ + // keeping initial conditions equal to zero + uData = zeros(lengthuCoeff,1) + yData = zeros(lengthyCoeff,1) + yCapData = zeros(lengthyCapCoeff-1,1) + uData = [uData;data(:,2)] + yData = [yData;data(:,1)] + sampleData = size(data,'r') + //pause + // reversing the coefficients + if ~size(uCoeff,'*') then + uCoeff = 0 + else + uCoeff = uCoeff(lengthuCoeff:-1:1) + end + if ~size(yCoeff) then + yCoeff = 0 + else + yCoeff = yCoeff(lengthyCoeff:-1:1) + end + if ~size(yCapCoeff,'*') then + yCapCoeff = 0 + else + yCapCoeff = -yCapCoeff(lengthyCapCoeff:-1:2) + end + //pause + for ii = 1:sampleData+1 + //pause + if ~size(uData(ii:ii+lengthuCoeff-1),'*') then + tempu = 0 + else + tempu = uCoeff*uData(ii:ii+lengthuCoeff-1); + end + if ~size(yData(ii:ii+lengthyCoeff-1),'*') + tempy = 0 + else + tempy = yCoeff*yData(ii:ii+lengthyCoeff-1); + end + if ~size(yCapData(ii:ii+lengthyCapCoeff-2),'*') then + tempyCap = 0 + else + tempyCap = yCapCoeff*yCapData(ii:ii+lengthyCapCoeff-2); + end + yCapData = [yCapData;tempu+tempy+tempyCap]; + end + // pause + extraSample = abs(size(yCapData,'r')-sampleData) + yCapData = yCapData(extraSample+1:$) + timeData = ((modelSampleTime:modelSampleTime:(size(yCapData,'r')*modelSampleTime))'); + if lhs == 1 then + clf() + plot(timeData,yCapData) + axisData = gca() + tempTimeUnit = 'Time('+modelTimeUnit+')' + xtitle('Predicted Response',tempTimeUnit,'y') + xgrid + varargout(1) = 0 + elseif lhs == 2 then + varargout(1) = yCapData + varargout(2) = 0 + elseif lhs == 3 then + varargout(1) = yCapData + varargout(2) = timeData + varargout(3) = 0 + end +endfunction diff --git a/rARX.mat b/rARX.mat new file mode 100644 index 0000000..22171f2 Binary files /dev/null and b/rARX.mat differ diff --git a/rData.mat b/rData.mat new file mode 100644 index 0000000..605635f Binary files /dev/null and b/rData.mat differ diff --git a/rarx.sci b/rarx.sci new file mode 100644 index 0000000..e2c291c --- /dev/null +++ b/rarx.sci @@ -0,0 +1,63 @@ +function varargout = rarx(varargin) + [lhs,rhs] = argn(0) +// + plantData = varargin(1) + orderData = varargin(2) + na = orderData(1);nb = orderData(2) + // arranging na ,nb,nk + if size(orderData,"*") == 2 then + nk = 1 + elseif size(orderData,'*') == 3 then + nk = orderData(3) + end + // storing the lambda value + if rhs == 3 then + lambda = varargin(3) + else + lambda = 0.95 + end + + nb1 = nb + nk - 1 + n = max(na, nb1) + // arranging the plant data + if typeof(plantData) == 'constant' then + Ts = 1;unitData = 'second' + elseif typeof(plantData) == 'iddata' then + Ts = plantData.Ts;unitData = plantData.TimeUnit + plantData = [plantData.OutputData plantData.InputData] + end + N = size(plantData,'r') + uIndex = nk:nb1 + yIndex = [] + if na <> 0 then + yIndex = 1:na + end + df = N - na - nb + Plast = 10^4 * (eye(na+nb,na + nb)) + theta = zeros(N + 1, na + nb) + yHat = plantData(:,1);yData = plantData(:,1) + tempData = zeros(N,na+nb) + for ii = 1:na + tempData(ii+1:ii+N,ii) = -plantData(:,1) + end + // arranging samples of u matrix + for ii = 1:nb + tempData(ii+nk:ii+N+nk-1,ii+na) = plantData(:,2) + end + //tempData = [zeros(1,na+nb);tempData] + tempData = tempData(1:N+1,:) + + for ii = 1:N + temp = tempData(ii,:) + yHat(ii) = temp*theta(ii,:)' + eps_i = yData(ii)-yHat(ii) + kappa_i = Plast * temp'/(lambda + temp * Plast * temp') + theta(ii+1,:) = ((theta(ii,:))' + eps_i * kappa_i)' + Plast = (eye(na + nb,na + nb) - kappa_i * (temp)) * Plast/lambda + + end + theta = theta(1:N,:) + yHat = yHat(1:N) + + varargout(1) = struct('theta',theta,'yHat',yHat) +endfunction diff --git a/sim.sci b/sim.sci new file mode 100644 index 0000000..9cd03f1 --- /dev/null +++ b/sim.sci @@ -0,0 +1,65 @@ +function varargout = sim(varargin) + [lhs,rhs] = argn(0) + //checking the number of inputs + if rhs < 2 || rhs > 3 then + error(msprintf(gettext("%s: Unexpected number of input arguments "),"sim")) + end + modelData = varargin(2) + inputData = varargin(1) + //checking the first input + if typeof(inputData) <> "constant" && typeof(inputData) <> "iddata" then + error(msprintf(gettext("%s: Plant input data must be ""iddata"" type or ""double vector"". "),"sim")) + end + if typeof(inputData) == "iddata" then + inputData = inputData.InputData + end + if ~iscolumn(inputData) then + error(msprintf(gettext("%s: Plant input data must be ""double column vector"". "),"sim")) + end + //checking the plant model type + if typeof(modelData) <> "idpoly" then + error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"sim")) + end + //adding noise + if rhs == 3 then + noiseFlag = varargin(3) + if typeof(noiseFlag) <> 'boolean' then + error(msprintf(gettext("%s: Last input data must be ""boolean"".type "),"sim")) + end + else + noiseFlag = %F + end + //adding noise of zero mean and 1 as standard deviation + if noiseFlag then + numberOfSamples = size(inputData,'r') + R = = grand(numberOfSamples,1,"nor",0,1) + inputData = inputData+R + end + z = poly(0,'z') + aPoly = poly(modelData.a(length(modelData.a):-1:1),'z','coeff') + bPoly = poly(modelData.b,'z','coeff') + fPoly = poly(modelData.f(length(modelData.f):-1:1),'z','coeff') + afPoly = aPoly*fPoly + + bCoeff = modelData.b + extra = 1 + if ~bCoeff(1,1) then + afCoeff = coeff(afPoly) + bLength = length(bCoeff);afLength = length(afCoeff) + if bLength == afLength then + extra = 1 + else + extra = z^-(bLength-afLength) + end + end + + bPoly = poly(modelData.b(length(modelData.b):-1:1),'z','coeff') + sys = syslin('d',bPoly,afPoly)*extra + outputData = (flts(inputData',sys))' + + if typeof(varargin(1)) == "iddata" then + outputData = iddata(outputData,[],varargin(1).Ts) + outputData.TimeUnit = varargin(1).TimeUnit + end + varargout(1) = outputData +endfunction diff --git a/spa.sci b/spa.sci new file mode 100644 index 0000000..4b91cc9 --- /dev/null +++ b/spa.sci @@ -0,0 +1,73 @@ +function varargout = spa(varargin) + [lhs , rhs] = argn(); + if ( rhs < 1 || rhs > 3 ) then + errmsg = msprintf(gettext("%s: Wrong number of input arguments" ),"spa"); + error(errmsg) + elseif typeof(varargin(1)) <> "iddata" then + error(msprintf(gettext("%s:Plant data must be ""iddata"" type.\n"),"spa")) + end + + plantData = varargin(1) + windowSize = %F + inputFreq = %F +//------------------------------------------------------------------------------ +// arranging the plant data + inputData = plantData.InputData; + if ~size(inputData,"*") then + error(msprintf(gettext("%s:Input data must be non-zero vector. \n"),"spa")) + end + outputData = plantData.OutputData + if ~size(outputData,"*") then + error(msprintf(gettext("%s:Output data must be non-zero vector. \n"),"spa")) + end +//------------------------------------------------------------------------------ + N = size(inputData,'r') + nout = size(outputData,'c') + nin = size(inputData,'c') + + if ~windowSize then + windowSize = min(N/10, 30) + end + + if inputFreq then + else + inputFreq = (1:128)/128 * %pi/plantData.Ts + end + + M = windowSize + Ryu = xcov(outputData,inputData,M,'biased') + Ruu = xcov(inputData,inputData,M,'biased') + Ryy = xcov(outputData,outputData,M,'biased') + G = [];spect = []; + for ii = 1:nout + phiY = spaCalculation(inputFreq,Ryy,M) + temp = phiY + for jj = 1:nin + phiYU = spaCalculation(inputFreq,Ryu,M)//sapply(freq, cov2spec, Ryu[i, j, ], M) + phiU = spaCalculation(inputFreq,Ruu,M) + G = [G phiYU./phiU] + //pause + temp = temp - phiYU .* conj(phiYU)./phiU + end + spect = [spect temp] + end +// disp(size(spect)) +// disp(spect) + frdData = frd(G',inputFreq',plantData.Ts,spect') + varargout(1) = frdData +endfunction + +function varargout = spaCalculation(varargin) + inputFreq = varargin(1) + xcovData = varargin(2) + Mnumber = varargin(3) + temp = [] + win_l=window('hn',2*Mnumber+1) + for ii = 1:size(inputFreq,'c') + seq1 = exp(-(%i) * (-Mnumber:Mnumber) * inputFreq(ii)) + data = (seq1.*win_l) + data2 = sum(data.*xcovData') + temp = [temp data2] + end + varargout(1) = temp +endfunction diff --git a/stem.sci b/stem.sci new file mode 100644 index 0000000..2bf5b0a --- /dev/null +++ b/stem.sci @@ -0,0 +1,19 @@ +// Stem plot +// Updated (1-12-06) +function stem(x,y,xy) +if argn(2) == 2 +xy = 5; +end; +set("figure_style","new") +plot2d3(x,y,style=2); +//plot2d4(x,y,style=-9); // default mark foreground colour: black +// Can be used instead of xpoly +// But default values cannot be changed +xpoly(x,y,"marks"); +p=get("hdl"); +p.mark_size_unit = 'point'; +p.mark_size = xy; +p.mark_style = 9; +p.mark_foreground = 2; +p.mark_background = 2; +endfunction; diff --git a/stepest.sci b/stepest.sci new file mode 100644 index 0000000..6310971 --- /dev/null +++ b/stepest.sci @@ -0,0 +1,71 @@ +function varargout = stepest(varargin) + [lhs,rhs] = argn(0) + //checking the number of inputs + if rhs > 2 then + error(msprintf(gettext("%s: Unexpected number of input arguments "),"stepest")) + end + modelData = varargin(1) + if typeof(modelData) <> "idpoly" then + error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"stepest")) + end + //adding noise + if rhs == 2 then + noiseFlag = varargin(2) + if typeof(noiseFlag) <> 'boolean' then + error(msprintf(gettext("%s: Last input data must be ""boolean"".type "),"stepest")) + end + else + noiseFlag = %F + end + z = poly(0,'z') + aPoly = poly(modelData.a(length(modelData.a):-1:1),'z','coeff') + bPoly = poly(modelData.b,'z','coeff') + fPoly = poly(modelData.f(length(modelData.f):-1:1),'z','coeff') + afPoly = aPoly*fPoly + + bCoeff = modelData.b + extra = 1 + if ~bCoeff(1,1) then + afCoeff = coeff(afPoly) + bLength = length(bCoeff);afLength = length(afCoeff) + if bLength == afLength then + extra = 1 + else + extra = z^-(bLength-afLength) + end + end + + bPoly = poly(modelData.b(length(modelData.b):-1:1),'z','coeff') + sys = syslin('d',bPoly,afPoly)*extra + if size(find(gsort(abs(roots(sys.den)))>1),'*') then + tempRoot = clean(roots(sys.den)) + [sorted index] = gsort(abs(real(tempRoot))) + tempRoot = tempRoot(index) + n = 26/log10(abs(real(tempRoot(1)))) + + else + finalValue = horner(sys,1)*1 + n = 10 + tempValue = sum(ldiv(sys.num,sys.den,n)) + if finalValue >= 0 then + while finalValue*0.99 >= tempValue || n >=1000 + n = n+10 + tempValue = sum(ldiv(sys.num,sys.den,n)) + end + elseif finalValue < 0 then + while finalValue*0.99 =1000 + n = n+10 + tempValue = sum(ldiv(sys.num,sys.den,n)) + end + end + end + uData = ones(1,n) + yData = flts(uData,sys) + timeData = (0:(n-1))*modelData.Ts + if noiseFlag then + varargout(1) = yData' + else + plot2d2(timeData,yData,2) + varargout(1) = [] + end +endfunction -- cgit