Home > wafo > spec > spec2mom.m

spec2mom

PURPOSE ^

Calculates spectral moments from spectrum

SYNOPSIS ^

[m,mtext] = spec2mom(S,nr,vari,even)

DESCRIPTION ^

 SPEC2MOM Calculates spectral moments from spectrum
 
  CALL:  [m,mtext] = spec2mom(S,nr,variables,even)
 
    m    = vector of moments
    mtext= a cell array of strings describing the elements of m, see below
    S    = spectrum (struct)
    nr   = order of moments (default 2, maximum 4)
    variables = variables in model, optional when two-dim.spectrum,
                string with x and/or y and/or t  (default 'xt')
    even = 0 for all moments, 1 for only even orders (default 1)
 
  Calculates spectral moments of up to order four by use of
  Simpson-integration.
 
        //
  m_jkl=|| k1^j*k2^k*w^l S(w,th) dw dth 
        //
 
  where k1=w^2/gravity*cos(th-phi),  k2=w^2/gravity*sin(th-phi)
  and phi is the angle of the rotation in S.phi. If the spectrum 
  has field .g, gravity is replaced by S.g.
   
  The strings in output mtext have the same position in the cell array
  as the corresponding numerical value has in output m
  Notation in mtext: 'm0' is the variance,
                     'mx' is the first-order moment in x,
                    'mxx' is the second-order moment in x,
                    'mxt' is the second-order cross moment between x and t,
                  'myyyy' is the fourth-order moment in y
                          etc.
  For the calculation of moments see Baxevani et al.
  Example:
       S=demospec('dir')
       [m,mtext]=spec2mom(S,2,'xyt')

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [m,mtext] = spec2mom(S,nr,vari,even)
002 %SPEC2MOM Calculates spectral moments from spectrum
003 %
004 % CALL:  [m,mtext] = spec2mom(S,nr,variables,even)
005 %
006 %   m    = vector of moments
007 %   mtext= a cell array of strings describing the elements of m, see below
008 %   S    = spectrum (struct)
009 %   nr   = order of moments (default 2, maximum 4)
010 %   variables = variables in model, optional when two-dim.spectrum,
011 %               string with x and/or y and/or t  (default 'xt')
012 %   even = 0 for all moments, 1 for only even orders (default 1)
013 %
014 % Calculates spectral moments of up to order four by use of
015 % Simpson-integration.
016 %
017 %       //
018 % m_jkl=|| k1^j*k2^k*w^l S(w,th) dw dth 
019 %       //
020 %
021 % where k1=w^2/gravity*cos(th-phi),  k2=w^2/gravity*sin(th-phi)
022 % and phi is the angle of the rotation in S.phi. If the spectrum 
023 % has field .g, gravity is replaced by S.g.
024 %  
025 % The strings in output mtext have the same position in the cell array
026 % as the corresponding numerical value has in output m
027 % Notation in mtext: 'm0' is the variance,
028 %                    'mx' is the first-order moment in x,
029 %                   'mxx' is the second-order moment in x,
030 %                   'mxt' is the second-order cross moment between x and t,
031 %                 'myyyy' is the fourth-order moment in y
032 %                         etc.
033 % For the calculation of moments see Baxevani et al.
034 % Example:
035 %      S=demospec('dir')
036 %      [m,mtext]=spec2mom(S,2,'xyt')
037 
038 % References
039 %  Baxevani A. et al. (2001) 
040 %  Velocities for Random Surfaces 
041 %
042 %
043 % Tested on: Matlab 6.0
044 % Tested on: Matlab 5.3
045 % History:
046 % Revised by I.R. 04.04.2001: Introducing the rotation angle phi.
047 % Revised by A.B. 23.05.2001: Correcting 'mxxyy' and introducing
048 % 'mxxyt','mxyyt' and 'mxytt'.
049 % Revised by A.B. 21.10.2001: Correcting 'mxxyt'.
050 % Revised by A.B. 21.10.2001: Adding odd-order moments.
051 % By es 27.08.1999
052 
053 
054 if nargin<2|isempty(nr)
055   nr=2;
056 end
057 if nargin<4|isempty(even)
058   even=1;
059 end
060 even=logical(even);
061 
062 if strcmpi(S.type,'freq')| strcmpi(S.type,'enc')|...
063       strcmpi(S.type,'k1d') %i.e. one-dim spectrum
064   if isfield(S,'w')
065     f=S.w(:);
066     vari='t';
067   elseif isfield(S,'k')
068     f=S.k(:);
069     vari='x';
070   else
071     f=2*pi*S.f(:);
072     S.S=S.S/2/pi;
073     vari='t';
074   end
075   S1=abs(S.S(:));
076   m=simpson(f,S1);
077   mtext={'m0'};
078   if nr>0
079     if ~even
080       m=[m simpson(f, f.*S1)];
081       mtext(end+1)={'mi'};
082     end
083     if nr>1
084       m=[m simpson(f,f.^2.*S1)];
085       mtext(end+1)={'mii'};
086       if nr>2
087         if ~even
088           m=[m simpson(f,f.^3.*S1)];
089           mtext(end+1)={'miii'};
090         end
091         if nr>3
092           m=[m simpson(f,f.^4.*S1)];
093           mtext(end+1)={'miiii'};
094         end
095       end
096     end
097     mtext(end+1)={strcat(' where i=',vari)};
098   end
099 
100 
101 
102 
103 
104 
105 else %two-dim spectrum
106   if (nargin<3|isempty(vari))&nr<=1;
107     vari='x';% set default value
108     Nv=1;
109     elseif nargin<3|isempty(vari);
110     vari='xt';
111     Nv=2;
112     else% secure the mutual order ('xyt')
113     vari=sort(lower(vari));
114     Nv=length(vari);
115     
116     if ( (vari(1)=='t') & (Nv>1)),
117       vari=[vari(2:Nv), 't'];
118     end 
119   end
120   S1=S;
121   %%%%to be removed
122   %%%%%%%
123   %if Nv==1
124    % switch vari
125     %  case 'x'
126         %S1=spec2spec(S,'k1d');
127      % case 'y'
128 %       S1=spec2spec(S,'k1d',-pi/2); % funkar detta??? -pi/2 ???
129     %  case 't'
130 %       S1=spec2spec(S,'freq');
131 
132 %end
133  %   [m,mtext]=spec2mom(S1,nr,vari,even);
134  % else % Nv>1
135    
136     %%to be removed
137     if ~strcmpi(S.type(end-2:end),'dir') 
138       S1=spec2spec(S,strcat(S.type(1:end-3),'dir'));
139     end
140     if ~isfield(S,'phi') | isempty(S.phi), S1.phi=0;
141       
142     end
143     th=S1.theta(:)-S1.phi;
144     Sw=simpson(th,S1.S).'; % integral S(w,th) dth
145     m=simpson(S1.w,Sw);    % the variance
146     mtext={'m0'};
147     
148     if nr>0
149       w=S1.w(:);
150       nw=length(w);
151       if strcmpi(vari(1),'x')
152         Sc=simpson(th,S1.S.*(cos(th)*ones(1,nw))).';
153         % integral S*cos(th) dth
154       end
155       if strcmpi(vari(1),'y')
156         Ss=simpson(th,S1.S.*(sin(th)*ones(1,nw))).';
157         % integral S*sin(th) dth
158         if strcmpi(vari(1),'x')
159         Sc=simpson(th,S1.S.*(cos(th)*ones(1,nw))).';
160         end
161       end
162       if ~isfield(S1,'g')
163         S1.g=gravity;
164       end
165       kx=w.^2/S1.g(1); % maybe different normalization in x and y => diff. g
166       ky=w.^2/S1.g(end);
167       
168       if Nv>=1
169         switch vari
170           case 'x'
171             vec=[kx.*Sc];
172             mtext(end+1)={'mx'};
173           case 'y'
174             vec=[ky.*Ss];
175             mtext(end+1)={'my'};
176           case 't'
177             vec=[w.*Sw];
178            mtext(end+1)={'mt'}; 
179         end 
180       else
181         vec=[kx.*Sc ky.*Ss w*Sw];
182         mtext(end+(1:3))={'mx', 'my', 'mt'};
183       end
184       if nr>1
185       if strcmpi(vari(1),'x')
186         Sc=simpson(th,S1.S.*(cos(th)*ones(1,nw))).';
187         % integral S*cos(th) dth
188         Sc2=simpson(th,S1.S.*(cos(th).^2*ones(1,nw))).';
189         % integral S*cos(th)^2 dth
190       end
191       if strcmpi(vari(1),'y')|strcmpi(vari(2),'y')
192         Ss=simpson(th,S1.S.*(sin(th)*ones(1,nw))).';
193         % integral S*sin(th) dth
194         Ss2=simpson(th,S1.S.*(sin(th).^2*ones(1,nw))).';
195         % integral S*sin(th)^2 dth
196         if strcmpi(vari(1),'x')
197           Scs=simpson(th,S1.S.*((cos(th).*sin(th))*ones(1,nw))).';
198           % integral S*cos(th)*sin(th) dth
199         end
200       end
201       if ~isfield(S1,'g')
202         S1.g=gravity;
203       end
204         
205       if Nv==2
206         switch vari
207           case 'xy'
208             vec=[kx.*Sc ky.*Ss kx.^2.*Sc2 ky.^2.*Ss2 kx.*ky.*Scs];
209             mtext(end+(1:5))={'mx','my','mxx', 'myy', 'mxy'};
210           case 'xt'
211             vec=[kx.*Sc w.*Sw kx.^2.*Sc2 w.^2.*Sw kx.*w.*Sc];
212             mtext(end+(1:5))={'mx','mt','mxx', 'mtt', 'mxt'};
213           case 'yt'
214             vec=[ky.*Ss w.*Sw ky.^2.*Ss2 w.^2.*Sw ky.*w.*Ss];
215             mtext(end+(1:5))={'my','mt','myy', 'mtt', 'myt'};
216         end
217       else
218         vec=[kx.*Sc ky.*Ss w.*Sw kx.^2.*Sc2 ky.^2.*Ss2  w.^2.*Sw kx.*ky.*Scs kx.*w.*Sc ky.*w.*Ss];
219         mtext(end+(1:9))={'mx','my','mt','mxx', 'myy', 'mtt', 'mxy', 'mxt', 'myt'};
220       end
221       if nr>3
222         if strcmpi(vari(1),'x')
223           Sc3=simpson(th,S1.S.*(cos(th).^3*ones(1,nw))).';
224           % integral S*cos(th)^3 dth
225           Sc4=simpson(th,S1.S.*(cos(th).^4*ones(1,nw))).';
226           % integral S*cos(th)^4 dth
227         end
228         if strcmpi(vari(1),'y')|strcmpi(vari(2),'y')
229           Ss3=simpson(th,S1.S.*(sin(th).^3*ones(1,nw))).';
230           % integral S*sin(th)^3 dth
231           Ss4=simpson(th,S1.S.*(sin(th).^4*ones(1,nw))).';
232           % integral S*sin(th)^4 dth
233           if strcmpi(vari(1),'x')  %both x and y
234             Sc2s=simpson(th,S1.S.*((cos(th).^2.*sin(th))*ones(1,nw))).';
235             % integral S*cos(th)^2*sin(th) dth
236             Sc3s=simpson(th,S1.S.*((cos(th).^3.*sin(th))*ones(1,nw))).';
237             % integral S*cos(th)^3*sin(th) dth
238             Scs2=simpson(th,S1.S.*((cos(th).*sin(th).^2)*ones(1,nw))).';
239             % integral S*cos(th)*sin(th)^2 dth
240             Scs3=simpson(th,S1.S.*((cos(th).*sin(th).^3)*ones(1,nw))).';
241             % integral S*cos(th)*sin(th)^3 dth
242             Sc2s2=simpson(th,S1.S.*((cos(th).^2.*sin(th).^2)*ones(1,nw))).';
243             % integral S*cos(th)^2*sin(th)^2 dth
244           end
245         end
246         if Nv==2
247           switch vari
248             case 'xy'
249               vec=[vec kx.^4.*Sc4 ky.^4.*Ss4 kx.^3.*ky.*Sc3s ...
250                     kx.^2.*ky.^2.*Sc2s2 kx.*ky.^3.*Scs3];
251               mtext(end+(1:5))={'mxxxx','myyyy','mxxxy','mxxyy','mxyyy'};
252             case 'xt'
253               vec=[vec kx.^4.*Sc4 w.^4.*Sw kx.^3.*w.*Sc3 ...
254                     kx.^2.*w.^2.*Sc2 kx.*w.^3.*Sc];
255               mtext(end+(1:5))={'mxxxx','mtttt','mxxxt','mxxtt','mxttt'};
256             case 'yt'
257               vec=[vec ky.^4.*Ss4 w.^4.*Sw ky.^3.*w.*Ss3 ...
258                     ky.^2.*w.^2.*Ss2 ky.*w.^3.*Ss];
259               mtext(end+(1:5))={'myyyy','mtttt','myyyt','myytt','myttt'};
260           end
261         else
262           vec=[vec kx.^4.*Sc4 ky.^4.*Ss4 w.^4.*Sw kx.^3.*ky.*Sc3s ...
263                kx.^2.*ky.^2.*Sc2s2 kx.*ky.^3.*Scs3 kx.^3.*w.*Sc3 ...
264                kx.^2.*w.^2.*Sc2 kx.*w.^3.*Sc ky.^3.*w.*Ss3 ...
265                ky.^2.*w.^2.*Ss2 ky.*w.^3.*Ss kx.^2.*ky.*w.*Sc2s ...
266                kx.*ky.^2.*w.*Scs2 kx.*ky.*w.^2.*Scs];
267           mtext(end+(1:15))={'mxxxx','myyyy','mtttt','mxxxy','mxxyy',...
268           'mxyyy','mxxxt','mxxtt','mxttt','myyyt','myytt','myttt','mxxyt','mxyyt','mxytt'};  
269 
270         end % if Nv==2 ... else ...
271       end % if nr>3      
272       end % if nr>1
273       m=[m simpson(w,vec)];
274     end % if nr>0
275   %  end; %%if Nv==1... else...    to be removed
276 end % ... else two-dim spectrum 
277

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