DSPEC2CHAR Evaluates directional spectral characteristics CALL: [ch,chtext] = dspec2char(S,fact); ch = a cell vector of spectral characteristics chtext = a cell vector of strings describing the elements of ch, see example. S = Directional spectral density struct with angular frequency fact = vector of factor integers or a string or a cellarray of strings, see below.(default [1]) DSPEC2CHAR assumes S is a directional spectrum S(w,theta)= S(w)*D(theta,w). If input spectrum is of wave number type, it is transformed into a directional spectrum before the calculations. For many of the parameters the Fourier series expansion of D(theta,w) is used: M-1 D(theta(i)) = {1 + 2*sum [an*cos(n*theta(i)) + bn*sin(n*theta(i))]}/(2*pi) n=1 where C1 = sqrt(a1^2+b1^2) and C2 = sqrt(a2^2+b2^2) Order of output is same as order in 'factors'. Input vector 'factors' correspondence: Factors calculated at every frequency (frequency dependent parameters): 1 FMdir = atan2(b1,a1) Mean wave direction. 2 FPdir = atan2(b2,a2)/2 Principal wave direction. 3 FSpr = sqrt(2*(1-C1)) Directional Spread of Mdir 4 FSkew = -C2*sin(2*(Pdir-Mdir))/Spr^3 Circular Skewness of Mdir 5 FMSpr = atan2(sqrt((0.5*b1.^2*(1+a2))-(a1*b1*b2)+... (0.5*a1.^2*(1-a2))),C1^2) Mean spreading angle 6 FLcrst = sqrt((1-C2)/(1+C2)) Long-Crestedness parameter 7 FS1 = C1/(1-C1) Cos^{2S} distribution dispersion parameter, S 8 FS2 = [1+3*C2+sqrt(1+(14+C2)*C2)]/(2*(1-C2)) Alternative estimate of S 9 FD1 = sqrt(-2*log(C1)) Wrapped Normal distribution parameter, D. 10 FD2 = sqrt(-log(C2)/2) Alternative estimate of D, see spreading.m. Factors calculated at the peak frequency, fp = 1/Tp: 11 TpMdir = Mean wave direction at the spectral peak 12 TpSpr = Directional Spread of TpMdir 13 TpSkew = Skewness of TpMdir 14 Wdir = {theta(i) | [y,i]=max(max(S.S,[],2))} Main wave direction Factors calculated from spectrally weighted averages of the Fourier coefficients a and b (frequency independent parameters): 15 Wdir2 = {theta(i) | D(theta(i))==max(D(theta))} Main wave direction 16 Mdir = atan2(b1,a1) Mean wave direction. 17 Pdir = atan2(b2,a2)/2 Principal wave direction. 18 Spr = sqrt(2*(1-C1)) Directional Spread of Mdir 19 Skew = -C2*sin(2*(Pdir-Mdir))/Spr^3 Circular Skewness of Mdir 20 MSpr = atan2(sqrt((0.5*b1.^2*(1+a2))-(a1*b1*b2)+... (0.5*a1.^2*(1-a2))),C1^2) Mean spreading angle 21 Lcrst = sqrt((1-C2)/(1-C2)) Long-Crestedness parameter 22 S1 = C1/(1-C1) Cos^{2s} distribution dispersion parameter, S 23 S2 = [1+3*C2+sqrt(1+(14+C2)*C2)]/(2*(1-C2)) Alternative estimate of S 24 D1 = sqrt(-2*log(C1)) Wrapped Normal distribution parameter, D. 25 D2 = sqrt(-log(C2)/2) Alternative estimate of D, see spreading.m. 26 TMdir = atan2(b,a) Mean wave direction Tucker's method. Note: All angles are given in radians. Examples: S = demospec('dir'); [ch,txt] = dspec2char(S,1:26); % fact a vector of integers ch0 = cell2struct(ch,txt,2); % Make a structure [txt{2,:}]=deal(','); txt{2,end}=' '; eval(['[' txt{:} '] = deal(ch{:})']) % Assign values to variables plot(S.w,FMdir) ch1 = dspec2char(S,'wdir'); % fact a string ch2 = dspec2char(S,{'mdir','pdir'}); % fact a cellarray of strings ch3 = dspec2char(S,'mdir','pdir'); % strings See also spec2char, spec2bw, spec2mom, spreading
Evaluates directional spectral characteristics | |
Returns Fourier coefficients. | |
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. | |
Display help text in Command Window. | |
Hold current graph. | |
Inverse discrete Fourier transform. | |
True for cell array. | |
True for character array (string). | |
True if field is in structure array. | |
Display legend. | |
Convert string to lowercase. | |
Wait for user response. | |
Linear plot. | |
Find possible matches for string. | |
Create axes in tiled positions. |
Evaluates directional spectral characteristics |
0001 function [ch,chtext] = dspec2chardspec2char(S,varargin) 0002 %DSPEC2CHAR Evaluates directional spectral characteristics 0003 % 0004 % CALL: [ch,chtext] = dspec2char(S,fact); 0005 % 0006 % ch = a cell vector of spectral characteristics 0007 % chtext = a cell vector of strings describing the elements of ch, see example. 0008 % S = Directional spectral density struct with angular frequency 0009 % fact = vector of factor integers or a string or 0010 % a cellarray of strings, see below.(default [1]) 0011 % 0012 % DSPEC2CHAR assumes S is a directional spectrum S(w,theta)= 0013 % S(w)*D(theta,w). If input spectrum is of wave number type, it is 0014 % transformed into a directional spectrum before the calculations. 0015 % For many of the parameters the Fourier series expansion of D(theta,w) is used: 0016 % M-1 0017 % D(theta(i)) = {1 + 2*sum [an*cos(n*theta(i)) + bn*sin(n*theta(i))]}/(2*pi) 0018 % n=1 0019 % where C1 = sqrt(a1^2+b1^2) and C2 = sqrt(a2^2+b2^2) 0020 % 0021 % Order of output is same as order in 'factors'. 0022 % Input vector 'factors' correspondence: 0023 % 0024 % Factors calculated at every frequency (frequency dependent parameters): 0025 % 1 FMdir = atan2(b1,a1) Mean wave direction. 0026 % 2 FPdir = atan2(b2,a2)/2 Principal wave direction. 0027 % 3 FSpr = sqrt(2*(1-C1)) Directional Spread of Mdir 0028 % 4 FSkew = -C2*sin(2*(Pdir-Mdir))/Spr^3 Circular Skewness of Mdir 0029 % 5 FMSpr = atan2(sqrt((0.5*b1.^2*(1+a2))-(a1*b1*b2)+... 0030 % (0.5*a1.^2*(1-a2))),C1^2) Mean spreading angle 0031 % 6 FLcrst = sqrt((1-C2)/(1+C2)) Long-Crestedness parameter 0032 % 7 FS1 = C1/(1-C1) Cos^{2S} distribution dispersion parameter, S 0033 % 8 FS2 = [1+3*C2+sqrt(1+(14+C2)*C2)]/(2*(1-C2)) Alternative estimate of S 0034 % 9 FD1 = sqrt(-2*log(C1)) Wrapped Normal distribution parameter, D. 0035 % 10 FD2 = sqrt(-log(C2)/2) Alternative estimate of D, see spreading.m. 0036 % 0037 % Factors calculated at the peak frequency, fp = 1/Tp: 0038 % 11 TpMdir = Mean wave direction at the spectral peak 0039 % 12 TpSpr = Directional Spread of TpMdir 0040 % 13 TpSkew = Skewness of TpMdir 0041 % 14 Wdir = {theta(i) | [y,i]=max(max(S.S,[],2))} Main wave direction 0042 % 0043 % Factors calculated from spectrally weighted averages of the 0044 % Fourier coefficients a and b (frequency independent parameters): 0045 % 15 Wdir2 = {theta(i) | D(theta(i))==max(D(theta))} Main wave direction 0046 % 16 Mdir = atan2(b1,a1) Mean wave direction. 0047 % 17 Pdir = atan2(b2,a2)/2 Principal wave direction. 0048 % 18 Spr = sqrt(2*(1-C1)) Directional Spread of Mdir 0049 % 19 Skew = -C2*sin(2*(Pdir-Mdir))/Spr^3 Circular Skewness of Mdir 0050 % 20 MSpr = atan2(sqrt((0.5*b1.^2*(1+a2))-(a1*b1*b2)+... 0051 % (0.5*a1.^2*(1-a2))),C1^2) Mean spreading angle 0052 % 21 Lcrst = sqrt((1-C2)/(1-C2)) Long-Crestedness parameter 0053 % 22 S1 = C1/(1-C1) Cos^{2s} distribution dispersion parameter, S 0054 % 23 S2 = [1+3*C2+sqrt(1+(14+C2)*C2)]/(2*(1-C2)) Alternative estimate of S 0055 % 24 D1 = sqrt(-2*log(C1)) Wrapped Normal distribution parameter, D. 0056 % 25 D2 = sqrt(-log(C2)/2) Alternative estimate of D, see spreading.m. 0057 % 26 TMdir = atan2(b,a) Mean wave direction Tucker's method. 0058 % 0059 % Note: All angles are given in radians. 0060 % 0061 % Examples: 0062 % S = demospec('dir'); 0063 % [ch,txt] = dspec2char(S,1:26); % fact a vector of integers 0064 % ch0 = cell2struct(ch,txt,2); % Make a structure 0065 % [txt{2,:}]=deal(','); txt{2,end}=' '; 0066 % eval(['[' txt{:} '] = deal(ch{:})']) % Assign values to variables 0067 % plot(S.w,FMdir) 0068 % ch1 = dspec2char(S,'wdir'); % fact a string 0069 % ch2 = dspec2char(S,{'mdir','pdir'}); % fact a cellarray of strings 0070 % ch3 = dspec2char(S,'mdir','pdir'); % strings 0071 % 0072 % See also spec2char, spec2bw, spec2mom, spreading 0073 0074 % References: 0075 % Krogstad, H.E., Wolf, J., Thompson, S.P., and Wyatt, L.R. (1999) 0076 % 'Methods for intercomparison of wave measurements' 0077 % Coastal Enginering, Vol. 37, pp. 235--257 0078 % 0079 % Krogstad, H.E. (1982) 0080 % 'On the covariance of the periodogram' 0081 % Journal of time series analysis, Vol. 3, No. 3, pp. 195--207 0082 % 0083 % Tucker, M.J. (1993) 0084 % 'Recommended standard for wave data sampling and near-real-time processing' 0085 % Ocean Engineering, Vol.20, No.5, pp. 459--474 0086 % 0087 % Young, I.R. (1999) 0088 % "Wind generated ocean waves" 0089 % Elsevier Ocean Engineering Book Series, Vol. 2, pp 239 0090 % 0091 % Benoit, M. and Goasguen, G. (1999) 0092 % "Comparative evalutation of directional wave analysis techniques applied 0093 % to field measurements", In Proc. 9'th ISOPE conference, Vol III, pp 87-94. 0094 0095 % Tested on: Matlab 5.2 0096 0097 % History: 0098 %revised pab 29.06.2001 0099 % - factors 15:25 are now calculated from the scaled Fourier coefficients of 0100 % D(theta) = int S(w,theta) dw instead of D(theta) = int D(w,theta) dw 0101 %revised pab 22.06.2001 0102 % - moved code from wspecplottest into dspec2char 0103 % - Fixed bugs: The parameters are now calculated from the correct 0104 % Scaled Fourier coefficients. S.phi is now taken into account. 0105 % - Added frequency dependent parameters, skewness and kurtosis 0106 % parameters and dispersion parameter for the wrapped normal spreading function. 0107 %by Vengatesan Venugopal [V.Venugopal@hw.ac.uk] 19.06.2001 0108 0109 % Options not implemented 0110 % FKurt = 2*(C2*cos(2*(Pdir-Mdir))-4*C1+3)/Spr^4 Circular kurtosis of FMdir 0111 % TpKurt = Kurtosis of TpMdir 0112 % Kurt = 2*(C2*cos(2*(Pdir-Mdir))-4*C1+3)/Spr^4 Circular kurtosis of Mdir 0113 % TSpr = sqrt(2*(1-C1)) Directional Spread of TMdir (Wrong?) 0114 0115 transform2degrees = 0; 0116 0117 switch nargin 0118 case 0, help dspec2char, return 0119 case 1, nfact = 1; 0120 case 2, 0121 fact = varargin{1}; 0122 if ischar(fact), fact = {fact}; end 0123 otherwise 0124 fact = varargin; 0125 end 0126 0127 tfact ={'FMdir','FPdir','FSpr','FSkew','FMSpr','FLcrst','FS1','FS2','FD1','FD2',... 0128 'TpMdir','TpSpr','TpSkew','Wdir','Wdir2', ... 0129 'Mdir','Pdir','Spr','Skew','MSpr','Lcrst','S1','S2','D1','D2','TMdir','TSpr'}; 0130 0131 0132 0133 if iscell(fact) 0134 N = length(fact(:)); 0135 nfact = zeros(1,N); 0136 ltfact = char(lower(tfact)); 0137 for ix=1:N, 0138 ind = strmatch(lower(fact(ix)),ltfact,'exact'); 0139 if length(ind)==1, 0140 nfact(ix)=ind; 0141 else 0142 error(['Not a valid factor: ' fact{ix}]) 0143 end 0144 end 0145 else 0146 nfact = fact; 0147 end 0148 if any(nfact>26 | nfact<1) 0149 error('Factor outside range (1,...,26)') 0150 end 0151 0152 vari = 'w'; 0153 if isfield(S,'k2d') 0154 S = spec2spec(S,'dir'); 0155 elseif isfield(S,'theta') 0156 S = ttspec(S,'w','r'); % Make sure it is directional spectrum given in radians 0157 else 0158 error('Directional spectra required!') 0159 end 0160 0161 phi = 0; 0162 if isfield(S,'phi') & ~isempty(S.phi) 0163 phi = S.phi; 0164 end 0165 0166 0167 0168 w = getfield(S,vari); 0169 w = w(:); 0170 theta = S.theta(:)-phi; 0171 Dtf = S.S; 0172 [Nt Nf] = size(Dtf); 0173 0174 if length(w)~=Nf, 0175 error('Length of frequency vector S.f or S.w must equal size(S.S,2)'), 0176 end 0177 if length(theta)~=Nt, 0178 error('Length of angular vector S.theta must equal size(S.S,1)'), 0179 end 0180 0181 Sf = simpson(S.theta,Dtf,1); 0182 ind = find(Sf); 0183 0184 %Directional distribution D(theta,w) = S(theta,w)/S(w) 0185 Dtf(:,ind) = Dtf(:,ind)./Sf(ones(Nt,1),ind); 0186 0187 if 0, 0188 Dtheta = simpson(w,Dtf,2); %Directional spreading, D(theta) = int D(w,theta) dw 0189 else 0190 Dtheta = simpson(w,S.S,2); %Directional spreading, D(theta) = int S(w,theta) dw 0191 end 0192 Dtheta = Dtheta/simpson(S.theta,Dtheta); 0193 0194 [y,ind] = max(Dtheta); 0195 Wdir = theta(ind); % main wave direction 0196 0197 0198 0199 % Find Fourier Coefficients of Dtheta and Dtf 0200 M = 3; % No of harmonics-1 0201 [a,b] = fourier(theta,Dtheta,2*pi,10); 0202 [aa,bb] = fourier(theta,Dtf,2*pi,10); 0203 if 0, % Alternatively 0204 Fcof = 2*ifft(Dtheta); 0205 Pcor = [1; exp(sqrt(-1)*[1:M-1].'*theta(1))]; % correction term to get 0206 % the correct integration limits 0207 Fcof = Fcof(1:M,:).*Pcor; 0208 a = real(Fcof(1:M)); 0209 b = imag(Fcof(1:M)); 0210 Fcof = 2*ifft(Dtf); 0211 Fcof = Fcof(1:M,:).*Pcor(:,ones(1,Nf)); 0212 aa = real(Fcof(1:M,:)); 0213 bb = imag(Fcof(1:M,:)); 0214 end 0215 0216 if 0, 0217 % Checking the components 0218 P1toM = zeros(Nt,M); 0219 P1toM(:,1) = a(1)/2; % mean level 0220 0221 for ix=2:M, 0222 % M-1 1st harmonic 0223 P1toM(:,ix) = a(ix)*cos((ix-1)*theta)+b(ix)*sin((ix-1)*theta); 0224 end 0225 Psum = sum(P1toM,2); 0226 subplot(2,1,1) 0227 plot(theta,Dtheta-Psum,'o','MarkerSize',2), legend('Dteta-Psum') 0228 subplot(2,1,2) 0229 plot(theta,Dtheta,'o','MarkerSize',2), hold on 0230 plot(theta,Psum,'m'), hold off,legend('Dteta','Psum') 0231 pause 0232 end 0233 0234 % The parameters below are calculated for 0235 a = pi*a; b = pi*b; 0236 aa = pi*aa; bb = pi*bb; 0237 0238 0239 %Fourier coefficients for D(theta) 0240 a0=a(1); a1=a(2);a2=a(3); 0241 b1=b(2); b2=b(3); 0242 0243 %Fourier coefficients for D(theta,w) 0244 aa0 = aa(1,:); aa1 = a(2,:);aa2 = aa(3,:); 0245 bb1 = bb(2,:); bb2 = bb(3,:); 0246 0247 FC1 = sqrt(aa1.^2+bb1.^2); 0248 FC2 = sqrt(aa2.^2+bb2.^2); 0249 0250 %plot(w,FC1,w,FC2),legend('FC1','FC2') 0251 0252 0253 FMdir = atan2(bb1,aa1); % Mean wave direction 0254 FPdir = 0.5*atan2(bb2,aa2); % Principal wave direction 0255 FSpr = sqrt(2*abs(1-FC1)); % Directional spread 0256 FSkew = -FC2.*sin(2*(FPdir-FMdir))./(FSpr.^3); % Skewness of Mdir 0257 FKurt = 2*(FC2.*cos(2*(FPdir-FMdir))-4.*FC1+3)./(FSpr.^4); % Kurtosis of Mdir 0258 0259 0260 % Mean spreading angle 0261 FMSpr = atan2(sqrt(0.5*bb1.^2.*(1+aa2)-aa1.*bb1.*bb2+0.5.*aa1.^2.*(1-aa2)),FC1.^2); 0262 0263 % Long-Crestedness parameter 0264 FLcrst = sqrt((1-FC2)./(1+FC2)); 0265 0266 %Estimates of Cos^{2s} distribution dispersion parameter, S 0267 FS1 = FC1./(1-FC1); 0268 FS2 = (1+3*FC2+sqrt(1+(14+FC2).*FC2))./(2*(1-FC2)); 0269 %Estimates of Wrapped Normal distribution parameter, D. 0270 FD1 = repmat(-inf,size(FC1)); 0271 ind = find(FC1); % avoid log(0) 0272 FD1(ind)= sqrt(-2*log(FC1(ind))); 0273 0274 FD2 = repmat(-inf,size(FC2)); 0275 ind = find(FC2); % avoid log(0) 0276 FD2(ind)= sqrt(-log(FC2(ind))/2); 0277 0278 [y,ind] = max(max(S.S,[],1)); % Index to spectral peak. 0279 0280 TpMdir = FMdir(ind); % Mean wave direction at the spectral peak 0281 TpSpr = FSpr(ind); % Directional Spread of TpMdir 0282 TpSkew = FSkew(ind); % Skewness of TpMdir 0283 TpKurt = FKurt(ind); % Kurtosis of TpMdir 0284 0285 0286 [y,ind] = max(max(S.S,[],2)); % Index to spectral peak. 0287 Wdir = mod(theta(ind)+pi,2*pi)-pi; % main wave direction 0288 [y,ind] = max(Dtheta); 0289 Wdir2 = mod(theta(ind)+pi,2*pi)-pi; % main wave direction 0290 0291 0292 0293 C1 = sqrt(a1.^2+b1.^2); 0294 C2 = sqrt(a2.^2+b2.^2); 0295 0296 0297 0298 Mdir = atan2(b1,a1); % Mean wave direction 0299 Pdir = 0.5*atan2(b2,a2); % Principal wave direction 0300 Spr = sqrt(2*abs(1-C1)); % Directional spread 0301 Skew = -C2.*sin(2*(Pdir-Mdir))./(Spr.^3); % Skewness of Mdir 0302 Kurt = 2*(C2.*cos(2*(Pdir-Mdir))-4.*C1+3)./(Spr.^4); % Kurtosis of Mdir 0303 0304 % Mean spreading angle 0305 MSpr = atan2(sqrt(0.5*b1.^2.*(1+a2)-a1.*b1.*b2+0.5.*a1.^2.*(1-a2)),C1.^2); 0306 0307 % Long-Crestedness parameter 0308 Lcrst = sqrt((1-C2)./(1+C2)); 0309 0310 %Estimates of Cos^{2s} distribution dispersion parameter, S 0311 S1 = C1./(1-C1); 0312 S2 = (1+3*C2+sqrt(1+(14+C2).*C2))./(2*(1-C2)); 0313 %Estimates of Wrapped Normal distribution parameter, D. 0314 D1 = sqrt(abs(2*log(C1))); 0315 D2 = sqrt(abs(log(C2)/2)); 0316 0317 0318 % --------------------------------------------------------------------------- 0319 % TUCKER PROCEDURE for MDIR and UI (pp202) 0320 0321 0322 [Sff,ind] = max(S.S(:,:)); 0323 Wwdir = theta(ind); % main wave direction 0324 0325 %thetam = (ind-1).'*2*pi/(Nt-1)-pi; 0326 thetam = theta(ind); 0327 Sff = Sff'; 0328 bot = sum(Sff); 0329 a = sum(Sff.*cos(thetam))/bot; 0330 b = sum(Sff.*sin(thetam))/bot; 0331 0332 TMdir = atan2(b,a); 0333 UI = sqrt(a.^2+b.^2); 0334 TSpr = sqrt(2*(1-UI)); %This is not correct 0335 0336 ch ={FMdir,FPdir,FSpr,FSkew,FMSpr,FLcrst,FS1,FS2,FD1,FD2,... 0337 TpMdir,TpSpr,TpSkew,Wdir,Wdir2,Mdir,Pdir,Spr,Skew,.... 0338 MSpr,Lcrst,S1,S2,D1,D2,TMdir}; 0339 0340 % Make sure the angles are between -pi<= theta < pi 0341 ind = [1:3 5 11:12 14:18 20 26]; 0342 for ix = ind, 0343 ch{ix} = mod(ch{ix}+pi,2*pi)-pi; 0344 end 0345 if transform2degrees, % Change from radians to degrees 0346 for ix = ind, 0347 ch{ix} = ch{ix}*180/pi; 0348 end 0349 % Old call: 0350 %r2d = ones(1,11); 0351 %r2d([1:3 5 7 ]) = 180/pi; 0352 %ch = ch.*r2d; 0353 end 0354 0355 % select the appropriate values 0356 ch = ch(nfact); 0357 chtext = tfact(nfact); 0358 0359 return 0360 0361
Comments or corrections to the WAFO group