Home > wafo > wsim > private > rfc2load_fat.m

rfc2load_fat

PURPOSE ^

Recontructs a load process given the frequency matrix (and residual).

SYNOPSIS ^

[X,res,comb,f]=rfc2load_fat(f,res,num_cc)

DESCRIPTION ^

 RFC2LOAD_FAT  Recontructs a load process given the frequency matrix (and residual).
 
  CALL:  X = rfc2load_fat(f,num_cc,residual);
 
   where
 
         X        = the reconstructed load,
         f        = the frequency matrix for the rainflow count,
         residual = the residual (optional input argument),
                    if left out the program will generate a
                    stationary residual using the frequency 
                    matrix.
         num_cc   = the expected number of cycles,
                    if num_cc=[] the program will continue
                    until the residual is empty,

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [X,res,comb,f]=rfc2load_fat(f,res,num_cc)
002 %RFC2LOAD_FAT  Recontructs a load process given the frequency matrix (and residual).
003 %
004 % CALL:  X = rfc2load_fat(f,num_cc,residual);
005 %
006 %  where
007 %
008 %        X        = the reconstructed load,
009 %        f        = the frequency matrix for the rainflow count,
010 %        residual = the residual (optional input argument),
011 %                   if left out the program will generate a
012 %                   stationary residual using the frequency 
013 %                   matrix.
014 %        num_cc   = the expected number of cycles,
015 %                   if num_cc=[] the program will continue
016 %                   until the residual is empty,
017 
018 %  Copyright 1993, Mats Frendahl & Igor Rychlik,
019 %  Dept. of Math. Stat., University of Lund.
020 
021 % Copyright (c) 2004 by Pär Johannesson
022 
023 % Tested  on Matlab  6.5
024 %
025 % History:
026 % Adapted to WAFO by PJ (Pär Johannesson) 16-Feb-2004
027 %   The function 'rfc2load' originally from FAT (Fatigue Analysis Toolbox)
028 %   FAT is a predecessor of WAFO
029 % Changed by PJ 19-Feb-2004
030 
031 %%%%
032 % Check input arguments
033 
034 ni = nargin;
035 no = nargout;
036 error(nargchk(1,3,ni));
037 
038 if ni<2, res=[]; end
039 if ni<3, num_cc=[]; end
040 
041 % Check if freq. matrix is interger matrix.
042 if ~(round(f)==f)
043     warning('The frequency matrix is not interger matrix.   The matrix will be scaled.')
044     f=floor(100*f);
045 end
046 
047 N=length(f);
048 
049 % If empty residual, then calculate a stat. residual
050 % using the frequency matrix.
051 if isempty(res)
052     res0=fr2res(f);
053     % Don't give a proper residual, modify it!
054     % Correction by PJ
055     [RFC,RFC1,res] = tp2rfc(res0(:));  % Ger proper residual
056     resWAFO = N-res+1;              % Convert to WAFO-def
057     RFMres = dtp2rfm(resWAFO,N,'CS'); % Rainflow matrix of residual
058     RFMres = fliplr(RFMres)';       % Convert to FAT-def
059     f = f-RFMres;                   % Remove cycles in res from rainflow matrix
060 end
061 
062 % Find the largest amplitude. Look up how many there are,
063 % generate this number of large amplitudes, reset the 
064 % frequency matrix to 0 for this kind of amplitude.
065 largest_res=sort([min(min(res)) max(max(res))]);
066 index_i=largest_res(1); index_j=N-largest_res(2)+1;
067 n_large_res=f(index_i,index_j); f(index_i,index_j)=0;
068 xtra_cc=[]; for i=1:n_large_res, xtra_cc=[xtra_cc largest_res]; end
069 
070 r=res(:)';
071 
072 % Find the larges amplitude in the residual and put the 
073 % extra large amplitudes here.
074 undone=1;  
075 for i=1:length(r)-1
076     if undone
077         if ( sum(largest_res==[r(i) r(i+1)])==2 )
078             r=[r(1:i-1) xtra_cc r(i:length(r))];
079             undone=0;
080         elseif ( sum(fliplr(largest_res)==[r(i) r(i+1)])==2 )
081             r=[r(1:i-1) fliplr(xtra_cc) r(i:length(r))];
082             undone=0;
083         end
084     end
085 end
086 
087 
088 % If the number of requested cycles are not empty then find
089 % the number of extected runs to give this number of cycles.
090 if num_cc~=[]
091     % Norm frequency matrix to probability matrix.
092     f=f/sum(sum(f));
093         
094     % Calculate the number of extected runs to get the requested number 
095     % of cycles.
096     num_cc_limit=2*num_cc+4;
097     
098     % Scale the frequency matrix to whole numbers and make it large to 
099     % be stationary.
100     f=floor(max([1000*N 1e6])*f);
101 end
102 
103 % Calculated the combination matrix
104 comb=fr2comb(f,r);
105 
106 % Initiate the residual.
107 X=r(1); 
108 
109 again=1; num_runs=0;
110 
111 while again
112     
113     num_runs=num_runs+1;
114         
115     h=0;
116     
117     if r(1)>r(2)
118         % disp('upslope')
119         
120         if r(1)>(r(2)+1)
121             m=r(1); M=r(2); j=N-m+1; i=m-1:-1:M+1; 
122             combinations=comb(i,j);
123             freq=f(i,j);
124             prob=comb2pro(combinations,freq);
125             h=min(find(rand<cumsum(prob)))-1;
126         end
127         
128         if h>0
129             C=m-h;
130             i=C:m; 
131             if f(C,j)>0
132                 f(C,j)=f(C,j)-1;
133                 i=(C+1):r(1);
134                 j=N+1-r(1);
135                 comb(i,j)=comb(i,j)-1;
136                 r=[m-1 C r]; 
137             else
138                 'Warning 1', break
139                 r=[r(1)-1 r(2:length(r))];
140             end
141         else
142             j=N+1-r(1);
143             if r(1)>(r(2)+1)
144                 i=(r(2)+1):(r(1)-1);
145                 comb(i,j)=comb(i,j)-1;
146             end
147             r(1)=r(1)-1;
148             
149         end
150         
151     else
152         
153         % disp('downslope')
154         if r(1)<r(2)-1
155             m=r(2); M=r(1); i=M; j=M+1:m-1; j=N-j+1;
156             combinations=comb(i,j)';
157             freq=f(i,j)';
158             prob=comb2pro(combinations,freq);
159             h=min(find(rand<cumsum(prob)))-1;
160         end   
161         
162         if h>0
163             j=M:M+h; j=N-j+1; 
164             c0=M+h;
165             c=N+1-c0;
166             if f(i,c)>0
167                 f(i,c)=f(i,c)-1;
168                 i=r(1);
169                 j=N+1-((i+h-1:-1:r(1)));
170                 comb(i,j)=comb(i,j)-1;
171                 r=[M+1 c0 r];
172             else
173                 'Warning 2', break
174                 r=[r(1)+1 r(2:length(r))];
175             end
176         else
177             i=r(1);
178             j=N+1-(r(1):(r(2)-1));
179             comb(i,j)=comb(i,j)-1;
180             r(1)=r(1)+1;
181                                     
182         end
183     end
184     
185     if r(1)==r(2), r=r(2:length(r)); end  
186     
187     if (length(r)>=2)
188         again=1;
189     else
190         disp('   The residual is empty. Program will terminate.')
191         again=0;
192     end
193     
194     X=[X r(1)];
195     
196     n_tmp=length(X);
197     
198     % This is an equivalent of shave, but faster
199     
200     if n_tmp>2
201         if (X(n_tmp-2)<X(n_tmp-1))&(X(n_tmp-1)<X(n_tmp))
202             X=[X(1:n_tmp-2) X(n_tmp)];
203         elseif (X(n_tmp-2)>X(n_tmp-1))&(X(n_tmp-1)>X(n_tmp))
204             X=[X(1:n_tmp-2) X(n_tmp)];
205         end
206     end
207     
208     comb=fliplr(triu(fliplr(comb),1));
209     
210     % If the number of requested cycles are not empty then check
211     % if the number of step (length of load) exceeds the expected 
212     % number.
213     if num_cc~=[]
214         if length(X)>num_cc_limit
215             disp('   The number of requested cycles has been reached.')
216             disp('   Program will terminate.')
217             again=0;
218         end
219     end
220     
221 end
222 
223 X = X';
224

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