summaryrefslogtreecommitdiff
path: root/sim.sci
blob: 9cd03f1579dd0489c22f8245241cfaa8d2001524 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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