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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
|
function varargout = sim(varargin)
// Simulate response of idpoly type polynomials
//
// Calling Sequence
// yData = sim(uData,sys)
// Parameters
// sys : idpoly type polynomial
// uData : plant input data in iddata class or nx1 matrix
// yData : simulated plant output data nx1 matrix
// Description
// sim function returnes the simulated response of idpoly type dynamic function for the given input data
//
// Examples
// uData = idinput(1024,'PRBS',[0 1/20],[-1 1])
// a = [1 0.2];b = [0 0.2 0.3];
// sys = idpoly(a,b,'Ts',0.1)
// yData = sim(uData,sys)
// Authors
// Ashutosh Kumar Bhargava
[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
|
