clc;
clear all;

format long;

fid1= fopen('E:\孫國偉\work\data61.txt','r');
data1=fscanf(fid1,'%f',[12,30]);

fid2= fopen('E:\孫國偉\work\data62.txt','r');
data2=fscanf(fid2,'%f',[12,30]);

fid3= fopen('E:\孫國偉\work\data63.txt','r');
data3=fscanf(fid3,'%f',[12,30]);

fid4= fopen('E:\孫國偉\work\data64.txt','r');
data4=fscanf(fid4,'%f',[12,30]);

fid5= fopen('E:\孫國偉\work\data65.txt','r');
data5=fscanf(fid5,'%f',[12,30]);

fid6= fopen('E:\孫國偉\work\data66.txt','r');
data6=fscanf(fid6,'%f',[12,30]);

fid7= fopen('E:\孫國偉\work\data67.txt','r');
data7=fscanf(fid7,'%f',[12,30]);

fid8= fopen('E:\孫國偉\work\data68.txt','r');
data8=fscanf(fid8,'%f',[12,30]);

fid9= fopen('E:\孫國偉\work\data69.txt','r');
data9=fscanf(fid9,'%f',[12,30]);

fid10= fopen('E:\孫國偉\work\data610.txt','r');
data10=fscanf(fid10,'%f',[12,30]);

fid11= fopen('E:\孫國偉\work\data611.txt','r');
data11=fscanf(fid11,'%f',[12,30]);

fid12= fopen('E:\孫國偉\work\data612.txt','r');
data12=fscanf(fid12,'%f',[12,30]);

%系統的觀測矩陣C陣

A=zeros(12,12,30);

A(1,:,:)=data1;
A(2,:,:)=data2;
A(3,:,:)=data3;
A(4,:,:)=data4;
A(5,:,:)=data5;
A(6,:,:)=data6;
A(7,:,:)=data7;
A(8,:,:)=data8;
A(9,:,:)=data9;
A(10,:,:)=data10;
A(11,:,:)=data11;
A(12,:,:)=data12;

%系統的觀測矩陣C陣

C=zeros(3,12);

C(1,4)=1;
C(2,5)=1;
C(3,3)=1;

B=eye(12);

a=zeros(12,12,30);

%離散化

for i=1:30
    
    [a(:,:,i),b]=c2d(A(:,:,i),B,2);

end

Q=[C;C*a(:,:,1);C*a(:,:,1)^2;C*a(:,:,1)^3;C*a(:,:,1)^4;C*a(:,:,1)^5;C*a(:,:,1)^6;C*a(:,:,1)^7;C*a(:,:,1)^8;C*a(:,:,1)^9;C*a(:,:,1)^10;C*a(:,:,1)^11];

m=rank(Q);

fprintf('%d\n\n',m);

y=[0.2 0.2 0.3*60]';

Y=zeros(36,30);

for i=1:30
    
Y(:,i)=[y;y;y;y;y;y;y;y;y;y;y;y];

end

Y1=Y(:,1);

for i=2:15
    
    Q=[Q;C;C*a(:,:,i);C*a(:,:,i)^2;C*a(:,:,i)^3;C*a(:,:,i)^4;C*a(:,:,i)^5;C*a(:,:,i)^6;C*a(:,:,i)^7;C*a(:,:,i)^8;C*a(:,:,i)^9;C*a(:,:,i)^10;C*a(:,:,i)^11];
    
    Y1=[Y1;Y(:,i)];
    
end

r1=rank(Q([1:6,7,9,4+36*13:6+36*13,9+36*13],:));
% r1=rank(Q([1:6,8:9],:));

fprintf('%d\n\n',r1);
% 
[u,s,v]=svd(Q);

r=rank(Q);

fprintf('%d\n\n',r);


for i=1:12
    
    fprintf('%20.14f\n',s(i,i));
    
end

% x01=u(:,1)'*Y1*v(:,1)/s(1,1);
% x02=u(:,2)'*Y1*v(:,2)/s(2,2);
% x03=u(:,3)'*Y1*v(:,3)/s(3,3);
% x04=u(:,4)'*Y1*v(:,4)/s(4,4);
% x05=u(:,5)'*Y1*v(:,5)/s(5,5);
% x06=u(:,6)'*Y1*v(:,6)/s(6,6);
% x07=u(:,7)'*Y1*v(:,7)/s(7,7);
% x08=u(:,8)'*Y1*v(:,8)/s(8,8);
% x09=u(:,9)'*Y1*v(:,9)/s(9,9);
% x10=u(:,10)'*Y1*v(:,10)/s(10,10);
% x11=u(:,11)'*Y1*v(:,11)/s(11,11);
% x12=u(:,12)'*Y1*v(:,12)/s(12,12);
% 
% %畫直方圖
% 
% figure(1)
% 
% subplot(3,3,1)
% bar(x01);
% legend('奇異值＝5574.39200840085230');
% subplot(3,3,2)
% bar(x02);
% legend('奇異值＝5574.32856625842490');
% subplot(3,3,3)
% bar(x03);
% legend('奇異值＝155.02313111350495');
% subplot(3,3,4)
% bar(x04);
% legend('奇異值＝155.02227952739554');
% subplot(3,3,5)
% bar(x05);
% legend('奇異值＝63.75926327704499');
% subplot(3,3,6)
% bar(x06);
% legend('奇異值＝2.59883066938513');
% subplot(3,3,7)
% bar(x07);
% legend('奇異值＝1.91237996824732');
% subplot(3,3,8)
% bar(x08);
% legend('奇異值＝1.91237124393970');
% subplot(3,3,9)
% bar(x09);
% legend('奇異值＝0.00000000000002');
% 
% figure(2)
% 
% subplot(3,3,1)
% bar(x10);
% legend('奇異值＝0.00000000000001');
% subplot(3,3,2)
% bar(x11);
% legend('奇異值＝0.00000000000000');
% subplot(3,3,3)
% bar(x12);
% legend('奇異值＝0.00000000000000');
% 
% 
% 
% 
% 
% 
% 
% 
% 
%