%% SMI
clc;
clear;
close all;
warning off;
%% system parameters
fs = 10e3;
sampleNumPerCycle = 69999;
TestNumber = 1;
dim = 200;% 维数
%
dimenSet = 110;
skipSample = 0*fs;% skip abnormal sample
fileName = './点二/19.3.7 7s_3-10_1_0/3周期输入输出.txt';
testData = load(fileName);
% 去掉第一行数据
testData(1,:) = [];
% 得到三组数据，分别保存
for kk = 1:TestNumber
% idx
startIdx = (kk-1)*sampleNumPerCycle+1 + skipSample;
stopIdx = (kk)*sampleNumPerCycle;
U{kk} = testData(startIdx:stopIdx,1);
% pre-processing
tmp = testData(startIdx:stopIdx,2);
Y{kk} = smooth(tmp,1);
% Y{kk} = U{kk}+randn(size(U{kk}));
end
% SMI
for kk = 1:TestNumber
Uk = U{kk};
Yk = Y{kk};
% SMI
[subSpace,ssfun,~] = moesp(Yk,Uk,dim);
% Therefore, the state space model is obtained by calling ssfun:
[A{kk},B{kk},C{kk},D{kk}] = ssfun(dimenSet);
%
end
% average
A_m = zeros(size(A{1}));
B_m = zeros(size(B{1}));
C_m = zeros(size(C{1}));
D_m = zeros(size(D{1}));
for kk = 1:TestNumber
A_m = A_m + A{kk};
B_m = B_m + B{kk};
C_m = C_m + C{kk};
D_m = D_m + D{kk};
end
A_m = A_m/TestNumber;
B_m = B_m/TestNumber;
C_m = C_m/TestNumber;
D_m = D_m/TestNumber;
% 使用该模型进行重构
sys = ss(A_m,B_m,C_m,D_m);
% x0 = zeros(dimenSet,1);
for kk = 1:TestNumber
% A_m = A{kk};
% B_m = B{kk};
% C_m = C{kk};
% D_m = D{kk};
Uk = U{kk};
Yk = Y{kk};
x0 = randn(dimenSet,1);
%
x_k = x0;
Errmat = zeros(dimenSet+1);
bufSize = 50;
powerMeritBuf = zeros(bufSize,1);
estMeritBuf = zeros(bufSize,1);
for sampleIdx = 1:length(Uk)
x_k_1 = A_m*x_k + B_m*Uk(sampleIdx);
y_k(sampleIdx) = C_m*x_k + D_m*Uk(sampleIdx);
% scale
powerMeritBuf(1:end-1,1) = powerMeritBuf(2:end,1);
powerMeritBuf(end,1) = Yk(sampleIdx);
scale = mean(abs(powerMeritBuf).^2);
% estimation merit
estMeritBuf(1:end-1,1) = estMeritBuf(2:end,1);
estMeritBuf(end,1) = y_k(sampleIdx);
weight = mean(abs(estMeritBuf).^2);
y_k(sampleIdx) = y_k(sampleIdx)*sqrt(scale/weight);
% update state
x_k = x_k_1;
% error vector
est_r = [x_k_1;
y_k(sampleIdx)];
input_r = [x_k;
Uk(sampleIdx)];
err = est_r - input_r;
%
Errmat = Errmat + err*err';
end
% noise detection
Errmat = Errmat/length(Uk);
Q = Errmat(1:dimenSet,1:dimenSet);
R = Errmat(1+dimenSet:end,1+dimenSet:end);
S = Errmat(1:dimenSet,1+dimenSet:end);
%
delay = 2*dim;
figure(kk);
%
tVec = [1:sampleNumPerCycle-skipSample]/fs;
plot(tVec,Yk,'k-');
xlabel('Time');
hold on;
plot(tVec(1:end-delay),y_k(delay+1:end),'r-');
legend('SMI拟合模型');
% error and relative error
errorVec = y_k(delay+1:end).' - Yk(1:end-delay);
%
relativeErr = errorVec./Yk(1:end-delay);
%
meanErr = mean(errorVec)
meanErr_relative = mean(relativeErr)
end
%% output A B C D
dimToDisplay = 3;
A_m(1:dimToDisplay,1:dimToDisplay)
B_m(1:dimToDisplay,:)
C_m(1,1:dimToDisplay)
D_m
%% save result
fp = fopen('./点二/19.3.7 7s_3-10_1_0/SMI_预测.txt','wb');
saveData = [tVec(1:end-delay)' y_k(delay+1:end)'];
for index = 1:size(saveData,1)
% fprintf(fp,'%f\n',sig(index));
fprintf(fp,'%f\n',saveData(index,2));
end
fclose(fp);
