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
returns the constant acceleration of gravity | |
Numerical integration with the Simpson method | |
Transforms between different types of spectra | |
Toggle Transform between angular frequency and frequency spectrum | |
Create character array (string). | |
Display message and abort function. | |
Get structure field contents. | |
1-D interpolation (table lookup) | |
True for cell array. | |
True for character array (string). | |
True if field is in structure array. | |
Convert string to lowercase. | |
Not-a-Number. | |
Compare strings. | |
Find possible matches for string. |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. | |
Joint (Vcf,Hd) PDF for lineare waves with Ochi-Hubble spectra. | |
Separates the linear component of the Spectrum | |
Simulates a Randomized 2nd order non-linear wave X(t) | |
Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. |
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
Comments or corrections to the WAFO group