网上找的代码有些问题，自己根据论文又改写了下。感觉效果还行。
测试案例：
生成了100个-5到5的随机点，连成载荷谱
输入的载荷谱：
提取出来的每个载荷循环
结果：
clear;clc%%C=xlsread(‘load_F’);%%%%%%%%%在此修改加载的文件名，数据格式一直才可正确运算%%%%C=[0;15;4;11;1;19;-10;3;-16;8;-14;-2;-9;16;0];%%例子%C=[5;-1;3;-4;4;-2;1;-3;0;-2;5]C=rand(100,1)*10-5;for i=1:length(C)Cx(i)=i;end;B=C;A=C;n=length(A);%% 步骤1%%%对载荷时间历程进行处理使它只包含峰谷峰谷交替出现for i=2:1:n-1if A(i-1)<=A(i)&&A(i)<=A(i+1)B(i)=NaN;elseif A(i-1)>=A(i)&&A(i)>=A(i+1)B(i)=NaN;endendB(isnan(B))=[];%%%%%%%%%%%%%%%%%%%%%%%% B为改造后载荷时间历程 n为B中波峰波谷的个数%% 步骤2，一次雨流计数%%% 幅值F 均值J F=[];J=[];f=0;ff=0;k=1;num=0;%%f为一次雨流计数标志位，ff为二次标志位while length(B)>3n=length(B);if n<=3||ff==2breakelseif n>3c=0;%%前两个点的幅值差if(f==1)%%步骤3，一次雨流计数结束，拼接载荷谱B=div(B,2);%%在最小值处分割n=length(B);ff=2;endfor j=1:n-2if j+2>nbreak;ends1=abs(B(j+1)-B(j));s2=abs(B(j+1)-B(j+2));e3=(B(j)+B(j+1))/2;if s1<=s2&&s1<=cF=[F;s1];J=[J;e3];D(k)=B(j);k=k+1;%%存取提取载荷点B(j)=[];%%剔除jD(k)=B(j);k=k+1;B(j)=[];%%剔除j+1D(k)=D(k-2);k=k+1;D(k)=NaN;k=k+1;num=num+1;c=s1;f=0;ff=0;n=length(B);elsef=1;c=s1;endendendend%%步骤4，计数剩下最后一个循环F=[F;(max(B)-min(B))];J=[J;(max(B)+min(B))/2];D(k)=max(B);D(k+1)=min(B);D(k+2)=max(B);num=num+1;for i=1:length(D)Dx(i)=i;end;%% 步骤5 %%%画图像 三维hist三维图像X=[J,F];figure;plot(Dx,D);title(‘提取载荷循环’);ylabel(‘幅值’);hist3(X,[30 30]);xlabel(‘均值’);title(‘雨流计数法’);ylabel(‘幅值’);zlabel(‘循环次数’);histogram(F);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%div函数function [out] = div(in,op)% 雨流计数法分割函数，输入必须为列向量% op=1从最大值分割，op=2从最小值分割if(op==1)[a,b]=max(in);else if(op==2)[a,b]=min(in);endendn=length(in);B1=in(b:n);B2=in(1:b);B=[B1;B2];A=B;for i=2:1:n-1if A(i-1)<=A(i)&&A(i)<=A(i+1)A(i)=NaN;elseif A(i-1)>=A(i)&&A(i)>=A(i+1)A(i)=NaN;endendA(isnan(A))=[];out=A;end