// Functions       : FitTF3, FitTF2, FitTF3i, FitTF2i
// Author          : Josef Schuetzenberger - Austria
// License         : Public Domain
// Version         : 1.0
// Date            : 2020-09-09
// Release         : http://fileexchange.scilab.org/toolboxes/585000

// Supported Scilab Version: 6.1

// DESCRIPTION
// Numerically computing (fitting) a Transfer Function to a measured response(y) 
// of an unknown system. The input function(x) could be a step, impulse or whatever.
//
// NOTES
// A step includes a wide range of frequencies and is therefore a good choice.
// You may have to subtract the steady state first and
// shift the data(y) backwards to compensate a time-delay in the response 
 
// PARAMETERS:
//   Inputs:
//     x: Vector of an impulse or step or slope or...
//     y: Vector of a response to the step-function or an impulse or... (Same length than x)
//
//   Output:
//     G: Transfer Function

clear() //Clear Variables
close() //Close graphic figures, progression bars or wait bars

function G=FitTF3(x,y) //third order transfer function (1+s+s^2+s^3)/(1+s+s^2+s^3)
    a=[];x1=0;y1=0;x2=0;y2=0;x3=0;y3=0
    for k=1:length(x), 
        x1=x1+x(k)
        x2=x2+x1
        x3=x3+x2
        y1=y1+y(k)
        y2=y2+y1
        y3=y3+y2
        a($+1,:)=[x3,x2,x1,x(k),y3,y2,y1]
    end
    X= a \ y'
    X=clean(X)
    kp= -X(1)/X(5)
    t1= X(2)/X(1)
    t2= X(3)/X(1)
    t3= X(4)/X(1)    
    t4= X(6)/X(5)
    t5= X(7)/X(5)
    t6= -1/X(5)
    G= syslin('c',kp*(1+t1*%s+t2*%s^2+t3*%s^3),1+t4*%s+t5*%s^2+t6*%s^3);
    //G= syslin('c',(X(1)+X(2)*%s+X(3)*%s^2+X(4)*%s^3),-X(5)-X(6)*%s-X(7)*%s^2+%s^3);
endfunction

function G=FitTF2(x,y) //second order transfer function (1+s+s^2)/(1+s+s^2)
    a=[];x1=0;y1=0;x2=0;y2=0
    for k=1:length(x), 
        x1=x1+x(k)
        x2=x2+x1
        y1=y1+y(k)
        y2=y2+y1
        a($+1,:)=[x2,x1,x(k),y2,y1]
    end
    X= a \ y'
    X=clean(X)
    kp= -X(1)/X(4)
    t1= X(2)/X(1)
    t2= X(3)/X(1)
    t3= X(5)/X(4)
    t4= -1/X(4)
    G= syslin('c',kp*(1+t1*%s+t2*%s^2),1+t3*%s+t4*%s^2);
    //G= syslin('c',(X(1)+X(2)*%s+X(3)*%s^2)*1,-X(4)-X(5)*%s+%s^2);
endfunction

function G=FitTF3I(x,y) //third order transfer function (Integral version) (1+s+s^2+s^3)/(s+s^2+s^3)
    a=[];x1=0;y1=0;x2=0;y2=0;x3=0;y3=0
    for k=1:length(x), 
        x1=x1+x(k)
        x2=x2+x1
        x3=x3+x2
        y1=y1+y(k)
        y2=y2+y1
        y3=y3+y2
        a($+1,:)=[x3,x2,x1,x(k),y2,y1]
    end
    X= a \ y'
    X=clean(X)
    t1= X(2)/X(1)
    t2= X(3)/X(1)
    t3= X(4)/X(1)    
    t4= X(5)/X(1)
    t5= X(6)/X(1)
    t6= 1/X(1)
    G= syslin('c',(1+t1*%s+t2*%s^2+t3*%s^3),-t4*%s-t5*%s^2+t6*%s^3);
   // G= syslin('c',(X(1)+X(2)*%s+X(3)*%s^2+X(4)*%s^3),-X(5)*%s-X(6)*%s^2+%s^3);
endfunction

function G=FitTF2I(x,y) //second order transfer function (Integral version) (1+s+s^2)/(s+s^2)
    a=[];x1=0;y1=0;x2=0;y2=0
    for k=1:length(x), 
        x1=x1+x(k)
        x2=x2+x1
        y1=y1+y(k)
        y2=y2+y1
        a($+1,:)=[x2,x1,x(k),y1]
    end
    X= a \ y'
    X=clean(X)
    t1= X(2)/X(1)
    t2= X(3)/X(1)
    t3= X(4)/X(1)
    t4=1/X(1)
    G= syslin('c',(1+t1*%s+t2*%s^2),-t3*%s+t4*%s^2);
    //G= syslin('c',(X(1)+X(2)*%s+X(3)*%s^2)*1,-X(4)*%s+%s^2);
endfunction

// Demo 

// Recorded data from a csv-file

M=csvRead("FitTF23.csv", [],[], "double")
x=M(:,3)         // a vector of inputs to the plant system  (Inflow to the tank)
o=x(1);x=x-o     // subtract the steady state
y=M(:,1)         // output of the plant system (level of the tank)
y=y-M(:,2)       // subtract the steady state  
G2= FitTF2I(x,y') // fit the second order transfer function (level of a tank)
G3= FitTF3I(x,y') // fit the third order transfer function (level of a tank)
H = time_id(3, x, y') // use the built in SISO least square identification
disp(G2,G3,H)
yr2= flts(x',dscr(G2,1)) // simulate response of 2 order transfer function 
yr3= flts(x',dscr(G3,1)) // simulate response of 3 order transfer function  
yrh= flts(x',H)          // simulate response of 3 order z transfer function  

Mp=[x,y,yr2',yr3',yrh']
plot(Mp)
title("Response of fitted second and third order transfer-functions", "fontsize",3)
xgrid(1)
legend("Input to the plant-system (Inflow)","Response of the plant-system (Level)", prettyprint(G2),prettyprint(G3),prettyprint(H),1,%t)
xlabel("Time (s)")
ylabel("Level/Inflow in (%)")
