用matlab从滞回曲线提取骨架曲线

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%程序功能：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('jgz3.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; 
%--------------------提取各个滞回环荷载骨架曲线点，储存在矩阵FramePointsPostive和FramePointsNegative中----------------------------- 
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)<ColMaxValueB(2); %比较最大荷载值，较小者删去 
FramePointsPostive(k,:)=[]; 
else 
FramePointsPostive(k-1,:)=[]; 
end; 
end; 
if abs(ColMinValueA(1)-ColMinValueB(1))<4; 
if ColMinValueA(2)<ColMinValueB(2); %比较最小荷载值，较小者删去 
FramePointsNegative(k-1,:)=[]; 
else 
FramePointsNegative(k,:)=[]; 
end; 
end; 
end; 
dlmwrite('FrameCurve.txt',FramePointsPostive,'-append','delimiter','\t','newline','pc','precision','%.2f') 
dlmwrite('FrameCurve.txt',FramePointsNegative,'-append','delimiter','\t','newline','pc','precision','%.2f') 

本贴内容来源于网络，转发供参考！

Author: Taochun Yang

佐雍得尝，与众分享，收获快乐！