structidentification.m

系统辨识中

%系统辨识实验（M序列）
%
%M序列 n＝4 Np＝2^4-1=15
clear;
clc;
load('expdata1.txt');
rawData=expdata1;%原始数据
sizeData=size(rawData);%获得数据大小
rawOutData=rawData(:,1);%原始数据输出
rawInData=rawData(:,2);%原始数据输入
%输出数据求和
rawOutAll=0;
for i=1:sizeData(1,1)
    rawOutAll=rawOutAll+rawOutData(i);
end
rawOutAve=rawOutAll/sizeData(1,1);%原始输出平均值
outData=rawOutData-rawOutAve;%扣除直流分量后的输出数据
inData=rawInData-45;%扣除直流分量后的输入数据
%
%不进行滤波处理
%
%%计算输出与输入之间的互相关函数
Rmz=zeros(1,sizeData(1,1));%init
Rmz=xcorr(outData,inData,'biased');   %%计算互相关函数   'biased'   - scales the raw cross-correlation by 1/M.
c=-Rmz(sizeData(1,1)-1); %c为补偿量 P.121
%%%%计算输入信号的自相关函数
Rm=zeros(1,sizeData(1,1));
Rm=xcorr(inData,'biased');   %%计算自相关函数 XCORR(A), when A is a vector, is the auto-correlation sequence.   
%%返回脉冲响应序列
G=zeros(1,sizeData(1,1));%inition g
for i=1:sizeData(1,1)
    G(i)=sizeData(1,1)*(Rm(i)+c)/(sizeData(1,1)+1)/5^2/10; % see P.115 4.5.32
end
t=1:15;
for i1=1:15
    g(i1)=G(i1+15);%从第二个周期开始取数
end
plot(t,g);
% Rmz1=Rmz';
% plot(t,Rmz1);
%H(l,k)=[g(k)     g(k+1) ...
%g(k+l-1);%k以1至dataNumber-2l+2取之间意选(P332).这我取k=50,l=10(保险起见)
%        g(k+1)   g(k+2) ...   g(k+1);
%        ...
%        g(k+l-1) g(k+l-1)...  g(k+2l-2)]
% case little noise
for l=1:15 %sizeData(1,1)
    %令|H(l)|=H(l)
    H(l)=0;
    for k=1:(15-2*l+2)  %(sizeData(1,1)-2*l+2)
        for i=1:l
            for j=1:l
                h(i,j)=G(k+j+i-2);
            end %end j
        end %end i
        %求detH(l,k)
        detH=det(h);
        H(l)=H(l)+detH;
    end %end k
    H(l)=H(l)/(15-2*l+2);   
end %end l
for n=1:8 %(sizeData(1,1)-1)%根据经验n一定小余8
    D(n)=H(n+1)/H(n+2);
end %end n
D %cout D
% rankh=rank(h);
% %rankSysMatrix=rank(h);%get system's real rank ;this is 3
% if det(h)==0
%     n=15;%the rank of the system
%     break
% end
% end
% x=1:sizeData(1,1);
% plot(x,G);