利用协方差驱动的稳定图绘制教程

%稳定图
clc;clear
[FileName,PathName]=uigetfile('*.txt','选择文件');  %显示选择时程文件的对话框
File=strcat(PathName,FileName);  %把内容横向连接成单个字符串  将括号内路径和文件名连成一条字符串/多个字符串连接成单个字符串
fip=fopen(File,'r');       %打开文件，r:读出   fip=+N表示文件打开成功，文件代号是N
TimeVy=load('-ascii',File);  %  load函数默认是读取.mat格式（二进制文件）用的，如果你是文本文件，请加上'-ascii'参数，用load的时候第一行注释要删除，直接从数据开始。 本句是载入时程信息
TimeVy=TimeVy(1:5000,1:5);  %得到添加白噪声的时间序列的行、列值
[TimePointNum,MesNodeTotalNum]=size(TimeVy);%[r,c]=size（A）将矩阵A的行数返回到第一个输出变量r，将矩阵的列数返回到第二个输出变量c  
MesNodeTotalNum=MesNodeTotalNum-1;      %将第一列设置为时间序列    （ 数据中第一列是时间序列 ）  测点数据就是除去第一列其他列        
TimePointNum=TimePointNum;     %重新赋值一下时间点，   如100HZ采样60S，一共有6000个时间点，所以有6000行

prompt={'TOEPLITZ矩阵块数M(>=采样频率/(2*结构基频)):','系统阶数N(<i*实际测点数):'};
dlg_title='Input Parameter of Calculate';
num_lines=1;
def={'50','2'};                                 %系统阶数N=2*n
answer=inputdlg(prompt,dlg_title,num_lines,def);
M=str2num(answer{1,1});
j=TimePointNum-2*M;
Hankel=zeros((2*M+1)*MesNodeTotalNum,j);
for ti=1:j
      m=0;
      for tii=1:(2*M+1)
           km=ti+tii-1;
            m=m+1;
             beginy=(m-1)*MesNodeTotalNum+1;
               endy=m*MesNodeTotalNum;
          Hankel(beginy:endy,ti)=(TimeVy(km,2:(MesNodeTotalNum+1)))';
      end
end