Home > wafo > cycles > rfcfilter.m

rfcfilter

PURPOSE ^

Rainflow filter a signal.

SYNOPSIS ^

[y] = rfcfilter(x,h,def)

DESCRIPTION ^

  RFCFILTER Rainflow filter a signal.
 
  CALL:  [y] = rfcfilter(x,h,def);
  
  Input:
    x   = Signal.   [nx1] OR [nx2]
    h   = Threshold for rainflow filter.
    def = 0: removes cycles with range < h. (default)
          1: removes cycles with range <= h.
  Output:
    y   = Rainflow filtered signal.
 
  Examples:  % 1. Filtered signal y is the turning points of x.
       y = rfcfilter(x,0,1);
             % 2. This removes all rainflow cycles with range less than 0.5.
       y = rfcfilter(x,0.5);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function [y] = rfcfilter(x,h,def)  
002 % RFCFILTER Rainflow filter a signal.
003 %
004 % CALL:  [y] = rfcfilter(x,h,def);
005 % 
006 % Input:
007 %   x   = Signal.   [nx1] OR [nx2]
008 %   h   = Threshold for rainflow filter.
009 %   def = 0: removes cycles with range < h. (default)
010 %         1: removes cycles with range <= h.
011 % Output:
012 %   y   = Rainflow filtered signal.
013 %
014 % Examples:  % 1. Filtered signal y is the turning points of x.
015 %      y = rfcfilter(x,0,1);
016 %            % 2. This removes all rainflow cycles with range less than 0.5.
017 %      y = rfcfilter(x,0.5);
018 
019 % Tested on Matlab 6.0
020 %
021 % History: 
022 % Updated by PJ 18-May-2000
023 %   Help text.
024 % Revised by PJ 12-Jan-2000
025 %   updated for WAFO
026 % Created by PJ (Pär Johannesson) 1999
027 
028   
029 % Check input and output
030 
031 ni = nargin;
032 no = nargout;
033 error(nargchk(2,3,ni));
034 
035 if ni < 3
036   def=[];
037 end
038 
039 % Set default values
040 
041 if isempty(def)
042   def = 0;
043 end
044 
045 % Initiation
046 
047 [n,m] = size(x);
048 y = zeros(n,m);
049 
050 y(1,:) = x(1,:);
051 j = 1;
052 if m == 1
053   t0 = 0;
054   y0 = y(1);
055 else
056   t0 = y(1,1);
057   y0 = y(1,2);
058 end
059 
060 z0 = 0;
061 
062 % The rainflow filter
063 
064 for i = 2:n
065   fpi = y0+h; %fp(y0,h);
066   fmi = y0-h; %fm(y0,h);
067   if m == 1
068     ti = 0;
069     xi = x(i);
070   else
071     ti = x(i,1);
072     xi = x(i,2);
073   end
074   if z0 == 0
075     if def == 0
076       test1 = (xi <= fmi);
077       test2 = (xi >= fpi);
078     else % def == 1
079       test1 = (xi < fmi);
080       test2 = (xi > fpi);
081     end
082     if test1
083       z1 = -1;
084     elseif test2
085       z1 = +1;
086     else
087       z1 = 0;
088     end
089     if z1 == 0
090       t1 = t0;
091       y1 = y0;
092     else
093       t1 = ti;
094       y1 = xi;
095     end
096   else
097     
098     % Which definition?
099     if def == 0
100       test1 = (((z0==+1) & (xi<=fmi))  | ((z0==-1) & (xi<fpi)));
101       test2 = (((z0==+1) & (xi>fmi)) | ((z0==-1) & (xi>=fpi)));
102     else % def == 1
103       test1 = (((z0==+1) & (xi<fmi))  | ((z0==-1) & (xi<=fpi)));
104       test2 = (((z0==+1) & (xi>=fmi)) | ((z0==-1) & (xi>fpi)));
105     end
106     
107     % Update z1
108     if     test1
109       z1 = -1;
110     elseif test2
111       z1 = +1;
112     else
113       warning(['Something wrong, i=' num2str(i)]);
114     end
115     
116     % Update y1
117     if     z1 ~= z0
118       t1 = ti;
119       y1 = xi;
120     elseif z1 == -1
121       % y1 = min([y0 xi]);
122       if y0<xi
123     t1 = t0;
124     y1 = y0;
125       else
126     t1 = ti;
127     y1 = xi;
128       end
129     elseif z1 == +1
130       % y1 = max([y0 xi]);
131       if y0 > xi
132     t1 = t0;
133     y1 = y0;
134       else
135     t1 = ti;
136     y1 = xi;
137       end
138     end
139     
140   end
141   
142   % Update y if y0 is a turning point
143   if abs(z0-z1) == 2
144     j=j+1;
145     if m ==1
146       y(j) = y0;
147     else
148       y(j,:) = [t0 y0];
149     end
150   end
151   
152   % Update t0, y0, z0
153   t0=t1;
154   y0=y1;
155   z0=z1;
156 end
157 
158 % Update y if last y0 is greater than (or equal) threshold
159 if def == 0
160   test = (abs(y0-y(j)) > h);
161 else % def == 1
162   test = (abs(y0-y(j)) >= h);
163 end
164 if test
165   j=j+1;
166   if m ==1
167     y(j) = y0;
168   else
169     y(j,:) = [t0 y0];
170   end
171 end
172 
173 % Truncate y
174 y=y(1:j,:);
175 
176 function uu = fp(u,h)
177   uu = u+h;
178   
179 function uu = fm(u,h)
180   uu = u-h;
181

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003