使用奇异值分解绘制Hankel和Toeplitz矩阵

第二步：奇异值分解
运行后要等待几分钟，待程序运行完后，弹出奇异值分解图；

% Toeplitz矩阵(diagonal-constant matrix)，指矩阵中每条自左上至右下的斜线上的元素相同
%Hankel矩阵：是指每一条副对角线上的元素都相等的方阵。
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的时候第一行注释要删除，直接从数据开始。 本句是载入时程信息
[TimePointNum,MesNodeTotalNum]=size(TimeVy); %[r,c]=size（A）%将矩阵A的行数返回到第一个输出变量r，将矩阵的列数返回到第二个输出变量c  
MesNodeTotalNum=MesNodeTotalNum-1;         %将第一列设置为时间序列    （ 数据中第一列是时间序列 ）  测点数据就是除去第一列其他列       
TimePointNum=TimePointNum;     %重新赋值一下时间点，   如100HZ采样60S，一共有6000个时间点，所以有6000行
%{}大括号，用于cell型的数组,元胞数组是MATLAB的一种特殊数据类型，可以将元胞数组看做一种无所不包的通用矩阵，或者叫做广义矩阵。
prompt={'TOEPLITZ矩阵块数M(>=采样频率/(2*结构基频)):'};  %j>20M   所谓结构基频就是结构本身最小的那个固有频率      建立一个提示字符串  
dlg_title='Input Parameter of Calculate';     %赋值 输入计算参数  一会用于对话框的名称
num_lines=1;  %赋值  一会用于对话框显示行数
def={'100'};    %赋值 用于对话框默认值
answer=inputdlg(prompt,dlg_title,num_lines,def);  %输入对话框  inputdlg
M=str2num(answer{1,1});%将字符串类型转换为数值类型
% 难点
j=TimePointNum-2*M;      % Hankel矩阵最后一个矩阵块对应坐标为（2*i，j）,即2*i+j需<=TimePointNum      相当于输入M=i   i需要人为确定
Hankel=zeros((2*M+1)*MesNodeTotalNum,j);     %建立一个零矩阵，后期填数，因为矩阵每一列比前一列递增一个值，所以到最后一个值的时候等于2*i+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;%计算每次赋值在所在行的结束位置