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)+...
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)+...
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

## CROSS-REFERENCE INFORMATION

This function calls:
 dspec2char Evaluates directional spectral characteristics fourier Returns Fourier coefficients. simpson Numerical integration with the Simpson method spec2spec Transforms between different types of spectra ttspec Toggle Transform between angular frequency and frequency spectrum char Create character array (string). error Display message and abort function. getfield Get structure field contents. help Display help text in Command Window. hold Hold current graph. ifft Inverse discrete Fourier transform. iscell True for cell array. ischar True for character array (string). isfield True if field is in structure array. legend Display legend. lower Convert string to lowercase. pause Wait for user response. plot Linear plot. strmatch Find possible matches for string. subplot Create axes in tiled positions.
This function is called by:
 dspec2char Evaluates directional spectral characteristics

## 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 %
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 %
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
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
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