// Functions       : FitPT1, FitPT2, FitPT3 
// Author          : Josef Schuetzenberger - Austria
// License         : Public Domain
// Version         : 1.0
// Date            : 2020-08-10
// Release         : http://fileexchange.scilab.org/toolboxes

// 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)
//
//   Outputs:
//     G: Transfer Function
//     Kp: Proportional gain
//     tau, t1, t2, t3: Time constants

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

function [G,kp,tau]=FitPT1(x,y)// Fits a PT1-Element from a response(y) to an impulse, step...
    A=[];x1=0;y1=0
    for k=1:length(x), 
        x1=x1+x(k)
        y1=y1+y(k)
        A($+1,:)=[x1,y1]
    end
    x= A \ y'
    tau= -1/x(2)
    kp= -x(1)/x(2)
    G= syslin('c',kp,1+%s*tau)
endfunction

function [G,kp,t1,t2]=FitPT2(x,y)// Fits a PT2-Element from a response(y) to an impulse, step...
    A=[];x1=0;y1=0;x2=0;y2=0
    for k=1:length(x), 
        y1=y1+y(k)        
        x1=x1+x(k)
        y2=y2+y1
        x2=x2+x1
        A($+1,:)=[x2,y2,y1]
    end
    x= A \ y'
    t2= -1/x(2)
    t1= x(3)/x(2)
    kp= -x(1)/x(2)
    G= syslin('c',kp,1+t1*%s+t2*%s^2)   
endfunction

function [G,kp,t1,t2,t3]=FitPT3(x,y) // Fits a PT3-Element from a response(y) to an impulse, step...
    A=[];x1=0;y1=0;x2=0;y2=0;x3=0;y3=0
    for k=1:length(x), 
        y1=y1+y(k)        
        x1=x1+x(k)
        y2=y2+y1
        x2=x2+x1
        y3=y3+y2
        x3=x3+x2
        A($+1,:)=[x3,y3,y2,y1]
    end    
    x= A \ y'  
    t1= x(3)/x(2)
    t2= x(4)/x(2)
    t3= -1/x(2)
    kp= -x(1)/x(2)
    G= syslin('c',kp,1+t1*%s+t2*%s^2+t3*%s^3)   
    //G = syslin('c',x(1),-x(2)-x(3)*%s-x(4)*%s^2+%s^3)    
endfunction

//Demo 

//Data generation to simulate a noisy system
len=990
x=ones(1,len)                            //Create step 
y=(1-%e^-(10*linspace(0,1,len)))*2^(1/3) //Create test data 
y=y.*y.*y                    // Distort test data
y=y+0.01*rand(x,'normal')    // Add randomness
y=y+grand(x, "nor",0,0.04)   // Add more randomness

//Fitting
[G1,kp,tau]=FitPT1(x,y)     //Fit a PT1 transfer-function 
PT1fit=(1-%e^(-len*(linspace(0,1,len))/tau))*kp //Calc step response
[G2,kp,t1,t2]=FitPT2(x,y)   //Fit a PT2 transfer-function (denominator will be a second degree polynomial)
[G3,kp,t1,t2,t3]=FitPT3(x,y)//Fit a PT3 transfer-function (denominator will be a third degree polynomial)
disp(G1,G2,G3)

//Plotting results
t = [0:1:len] //Time range
figure(1)
subplot(2,1,1)
title("step response of fitted PT1,PT2,PT3 transfer-functions", "fontsize",3)
a=gcf()
a.background=8
plot(y,'r+-') //original response to step
plot(PT1fit,'blue')
plot(csim('step',t,G2),'black')
plot(csim('step',t,G3),'green')
legend("Original Data to fit", prettyprint(G1),prettyprint(G2),prettyprint(G3),4,%t)
xlabel("Discrete Time Index n−−−−>")
ylabel("Step Response −−−−−>")

frq=0:0.002:0.15     // Frequency range
subplot(2,1,2)
title("frequency response of fitted PT1,PT2,PT3 transer-functions", "fontsize",3)
xgrid(1)
rep=repfreq(G1,frq)  // Complex frequency response
[db,phi]=dbphi(rep)           
plot("ln",frq,db,'blue')  // Plot magnitude of frequency response in db 
rep=repfreq(G2,frq)
[db,phi]=dbphi(rep)
plot("ln",frq,db,'bla')
rep=repfreq(G3,frq)
[db,phi]=dbphi(rep)
plot("ln",frq,db,'g')
xlabel(' Frequency −−−−>')
ylabel( ' Amplitude in db --->')
