clc;
clear all;

format long;

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

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

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

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

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

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

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

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

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

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

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

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

fid13= fopen('E:\孫國偉\work\data313.txt','r');
data13=fscanf(fid13,'%f',[18,30]);

fid14= fopen('E:\孫國偉\work\data314.txt','r');
data14=fscanf(fid14,'%f',[18,30]);

fid15= fopen('E:\孫國偉\work\data315.txt','r');
data15=fscanf(fid15,'%f',[18,30]);

fid16= fopen('E:\孫國偉\work\data316.txt','r');
data16=fscanf(fid16,'%f',[18,30]);

fid17= fopen('E:\孫國偉\work\data317.txt','r');
data17=fscanf(fid17,'%f',[18,30]);

fid18= fopen('E:\孫國偉\work\data318.txt','r');
data18=fscanf(fid18,'%f',[18,30]);


%系統的觀測矩陣C陣

A=zeros(18,18,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;
A(13,:,:)=data13;
A(14,:,:)=data14;
A(15,:,:)=data15;
A(16,:,:)=data16;
A(17,:,:)=data17;
A(18,:,:)=data18;

%系統的觀測矩陣C陣

C=zeros(5,18);

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

B=eye(18);

a=zeros(18,18,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;C*a(:,:,1)^12;C*a(:,:,1)^13;C*a(:,:,1)^14;C*a(:,:,1)^15;C*a(:,:,1)^16;C*a(:,:,1)^17];

y=[0.2 0.2 0.1*60 0.1*60 0.3*60]';

Y=zeros(90,30);

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

end

Y1=Y(:,1);

for i=2:30
    
    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;C*a(:,:,i)^12;C*a(:,:,i)^13;C*a(:,:,i)^14;C*a(:,:,i)^15;C*a(:,:,i)^16;C*a(:,:,i)^17];
    
    Y1=[Y1;Y(:,i)];
    
end

% r1=rank(Q([1:10,96+90*3:97+90*3],:));
% 
% fprintf('%d\n\n',r1);

[u,s,v]=svd(Q);

r=rank(Q);

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


for i=1:18
    
    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);
x13=u(:,13)'*Y1*v(:,13)/s(13,13);
x14=u(:,14)'*Y1*v(:,14)/s(14,14);
x15=u(:,15)'*Y1*v(:,15)/s(15,15);
x16=u(:,16)'*Y1*v(:,16)/s(16,16);
x17=u(:,17)'*Y1*v(:,17)/s(17,17);
x18=u(:,18)'*Y1*v(:,18)/s(18,18);


%畫直方圖

figure(1)

subplot(3,3,1)
bar(x01);
subplot(3,3,2)
bar(x02);
subplot(3,3,3)
bar(x03);
subplot(3,3,4)
bar(x04);
subplot(3,3,5)
bar(x05);
subplot(3,3,6)
bar(x06);
subplot(3,3,7)
bar(x07);
subplot(3,3,8)
bar(x08);
subplot(3,3,9)
bar(x09);

figure(2)

subplot(3,3,1)
bar(x10);
subplot(3,3,2)
bar(x11);
subplot(3,3,3)
bar(x12);
subplot(3,3,4)
bar(x13);
subplot(3,3,5)
bar(x14);
subplot(3,3,6)
bar(x15);
subplot(3,3,7)
bar(x16);
subplot(3,3,8)
bar(x17);
subplot(3,3,9)
bar(x18);