%----------信號提取-------------%
clear off ;
clc
f1=fopen('shiyan101.DTS','r');

[a1,counta1]=fread(f1,16384,'uint8');%讀取8進制的信號，共16384個字節
fclose(f1);
j1=0;
for i11=1:2:counta1-4
    j1=j1+1;
    aaa1(j1)=a1(i11)+256*a1(i11+1)-2048;
    aaaa1(j1)=aaa1(j1)*5.0/2048; %將8進制的信號轉化為10進制的信號
end

f2=fopen('shiyan102.DTS','r');
[a2,counta2]=fread(f2,16384,'uint8');
fclose(f2);
j2=0;
for i22=1:2:counta2-4
    j2=j2+1;
    aaa2(j2)=a2(i22)+256*a2(i22+1)-2048;
    aaaa2(j2)=aaa2(j2)*5.0/2048; 
end
 
f3=fopen('shiyan103.DTS','r');
[a3,counta3]=fread(f3,16384,'uint8');
fclose(f3);
j3=0;
for i33=1:2:counta3-4
    j3=j3+1;
   aaa3(j3)=a3(i33)+256*a3(i33+1)-2048;
   aaaa3(j3)=aaa3(j3)*5.0/2048; 
    
end

figure(1)
plot(aaaa1)
xlabel('時間/s')
ylabel('電壓/v')
title('shiyan101信號圖')
set(gcf,'color','w')

figure(2)
plot(aaaa2)
xlabel('時間/s')
ylabel('電壓/v')
title('shiyan102的信號圖')
set(gcf,'color','w')

figure(3)
plot(aaaa3)
xlabel('時間/s')
ylabel('電壓/v')
title('shiyan103的信號圖')
set(gcf,'color','w')
save('shiyan101','aaaa1')
save('shiyan102','aaaa2')
save('shiyan103','aaaa3')
%------------FFT變換----------%
n=512; %設定每周期數據采集點數,即采樣頻率
x=aaaa1;
T=length(x)/512;                    %采樣周期數
t=0:2*pi/(n-1):2*T*pi;              %采樣數
N=length(t)

x=aaaa1;
y=fft(x);
m=abs(y);
p=unwrap(angle(y));
f=(0:length(y)-1)*512/length(y);
figure(4)
subplot(211)
plot(f,m)
xlabel('頻率/Hz')
ylabel('幅值/dB')
title('shiyan101的幅頻圖')
grid on
set(gcf,'color','w')
subplot(212)
plot(f,p*180/pi)
xlabel('頻率/Hz')
ylabel('相角/rads')
title('shiyan101的相頻圖')
grid on
set(gcf,'color','w')%shiyan101
tf=(0:1/512:1);
x=aaaa2;
y=fft(x);
m=abs(y);
p=unwrap(angle(y));
f=(0:length(y)-1)*512/length(y);
figure(5)
subplot(211)
plot(f,m)
xlabel('頻率/Hz')
ylabel('幅值/dB')
title('shiyan102的幅頻圖')
grid on
set(gcf,'color','w')
subplot(212)
plot(f,p*180/pi)
xlabel('頻率/Hz')
ylabel('相角/rads')
title('shiyan102的相頻圖')
grid on
set(gcf,'color','w')%shiyan102
tf=(0:1/512:1);
x=aaaa3;
y=fft(x);
m=abs(y);
p=unwrap(angle(y));
f=(0:length(y)-1)*512/length(y);
figure(6)
subplot(211)
plot(f,m)
xlabel('頻率/Hz')
ylabel('幅值/dB')
title('shiyan103的幅頻圖')
grid on
set(gcf,'color','w')
subplot(212)
plot(f,p*180/pi)
xlabel('頻率/Hz')
ylabel('相角/rads')
title('shiyan103的相頻圖')
grid on
set(gcf,'color','w')%shiyan103

%-----------互相關分析-------------%
x=aaaa2;
y=aaaa3;
figure(7)
[c_ww,lags] = xcorr(x,y,'coeff');
[y_max,t_max]=max(c_ww)              %找出c_ww的最大值及對應的t_max

sample_delay=t_max-N;              %計算與中心點相差的樣點數
T_sig=(N-1)/6;                     %求信號一個周期的樣點數,一個周期相對應于2*pi
delay1=pi/6                        %pi/6的弧度
delay2=2*pi*sample_delay/T_sig     %計算與中心點相差的樣點數所對應的弧度值
%%%%上面四步主要是看看能不能由互相關函數求相位差

stem(lags,c_ww)
xlabel('時間/s')
ylabel('互相關函數值Rxy(t)');
title('shiyan1的互相關估計')
grid on
set(gcf,'color','w')
Cxy=corrcoef(x,y)
% ------------模態分析--------------%
Y1=aaaa2;
Y2=aaaa3;
X= aaaa1;
IDataXY=[Y1' Y2'  X'];
figure(8);
plot(IDataXY);
grid on
set(gcf,'color','w')
title('參與辯識的數據');
m1=n4sid(IDataXY,[1:10],'ssp','can');%用n4sid系統辨識函數辨識系統
mss=idss(m1)
A=m1.a;
B=m1.b;
C=m1.c;
D=m1.d;
[num,den]=ss2tf(A,B,C,D,1)%狀態空間-傳遞函數系數
figure(10)
bode(num,den)
grid on
set(gcf,'color','w')