Home > wafo > trgauss > mctp2tc.m

mctp2tc

PURPOSE ^

Frequencies of upcrossing troughs and crests using Markov chain of turning points.

SYNOPSIS ^

F=mctp2tc(freqPVR,utc,param,freqPVL)

DESCRIPTION ^

 MCTP2TC  Frequencies of upcrossing troughs and crests using Markov chain of turning points. 
  
   CALL: f_tc = mctp2tc(f_Mm,utc,param); 
  
   where 
  
         f_tc  = the matrix with frequences of upcrossing troughs and crests, 
         f_Mm  = the frequency matrix for the Max2min cycles, 
         utc   = the reference level, 
         param = a vector defining the discretization used to compute f_Mm, 
                 note that f_mM has to be computed on the same grid as f_mM.  
  
   optional call: f_tc = mctp2tc(f_Mm,utc,param,f_mM); 
  
         f_mM  = the frequency matrix for the min2Max cycles.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function F=mctp2tc(freqPVR,utc,param,freqPVL) 
002 %MCTP2TC  Frequencies of upcrossing troughs and crests using Markov chain of turning points. 
003 % 
004 %  CALL: f_tc = mctp2tc(f_Mm,utc,param); 
005 % 
006 %  where 
007 % 
008 %        f_tc  = the matrix with frequences of upcrossing troughs and crests, 
009 %        f_Mm  = the frequency matrix for the Max2min cycles, 
010 %        utc   = the reference level, 
011 %        param = a vector defining the discretization used to compute f_Mm, 
012 %                note that f_mM has to be computed on the same grid as f_mM.  
013 % 
014 %  optional call: f_tc = mctp2tc(f_Mm,utc,param,f_mM); 
015 % 
016 %        f_mM  = the frequency matrix for the min2Max cycles. 
017 % 
018  
019  
020 if nargin<4 
021   freqPVL=freqPVR; 
022 end 
023 if nargin<3   
024    display('too few inputs parameters, stop') 
025    break 
026 end 
027  
028 u=levels(param); 
029 udisc=fliplr(u); 
030 ntc=sum(udisc>=utc); 
031 %ntc1=sum(u<=utc); 
032  
033 %if udisc(ntc) > utc 
034 %   ntc=ntc+1 
035 %end 
036  
037 %if u(ntc1) >= utc 
038 %   ntc1=ntc1-1 
039 %end 
040  
041  
042 if ntc < 2 
043   display('the reference level out of range, stop') 
044     break 
045     end 
046 if ntc > (param(3)-1) 
047   display('the reference level out of range, stop') 
048     break 
049     end 
050  
051 % normalization of frequency matrices 
052  
053 nP=length(freqPVR); 
054 for i=1:nP, 
055     rowsum=sum(freqPVR(i,:)); 
056     if rowsum~=0 
057        freqPVR(i,:)=freqPVR(i,:)/rowsum; 
058     end 
059 end 
060 P=fliplr(freqPVR);  
061  
062 Ph=rot90(fliplr(freqPVL),-1); 
063 for i=1:nP, 
064     rowsum=sum(Ph(i,:)); 
065     if rowsum~=0 
066        Ph(i,:)=Ph(i,:)/rowsum; 
067     end 
068 end 
069 Ph=fliplr(Ph); 
070  
071  
072 n=nP; F=zeros(n,n);  
073  
074 if ntc > n-1 
075   display('index for mean-level out of range, stop') 
076     break 
077     end 
078  
079   
080 %F(1:ntc,1:ntc1)=freqPVL(1:ntc,1:ntc1); 
081 F(1:ntc-1,1:(n-ntc))=freqPVL(1:ntc-1,1:(n-ntc)); 
082  
083  
084 F=cmat2nt(F); 
085  
086 for i=2:ntc, 
087     for j=ntc:n-1, 
088  
089     if i<ntc 
090    
091        Ap=P(i:ntc-1,i+1:ntc); Bp=Ph(i+1:ntc,i:ntc-1); 
092        dim_p=ntc-i; 
093        tempp=zeros(dim_p,1); 
094        I=eye(size(Ap)); 
095        if i==2 
096            e=Ph(i+1:ntc,1); 
097          else 
098            e=sum(Ph(i+1:ntc,1:i-1)')'; 
099        end 
100        if max(abs(e))>1e-10 
101           if dim_p==1 
102          tempp(1,1)=(Ap/(1-Bp*Ap)*e); 
103           else 
104         tempp=Ap*((I-Bp*Ap)\e); 
105           end 
106         end 
107      end 
108  
109      if j>ntc 
110  
111        Am=P(ntc:j-1,ntc+1:j); Bm=Ph(ntc+1:j,ntc:j-1);   
112        dim_m=j-ntc;  
113        tempm=zeros(dim_m,1); 
114        Im=eye(size(Am)); 
115        if j==n-1        
116          em=P(ntc:j-1,n); 
117        else 
118          em=sum(P(ntc:j-1,j+1:n)')'; 
119        end 
120        if max(abs(em))>1e-10 
121           if dim_m==1 
122         tempm(1,1)=(Bm/(1-Am*Bm)*em); 
123           else 
124         tempm=Bm*((Im-Am*Bm)\em); 
125           end 
126        end 
127     end 
128       
129      if (j>ntc)*(i<ntc) 
130     F(i,n-j+1)=F(i,n-j+1)+tempp'*freqPVL(i:ntc-1,n-ntc:-1:n-j+1)*tempm; 
131     F(i,n-j+1)=F(i,n-j+1)+tempp'*freqPVL(i:ntc-1,n-j:-1:1)*ones(n-j,1); 
132     F(i,n-j+1)=F(i,n-j+1)+ones(1,i-1)*freqPVL(1:i-1,n-ntc:-1:n-j+1)*tempm; 
133      end 
134      if (j==ntc)*(i<ntc) 
135     F(i,n-j+1)=F(i,n-j+1)+tempp'*freqPVL(i:ntc-1,n-j:-1:1)*ones(n-j,1); 
136     for k=1:ntc  
137         F(i,n-k+1)=F(i,n-ntc+1); 
138     end 
139      end 
140      if (j>ntc)*(i==ntc) 
141     F(i,n-j+1)=F(i,n-j+1)+ones(1,i-1)*freqPVL(1:i-1,n-ntc:-1:n-j+1)*tempm; 
142     for k=ntc:n 
143     F(k,n-j+1)=F(ntc,n-j+1); 
144     end 
145      end 
146  
147   end 
148 end  
149 F;    
150 F=nt2cmat(F); 
151  
152  
153 %fmax=max(max(F)); 
154  
155 %  contour (u,u,flipud(F),... 
156 %fmax*[0.005 0.01 0.02 0.05 0.1 0.2 0.4 0.6 0.8]) 
157 %  axis([param(1) param(2) param(1) param(2)]) 
158  
159 %  title('Crest-trough density') 
160 %  ylabel('crest'), xlabel('trough')                     
161 %  axis('square') 
162 %if mlver>1, commers, end 
163  
164  
165  
166  
167  
168  
169  
170  
171

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