滞回环的骨架曲线matlab——可直接复制

本代码源自网络，我仅略微更改后分享。https://www.bilibili.com/video/av540086839/

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%程序功能：1.生成并输出骨架曲线点至FrameCurve.txt，提取后需要对骨架曲线点按照位移值排序；
%默认格式：txt文件为两列，顺序依次为：、
%读取的曲线默认：从正向开始加载，若为负向开始加载，应先进行调整。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
fid=fopen('FrameCurve.txt','wt');    %%创建输出文件FrameCurve.txt，初次写入骨架曲线点[0,0]，应覆盖内容
coordinate=[0.0 0.0];
dlmwrite('FrameCurve.txt',coordinate,'delimiter','\t','newline','pc','precision','%.2f')
sta=fclose(fid);
outdata=load('SJ.txt');                      %%读取文件数据，默认是为,，outdata矩阵,
LineNum=size(outdata,1);
EndCircle=1;
%--------------------找各级循环分割点位置----------------------负->正or负->零
for i=3:LineNum
    if (outdata(i,1)*outdata(i-1,1)<=0)&&(outdata(i-1,1)<0);%--默认初始加载从正半周开始
        EndCirclePoints(EndCircle)=i;                       %分割点的行号EndCirclePoints
        EndCircle=EndCircle+1;                              %中间分割点个数为EndCircle-1
    end;
end;

%--------------------检查最后是否存在不完整滞回环----------------------
if  EndCirclePoints(EndCircle-1)<LineNum
    LoopNum=length(EndCirclePoints)+1;%LoopNum仅为完整循环的个数，应还有最后一个不完整循环
else    
    LoopNum=length(EndCirclePoints);
end;
%-------------------------分割各个滞回环---------------------------
for k=1:LoopNum
    if k==1
        LoopCircles(k)={outdata(1:EndCirclePoints(k),:)};    
    elseif k<LoopNum        
        LoopCircles(k)={outdata(EndCirclePoints(k-1):EndCirclePoints(k),:)};    
    else        
        LoopCircles(k)={outdata(EndCirclePoints(k-1):LineNum,:)};%后续步骤中，最后一段应区分是少于半周、少于整周
    end;
end;
%--------------------提取各个滞回环荷载骨架曲线点---------------
for k=1:LoopNum    
    A=LoopCircles{k};        
    if k<LoopNum
         [ColMaxValue,LineMax]=max(A);%获取各列最大值所在行号LineMax(1*n向量)，ColMaxValue为各列的最大值(1*n向量)，n为A的列数
         [ColMinValue,LineMin]=min(A);%获取各列最小值所在行号LineMin        
    FramePointMax=A(LineMax(2),:);  %第2列为荷载列
    FramePointMin=A(LineMin(2),:);
    FramePointsPostive(k,:)=FramePointMax;      %骨架曲线点储存在矩阵FramePoints中  
    FramePointsNegative(k,:)=FramePointMin;
%-------------------区分最后一段是完整周，少于半周，还是多于半周---------------    
else
    if EndCirclePoints(EndCircle-1)==LineNum    %完整周-------------------------------------------------- 
    [ColMaxValue,LineMax]=max(A);    %获取各列最大值所在行号LineMax(1*n向量)，ColMaxValue为各列的最大值(1*n向量)，n为A的列数            
    [ColMinValue,LineMin]=min(A);          %获取各列最小值所在行号LineMin 
    FramePointMax=A(LineMax(2),:);         %第2列为荷载列 
    FramePointMin=A(LineMin(2),:);
    FramePointsPostive(k,:)=FramePointMax; %骨架曲线点储存在矩阵FramePoints中   
    FramePointsNegative(k,:)=FramePointMin;
    elseif 
        outdata(LineNum,1)>=0               %仅正半周，仅一个最大值
        [ColMaxValue,LineMax]=max(A);          %获取最大值行号LineMax
        FramePointMax=A(LineMax(2),:);           
        FramePointsPostive(k,:)=FramePointMax;        
    else                                %有负半周，有最大值和最小值
        [ColMinValue,LineMin]=min(A);          %获取最小值行号LineMin  
        FramePointMin=A(LineMin(2),:);       
        FramePointsPostive(k,:)=FramePointMax;  
        FramePointsNegative(k,:)=FramePointMin;
     end;   
 end;
 end;
 %------------比较最大位移，删除骨架曲线中同一荷载级的较小骨架点---------------------------
 for k=LoopNum:-1:2    
     A=LoopCircles{k};   
     B=LoopCircles{k-1};   
     [ColMaxValueA,LineMaxA]=max(A);        %获取各列最大值所在行号LineMax(1*n向量)，ColMaxValue为各列的最大值(1*n向量)，n为A的列数    
     [ColMinValueA,LineMinA]=min(A);          %获取各列最小值所在行号LineMin   
     [ColMaxValueB,LineMaxB]=max(B);          %获取各列最大值所在行号LineMax(1*n向量)，ColMaxValue为各列的最大值(1*n向量)，n为A的列数
     [ColMinValueB,LineMinB]=min(B);          %获取各列最小值所在行号LineMin     
     if abs(ColMaxValueA(1)-ColMaxValueB(1))<4;         %判断为同一荷载级，容差取为4mm  
         if ColMaxValueA(2)2);            %比较最大荷载值，较小者删去  
             FramePointsPostive(k,:)=[];       
         else          
            FramePointsPostive(k-1,:)=[];        
         end;   
     end;    
     if abs(ColMinValueA(1)-ColMinValueB(1))<4;     
         if ColMinValueA(2)2);            %比较最小荷载值，较小者删去
             FramePointsNegative(k-1,:)=[];        
         else            
             FramePointsNegative(k,:)=[];       
         end;   
     end;     
end;
C=[FramePointsPostive;FramePointsNegative]
G=sortrows(C,1)dlmwrite('FrameCurve.txt',G)          %骨架曲线数值存储在FrameCurve.txt文档中

d =G(:,1);e =G(:,2);
a=outdata(:,1);c=outdata(:,2);plot(a,c,'r');
hold on;
plot(d,e,'s-k');   %画出滞回曲线与骨架曲线