MATLAB雨流计数法程序
一、仿真随机过程
sw=2;         
A=sqrt(2*sw*detaw)
（1）中心频率为10pi，detaw=0.01，带宽为2pi的低频窄带随机过程
for n=9:0.01:11;     
    xn=A*sin(n*pi*t+rand(1,1)*2*pi);   
x=x+xn;
end
（2）中心频率为100pi，detaw=0.01，带宽为2pi的高频窄带随机过程
for n=99:0.01:101;   
    yn=A*sin(n*pi*t+rand(1,1)*2*pi);
    y=y+yn;
end
（3）低频与高频组合成的宽带随机过程
z=x+y
问题一：关于带宽。带宽为2pi似乎太宽了，若带宽改为0.2pi，detaw依然不变，得到图如下。
二、用雨流计数法截取循环
（1）得到随机过程x的极点值（储存在矩阵S中）。
n=length(x);
s=x(1);
for i=2:n-1;
    if (x(i)>x(i-1))&&(x(i)>x(i+1))||((x(i)<x(i-1))&&(x(i)<x(i+1)));
        s=[s,x(i)];
    end
end
S=[s,x(n)];
（2）对极值点间变程进行比较分析，截取循环
while i+3<length(S) 
    i=i+1;
    detaS1=abs(S(i+1)-S(i));  %获得四个极值点间的三段变程
    detaS2=abs(S(i+2)-S(i+1));
    detaS3=abs(S(i+3)-S(i+2));
    if (detaS2<=detaS1)&&(detaS2<=detaS3) %将变程满足截取条件的循环截取出
        Ba=[Ba,(S(i+2)-S(i+1))/2];    %截取出循环的幅值
        Bm=[Bm,(S(i+2)+S(i+1))/2];    %截取出循环的均值
        Q=[Q;S(i+1),S(i+2)];          %构成截取出循环的极值点存于Q中，即原随机过程抛弃的点
        S(i+1)=[];
        S(i+2)=[];
        i=i-1;
    end
end
Sr=S;  %剩余的极值点存于Sr中
（3）分别对低频窄带、高频窄带、宽带进行雨流处理
低频窄带
高频窄带
宽带
问题二：由于观察到用雨流法处理这三种随机过程的效果差别比较大，就编程统计计算了雨流处理这三种随机过程抛弃的极值点数（即截取的循环）与原载荷极值点数之比rate。因为我理解为rate越大，雨流处理效果越好（不知道对不对？）得到的统计结果如下：
首先，每种随机过程都对一个信号雨流处理的结果比较
html特效代码雨代码
雨流
低频窄带随机过程
高频窄带随机过程
组合宽带随机过程
采样点数N1
1001
1001
1001
极值点数N2
943
82
941
剩余点数N3
783
38
277
截取点数N4
160
44
624
截取率rate=N4/N2
0.1697
0.5366
0.6992
其次，每种随机过程都对很多个信号雨流处理的统计结果（对100组信号进行统计）:
                                                  rate
低频窄带随机过程  0.0814
频窄带随机过程    0.3006
组合宽带随机过程  0.3178
这样看的话这种雨流计数法对这三种随机过程的适用性：宽带>高频窄带>低频窄带（自己这么想的，不知是否行的通）
三、对剩余载荷的处理
re=[ZSr,ZSr];
[reside,Ba,Bm,Q]=Rainflow(re);
问题三：我的理解是对于剩余载荷residue，在复制一个一样的residue，让它们首尾相接作为一个新的过程再进行雨流处理，得到{residue2+cycles}。我对文章的理解是residue2要和原本的residue一致，但这又如何保证呢？总之，对这里理解的不是很清晰。
四、载荷的重构（将循环插入剩余载荷中）
（1）筛选出要插入循环Q中高于门槛值threshold的循环存于S中
[c,r]=size(Q);     
S=[];
Abs=[];
for i=1:c                         
    if abs(Q(i,2)-Q(i,1))>threshold
        S=[S;Q(i,1),Q(i,2)];
        Abs=[Abs,abs(Q(i,2)-Q(i,1))];
    end
end
（2）对S中循环按幅值大小降序排列
[Abs,index]=sort(Abs);
index=fliplr(index);
k1=S(:,1);S(:,1)=k1(index);
k2=S(:,2);S(:,2)=k2(index);
（3）按幅值由大到小插入载荷序列中（这部分编出来程序运行没有问题，但结果显示似乎并没有插进去，还在寻问题中…）
[c,r]=size(S);   
no1=[];no2=[];
while m<=c
    if S(m,1)<S(m,2)      %对于第m个循环
        n=length(reside);
        rde=[];indexrde=[];
        for i=1:n-1
            if (reside(i)>=S(m,2))&&(reside(i+1)<=S(m,1))      %筛选出第m个循环能够插入的序列        rde=[rde;reside(i),reside(i+1)];
                indexrde=[indexrde;i,i+1];
            end
        end
        [k,g]=size(rde);
        if k==0
          m=m+1;
        else           

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。