added scilabs files

1 files changed, 75 insertions, 0 deletions

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