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
|
//Author: Parthasarathi Panda
//parthasarathipanda314@gmail.com
function [sos,g]=ss2sos(A,B,C,D)
//not taking if, order and scale as input since they do not seem useful
if (type(A)~=1 | type(B)~=1 | type(C)~=1 | type(D)~=1) then
error('check input types');
end
if (length(size(A))~=2 | length(size(B))~=2 | length(size(C))~=2 | length(size(D))~=2) then
error('input must be 2d matrices');
end
if (norm(imag(A))~=0 | norm(imag(B))~=0 | norm(imag(C))~=0 | norm(imag(D))~=0) then
error('input must be real matrices');
end
//cross checking dimensions of matrix
if size(D)~=[1,1] then
error('for single input single output, D must be 1 by 1');
end
[n,k]=size(B);
if k~=1 then
error('for single input single output, B must be column matrix');
end
[n,k]=size(C);
if n~=1 then
error('for single input single output, C must be row matrix');
end
if size(A)~=[1,1] then
error('A must be square matrix');
end
//obtaining the transfer function(continuous)
tf=ss2tf(syslin('c',A,B,C,D));
//factorising the numerator and the denominator into second order systems
[zero,gn]=sosbreak(numer(tf));//function is defined in the same folder
[pole,gd]=sosbreak(denom(tf));
//reducing each pair of second order in the necessary form
sos=[];
for i=[1:length(pole)]
den=pole(i);
// if no polynomial is left in the numerator
if i>length(zero) then
//z^(-1) is in the numerator if the denominator in linear
if degree(den)==1 then
b=[0,1,0];
//z^(-2) is in the numerator if the denominator is quadratic
else
b=[0,0,1];
end
//if polynomial is left
else
b=coeff(zero(i));//obtain coefficient of the polynomial
b=[b,zeros(1,degree(den)+1-length(b))];//add a zero in the end if the numerator is linear and denominator is quadratic
b=b($:-1:1);//reverse the coefficients
b=[b,zeros(1,3-length(b))];//add a zero if both numerator and denominator are linear
end
a=coeff(den);//coefficient of the denominator
a=a($:-1:1);//reversing the coefficients
a=[a,zeros(1,3-length(a))];//adding zeros if the denominator is linear
v=[b,a];
sos=[sos;v]//adding the second order sub-system
end
g=gn/gd;//computing the gain
endfunction
|
