SMI.zip_matlab_moesp

%% 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);