I/O描述到状态空间两种方法（附matlab代码）

作业要求：根据IO描述到状态空间的两种不同算法，编程实现给出能控能观标准型。输入为传递函数系数（a,b）输出为能控能观标准型（A B C D）。
要求任意阶次均可检验。代码如下：

方法一：

clc
clear
b=input('请输入系数矩阵b:');
a=input('请输入系数矩阵a:');

[mb,nb]=size(b);[ma,na]=size(a);
if nb==na+1
    A=zeros(na);B=zeros(na,1);C=zeros(1,na);D=b(1,1);e=eye(na-1);
    A(1:na-1,2:na)=e;
    A(na,:)=-rot90(rot90(a));
    B(na,1)=1;
    for i=1:na
        C(1,i)=b(1,na+2-i)-b(1,1)*a(1,na-i+1);
    end
end
if nb<na+1
     A=zeros(na);B=zeros(na,1);C=zeros(1,na);e=eye(na-1);
    A(1:na-1,2:na)=e;
    C(1,1:nb)=rot90(rot90(b));
    A(na,:)=-rot90(rot90(a));
    B(na,1)=1;
    D=[];
end
A
B
C
D

方法二：

clc
clear
b=input('请输入系数矩阵b:');
a=input('请输入系数矩阵a:');
len_b=length(b);  %b=[bm,bm-1,...,b1,b0]
m=len_b-1;
len_a=length(a);  %a=[an-1,...,a1,a0]
n=len_a;
A1=zeros(n-1,1);
A2=eye(n-1);
A=[A1 A2];
fb=fliplr(a);   %对b向量做水平翻转
A(n,:)=-fb(:)   %输出A
if m~=0
    if m<n
        x=n-m;X=zeros(1,x);b=[X b];
    end
     if m>=n
       b=b(m-n+1:m+1);
     end
   ee=eye(n+1);
   aa=a';
   for i=1:n
       ee(i,i+1:n+1)=aa(1:n-i+1);
   end
   V=b';
   T=ee';
   B=inv(T)*V; B=B(2:n+1);%输出B
   MN=zeros(1,n-1);
    C=[1 MN] ; %输出C
    D=[b(1)];   %输出D
end
if  m==0
    k=zeros(1,n-1);
    B=[k,b(len_b)];
    B=B';
    p=zeros(1,n-1);
    C=[1,p];
    D=[];
end
A
B
C
D