Home > wafo > wsim > private > fr2res.m

fr2res

PURPOSE ^

Generates a stationary residual from the frequency matrix.

SYNOPSIS ^

r=fr2res(f)

DESCRIPTION ^

 FR2RES Generates a stationary residual from the frequency matrix.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function r=fr2res(f)
002 %FR2RES Generates a stationary residual from the frequency matrix.
003 
004 %  Copyright 1993, Mats Frendahl & Igor Rychlik,
005 %  Dept. of Math. Stat., University of Lund.
006 
007 % Changed by PJ (Pär Johannesson) 16-Feb-2004
008 
009 N=length(f); 
010 
011 % Get index for upper left element > 0 in frequency matrix.
012 
013 [index_i,index_j,number]=getulcc(f);
014 
015 if N+1-index_j-index_i < 4 
016     r=[N+1-index_j index_i];
017     return
018 end
019 
020 left_done=0; right_done=0; success=0; direction=-1;
021 
022 while ~success
023     
024     % Let the residual be (min,max) with levels found above.
025     
026     r=[N+1-index_j index_i];
027     
028     comb=fr2comb(f,r);
029     
030     X=r(1); again=1; number_of_times=0;
031     
032     while again
033         
034         h=0;
035         
036         if r(1)>r(2)
037             % disp('upslope')
038             
039             if r(1)>(r(2)+1)
040                 m=r(1); M=r(2); j=N-m+1; i=m-1:-1:M+1; 
041                 combinations=comb(i,j);
042                 freq=f(i,j);
043                 prob=comb2pro(combinations,freq);
044                 h=min(find(rand<cumsum(prob)))-1;
045             end
046             
047             if h>0
048                 C=m-h;
049                 i=C:m; 
050                 if f(C,j)>0
051                     r=[m-1 C r]; 
052                 else
053                     'Warning 1', break
054                     r=[r(1)-1 r(2:length(r))];
055                 end
056             else
057                 r(1)=r(1)-1;
058             end
059             
060         else
061             
062             % disp('downslope')
063             if r(1)<r(2)-1
064                 m=r(2); M=r(1); i=M; j=M+1:m-1; j=N-j+1;
065                 combinations=comb(i,j)';
066                 freq=f(i,j)';
067                 prob=comb2pro(combinations,freq);
068                 h=min(find(rand<cumsum(prob)))-1;
069             end   
070             
071             if h>0
072                 j=M:M+h; j=N-j+1; 
073                 c0=M+h;
074                 c=N+1-c0;
075                 if f(i,c)>0
076                     r=[M+1 c0 r];
077                 else
078                     'Warning 2', break
079                     r=[r(1)+1 r(2:length(r))];
080                 end
081             else
082                 r(1)=r(1)+1;
083             end
084         end
085         
086         if r(1)==r(2), r=r(2:length(r)); end  
087         
088         if (length(r)>=2), again=1; else again=0; end
089         
090         X=[X r(1)];
091         
092         number_of_times=number_of_times+1;
093         
094         if length(r)<2
095             r=[N+1-index_j index_i];
096             number_of_times=0;
097         end
098         
099         again=number_of_times<100;    
100         
101     end
102     
103     if ~isempty(r)
104         if direction==-1
105             r1=r;
106             left_done=1;
107             direction=1;
108         elseif direction==1
109             r2=r;
110             right_done=1;
111         end
112     end
113     
114     success=left_done&right_done;
115     
116 end
117 
118 r=[r1(1:length(r1)) N+1-index_j  fliplr(r2)];
119 
120 %subplot(2,2,1)
121 %plot(r1)
122 %subplot(2,2,2)
123 %plot(fliplr(r2))
124 %subplot(2,2,3)
125 %plot(r)
126 
127

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