Home > wafo > spec > dspec2char.m

dspec2char

PURPOSE ^

Evaluates directional spectral characteristics

SYNOPSIS ^

[ch,chtext] = dspec2char(S,varargin)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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