1 files changed, 71 insertions, 0 deletions
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 <tempValue || n >=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