Home > wafo > spec > spec2char.m

spec2char

PURPOSE ^

Evaluates spectral characteristics and their covariance

SYNOPSIS ^

[ch,R1,chtext,R]=spec2char(S,fact,T)

DESCRIPTION ^

 SPEC2CHAR  Evaluates spectral characteristics and their covariance
 
  CALL: [ch R,chtext] = spec2char(S,fact,T)
 
        ch = vector of spectral characteristics
        R  = matrix of the corresponding covariances given T
    chtext = a cellvector of strings describing the elements of ch, see example.
        S  = spectral struct with angular frequency
      fact = vector with factor integers or a string or
             a cellarray of strings, see below.(default [1])
        T  = recording time (sec) (default 1200 sec = 20 min)
 
  If input spectrum is of wave number type, output are factors for
  corresponding 'k1D', else output are factors for 'freq'.
  Input vector 'factors' correspondence:
     1 Hm0   = 4*sqrt(m0)                              Significant wave height
     2 Tm01  = 2*pi*m0/m1                              Mean wave period
     3 Tm02  = 2*pi*sqrt(m0/m2)                        Mean zero-crossing period
     4 Tm24  = 2*pi*sqrt(m2/m4)                        Mean period between maxima
     5 Tm_10 = 2*pi*m_1/m0                             Energy period
     6 Tp    = 2*pi/{w | max(S(w))}                    Peak period  
     7 Ss    = 2*pi*Hm0/(g*Tm02^2)                     Significant wave steepness
     8 Sp    = 2*pi*Hm0/(g*Tp^2)                       Average wave steepness
     9 Ka    = abs(int S(w)*exp(i*w*Tm02) dw ) /m0     Groupiness parameter
    10 Rs    = (S(0.092)+S(0.12)+S(0.15)/(3*max(S(w))) Quality control parameter
    11 Tp1   = 2*pi*int S(w)^4 dw                      Peak Period (robust estimate for Tp) 
               ------------------
               int w*S(w)^4 dw 
 
    12 alpha = m2/sqrt(m0*m4)                          Irregularity factor
    13 eps2  = sqrt(m0*m2/m1^2-1)                      Narrowness factor
    14 eps4  = sqrt(1-m2^2/(m0*m4))=sqrt(1-alpha^2)    Broadness factor
    15 Qp    = (2/m0^2)int_0^inf w*S(w)^2 dw           Peakedness factor
 
  Order of output is same as order in 'factors'
  The covariances are computed with a Taylor expansion technique
  and is currently only available for factors 1, 2, and 3. Variances
  are also available for factors 4,5,7,12,13,14 and 15 
  
  Quality control:
   Critical value for quality control parameter Rs is Rscrit = 0.02
   for surface displacement records and Rscrit=0.0001 for records of
   surface acceleration or slope. If Rs > Rscrit then probably there 
   are something wrong with the lower frequency part of S.
 
   Ss may be used as an indicator of major malfunction, by checking that
   it is in the range of 1/20 to 1/16 which is the usual range for
   locally generated wind seas. 
 
  Examples:
    S      = demospec;
    [ch R,txt] = spec2char(S,[1 2 3]);    % fact a vector of integers
    cat(2,char(txt),repmat(' = ',3,1),num2str(ch')) 
    ch = num2cell(ch);
    ch0 = cell2struct(ch,txt,2);          % Make a structure 
    [txt{2,:}]=deal(','); txt{2,end}=' ';
    eval(['[' txt{:} '] = deal(ch{:})'])  % Assign values to variables
    ch1    = spec2char(S,'Ss');           % fact a string
    ch2    = spec2char(S,{'Tp','Tp1'});   % fact a cellarray of strings
 
  See also  dspec2char, spec2bw, spec2mom, simpson

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [ch,R1,chtext,R]=spec2char(S,fact,T)
002 %SPEC2CHAR  Evaluates spectral characteristics and their covariance
003 %
004 % CALL: [ch R,chtext] = spec2char(S,fact,T)
005 %
006 %       ch = vector of spectral characteristics
007 %       R  = matrix of the corresponding covariances given T
008 %   chtext = a cellvector of strings describing the elements of ch, see example.
009 %       S  = spectral struct with angular frequency
010 %     fact = vector with factor integers or a string or
011 %            a cellarray of strings, see below.(default [1])
012 %       T  = recording time (sec) (default 1200 sec = 20 min)
013 %
014 % If input spectrum is of wave number type, output are factors for
015 % corresponding 'k1D', else output are factors for 'freq'.
016 % Input vector 'factors' correspondence:
017 %    1 Hm0   = 4*sqrt(m0)                              Significant wave height
018 %    2 Tm01  = 2*pi*m0/m1                              Mean wave period
019 %    3 Tm02  = 2*pi*sqrt(m0/m2)                        Mean zero-crossing period
020 %    4 Tm24  = 2*pi*sqrt(m2/m4)                        Mean period between maxima
021 %    5 Tm_10 = 2*pi*m_1/m0                             Energy period
022 %    6 Tp    = 2*pi/{w | max(S(w))}                    Peak period  
023 %    7 Ss    = 2*pi*Hm0/(g*Tm02^2)                     Significant wave steepness
024 %    8 Sp    = 2*pi*Hm0/(g*Tp^2)                       Average wave steepness
025 %    9 Ka    = abs(int S(w)*exp(i*w*Tm02) dw ) /m0     Groupiness parameter
026 %   10 Rs    = (S(0.092)+S(0.12)+S(0.15)/(3*max(S(w))) Quality control parameter
027 %   11 Tp1   = 2*pi*int S(w)^4 dw                      Peak Period (robust estimate for Tp) 
028 %              ------------------
029 %              int w*S(w)^4 dw 
030 %
031 %   12 alpha = m2/sqrt(m0*m4)                          Irregularity factor
032 %   13 eps2  = sqrt(m0*m2/m1^2-1)                      Narrowness factor
033 %   14 eps4  = sqrt(1-m2^2/(m0*m4))=sqrt(1-alpha^2)    Broadness factor
034 %   15 Qp    = (2/m0^2)int_0^inf w*S(w)^2 dw           Peakedness factor
035 %
036 % Order of output is same as order in 'factors'
037 % The covariances are computed with a Taylor expansion technique
038 % and is currently only available for factors 1, 2, and 3. Variances
039 % are also available for factors 4,5,7,12,13,14 and 15 
040 % 
041 % Quality control:
042 %  Critical value for quality control parameter Rs is Rscrit = 0.02
043 %  for surface displacement records and Rscrit=0.0001 for records of
044 %  surface acceleration or slope. If Rs > Rscrit then probably there 
045 %  are something wrong with the lower frequency part of S.
046 %
047 %  Ss may be used as an indicator of major malfunction, by checking that
048 %  it is in the range of 1/20 to 1/16 which is the usual range for
049 %  locally generated wind seas. 
050 %
051 % Examples:
052 %   S      = demospec;
053 %   [ch R,txt] = spec2char(S,[1 2 3]);    % fact a vector of integers
054 %   cat(2,char(txt),repmat(' = ',3,1),num2str(ch')) 
055 %   ch = num2cell(ch);
056 %   ch0 = cell2struct(ch,txt,2);          % Make a structure 
057 %   [txt{2,:}]=deal(','); txt{2,end}=' ';
058 %   eval(['[' txt{:} '] = deal(ch{:})'])  % Assign values to variables
059 %   ch1    = spec2char(S,'Ss');           % fact a string
060 %   ch2    = spec2char(S,{'Tp','Tp1'});   % fact a cellarray of strings
061 %
062 % See also  dspec2char, spec2bw, spec2mom, simpson
063 
064 % References:
065 % Krogstad, H.E., Wolf, J., Thompson, S.P., and Wyatt, L.R. (1999)
066 % 'Methods for intercomparison of wave measurements'
067 % Coastal Enginering, Vol. 37, pp. 235--257
068 %
069 % Krogstad, H.E. (1982)
070 % 'On the covariance of the periodogram'
071 % Journal of time series analysis, Vol. 3, No. 3, pp. 195--207
072 %
073 % Tucker, M.J. (1993)
074 % 'Recommended standard for wave data sampling and near-real-time processing'
075 % Ocean Engineering, Vol.20, No.5, pp. 459--474
076 %
077 % Young, I.R. (1999)
078 % "Wind generated ocean waves"
079 % Elsevier Ocean Engineering Book Series, Vol. 2, pp 239
080 
081 % Tested on: Matlab 5.2
082 % History: 
083 % revised pab jan2004
084 %  -added todo comments
085 % revised pab 25.06.2001
086 % -added chtext to output+more examples
087 % revised pab 07.07.2000
088 %  -fixed a bug for variance of Tm02
089 %  - added variance for Tm24
090 % revised pab 22.06.2000
091 %  - added alpha, eps2,eps4,Qp
092 %  - added the possibility that fact is a string or a cellarray of strings 
093 % revised pab 23.05.2000
094 %  - added ttspec call
095 %  revised by es 23.05.00, do not call spec2spec if already .type='freq'
096 % Revised pab 06.03.2000
097 % -updated header info
098 %by pab 16.02.2000
099   
100   
101 % TODO % Need more checking on computing the variances for Tm24,alpha, eps2 and eps4 
102 % TODO % Covariances between Tm24,alpha, eps2 and eps4 variables are also needed
103 
104 tfact = {'Hm0', 'Tm01', 'Tm02', 'Tm24', 'Tm_10','Tp','Ss', 'Sp', 'Ka', ...
105       'Rs', 'Tp1','alpha','eps2','eps4','Qp'} ;
106 
107 if nargin<2|isempty(fact)
108   nfact = 1;
109 elseif iscell(fact)|ischar(fact)
110   if ischar(fact), fact = {fact}; end
111   N     = length(fact(:));
112   nfact = zeros(1,N);
113   ltfact = char(lower(tfact));
114   for ix=1:N,
115     ind = strmatch(lower(fact(ix)),ltfact,'exact');
116     if length(ind)==1,
117       nfact(ix)=ind;
118     else
119       error(['Not a valid factor: ' fact{ix}])
120     end
121   end 
122 else
123   nfact = fact;
124 end
125 if any(nfact>15 | nfact<1)
126   error('Factor outside range (1,...,15)')
127 end
128 
129 if nargin<3|isempty(T)
130   T = 1200; % recording time default 1200 secs (=20 minutes)
131 end
132 
133 if isfield(S,'k')
134   S=spec2spec(S,'k1d');
135   vari='k';
136 else 
137   if ~strcmp(lower(S.type),'freq'),S = spec2spec(S,'freq');end
138   S = ttspec(S,'w');
139   vari = 'w';
140 end
141 
142 f  = getfield(S,vari);
143 f  = f(:);
144 S1 = S.S(:);
145 %m=spec2mom(S,4,[],0)./[ (2*pi).^[0:4] ]; % moments corresponding to freq
146 % in Hz
147 m = simpson(f,[S1 f.*S1 f.^2.*S1 f.^3.*S1 f.^4.*S1])./[ (2*pi).^[0:4] ];
148 
149 ind  = find(f>0);
150 m(6) = simpson(f(ind),S1(ind)./f(ind))*2*pi;  % = m_1
151 m_10 = simpson(f(ind),S1(ind).^2./f(ind))*(2*pi)^2/T;    % = COV(m_1,m0|T=t0)
152 m_11 = simpson(f(ind),S1(ind).^2./f(ind).^2)*(2*pi)^3/T; % = COV(m_1,m_1|T=t0)
153 
154 
155 %      Hm0        Tm01        Tm02             Tm24         Tm_10
156 Hm0  = 4*sqrt(m(1)); 
157 Tm01 = m(1)/m(2); 
158 Tm02 = sqrt(m(1)/m(3)); 
159 Tm24 = sqrt(m(3)/m(5)); 
160 Tm_10= m(6)/m(1);
161 
162 Tm12 = m(2)/m(3);
163 
164 g    = gravity;
165 [maxS ind] = max(S1);
166 Tp   = 2*pi/f(ind);                                   % peak period /length
167 Ss   = 2*pi*Hm0/g/Tm02^2;                             % Significant wave steepness
168 Sp   = 2*pi*Hm0/g/Tp^2;                               % Average wave steepness 
169 Ka   = abs(simpson(f,S1.*exp(sqrt(-1)*f*Tm02)))/m(1); % groupiness factor
170 
171 % Quality control parameter 
172 % critical value is approximately 0.02 for surface displacement records
173 % If Rs>0.02 then there are something wrong with the lower frequency part 
174 % of S.
175 Rs   = sum(interp1(f,S1,[0.0146 0.0195 0.0244]*2*pi,'linear'))/3/maxS; 
176 Tp2  = 2*pi*simpson(f,S1.^4)/simpson(f,f.*S1.^4);
177 
178 
179 alpha1 = Tm24/Tm02;                 % m(3)/sqrt(m(1)*m(5));
180 eps2   = sqrt(Tm01/Tm12-1);         % sqrt(m(1)*m(3)/m(2)^2-1);
181 eps4   = sqrt(1-alpha1^2);          % sqrt(1-m(3)^2/m(1)/m(5));
182 Qp     = 2/m(1)^2*simpson(f,f.*S1.^2);
183 
184 
185 
186 
187 ch = [Hm0 Tm01 Tm02 Tm24 Tm_10 Tp Ss Sp Ka Rs Tp2 alpha1 eps2 eps4 Qp];
188 
189 
190 % Select the appropriate values
191 ch     = ch(nfact);
192 chtext = tfact(nfact);
193 
194 if nargout>1,
195   % covariance between the moments:
196   %COV(mi,mj |T=t0) = int f^(i+j)*S(f)^2 df/T
197   ONE = ones(size(f));
198   S2 = S1.^2;
199   mij = simpson(f,[ONE f f.^2  f.^3 f.^4 f.^5 f.^6 f.^7 f.^8].*S2(:,ones(1,9)))/T./[ (2*pi).^[-1:7] ];
200 %   mij = simpson(f,[S1.^2 f.*S1.^2 (f.*S1).^2  f.^3.*S1.^2 (f.^2.*S1).^2 ...
201 %    f.*(f.^2.*S1).^2 (f.^3.*S1).^2 f.*(f.^3.*S1).^2 (f.^4.*S1).^2])/T./[ (2*pi).^[-1:7] ]
202 
203 % and the corresponding variances for
204 %{'hm0', 'tm01', 'tm02', 'tm24', 'tm_10','tp','ss', 'sp', 'ka', 'rs', 'tp1','alpha','eps2','eps4','qp'}  
205   R = [4*mij(1)/m(1) ...
206     mij(1)/m(2)^2-2*m(1)*mij(2)/m(2)^3+m(1)^2*mij(3)/m(2)^4 ...
207     0.25*(mij(1)/(m(1)*m(3))-2*mij(3)/m(3)^2+m(1)*mij(5)/m(3)^3) ...
208     0.25*(mij(5)/(m(3)*m(5))-2*mij(7)/m(5)^2+m(3)*mij(9)/m(5)^3) ...
209     m_11/m(1)^2+(m(6)/m(1)^2)^2*mij(1)-2*m(6)/m(1)^3*m_10,...
210     NaN,...
211     (8*pi/g)^2*(m(3)^2/(4*m(1)^3)*mij(1)+mij(5)/m(1)-m(3)/m(1)^2*mij(3)),...
212     NaN*ones(1,4),...
213     m(3)^2*mij(1)/(4*m(1)^3*m(5))+mij(5)/(m(1)*m(5))+mij(9)*m(3)^2/(4*m(1)*m(5)^3)-...
214     m(3)*mij(3)/(m(1)^2*m(5))+m(3)^2*mij(5)/(2*m(1)^2*m(5)^2)-m(3)*mij(7)/m(1)/m(5)^2,...
215     (m(3)^2*mij(1)/4+(m(1)*m(3)/m(2))^2*mij(3)+m(1)^2*mij(5)/4-m(3)^2*m(1)*mij(2)/m(2)+...
216         m(1)*m(3)*mij(3)/2-m(1)^2*m(3)/m(2)*mij(4))/eps2^2/m(2)^4,...
217     (m(3)^2*mij(1)/(4*m(1)^2)+mij(5)+m(3)^2*mij(9)/(4*m(5)^2)-m(3)*mij(3)/m(1)+....
218     m(3)^2*mij(5)/(2*m(1)*m(5))-m(3)*mij(7)/m(5))*m(3)^2/(m(1)*m(5)*eps4)^2,...
219     NaN];
220  
221   % and covariances by a taylor expansion technique:
222   % Cov(Hm0,Tm01) Cov(Hm0,Tm02) Cov(Tm01,Tm02)
223   S0 = [ 2/(sqrt(m(1))*m(2))*(mij(1)-m(1)*mij(2)/m(2)),...
224     1/sqrt(m(3))*(mij(1)/m(1)-mij(3)/m(3)),......
225     1/(2*m(2))*sqrt(m(1)/m(3))*(mij(1)/m(1)-mij(3)/m(3)-mij(2)/m(2)+m(1)*mij(4)/(m(2)*m(3)))];
226   tmp = NaN;
227   R1  = tmp(ones(15,15));
228   
229   for ix=1:length(R), R1(ix,ix) = R(ix);  end
230   
231   
232   R1(1,2:3)   = S0(1:2);
233   R1(2,3)     = S0(3);
234   for ix = 1:2, %make lower triangular equal to upper triangular part
235     R1(ix+1:3,ix) = R1(ix,ix+1:3).';
236   end
237   
238   R = R(nfact);
239   R1= R1(nfact,nfact);
240 end
241 
242  % Needs further checking:
243  % Var(Tm24)= 0.25*(mij(5)/(m(3)*m(5))-2*mij(7)/m(5)^2+m(3)*mij(9)/m(5)^3) ...
244 
245 
246 
247

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