Home > wafo > spec > wnormspec.m

wnormspec

PURPOSE ^

Normalize a spectral density such that m0=m2=1

SYNOPSIS ^

[S,mn4,m0,m2,m4,m1]=wnormspec(spectrum,m0,m2,plotflag)

DESCRIPTION ^

 WNORMSPEC Normalize a spectral density such that m0=m2=1
 
  CALL:  [Sn,mn4,m0,m2,m4,m1]=wnormspec(S,m0,m2,plotflag);
 
         Sn  = the normalized spectrum struct
         mn4 = the fourth spectral moment of Sn.
         mi  = the i'th spectral moment of S.
         S   = the original spectrum struct
    plotflag = 0, do not plot the normalized spectrum (default).
               1, plot the normalized spectrum.
            
  Normalization performed such that
     INT S(freq) dfreq = 1       INT freq^2  S(freq) dfreq = 1
  where integration limits are given by  freq  and  S(freq)  is the 
  spectral density; freq can be frequency or wave number.
  The normalization is defined by
  A=sqrt(m0/m2); B=1/A/m0; freq'=freq*A; S(freq')=S(freq)*B;
 
  If S is a directional spectrum then a normalized gravity (.g) is added
  to Sn, such that mxx normalizes to 1, as well as m0 and mtt.
  (See spec2mom for notation of moments)
 
  If S is complex-valued cross spectral density which has to be
  normalized, then m0, m2 (suitable spectral moments) should be given.
 
  Example: 
    S = jonswap;
    [Sn,mn4] = wnormspec(S);
    mts = spec2mom(S,2)     % Should be equal to one!
  
  See also  spec2mom, spec2spec

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [S,mn4,m0,m2,m4,m1]=wnormspec(spectrum,m0,m2,plotflag)
002 %WNORMSPEC Normalize a spectral density such that m0=m2=1
003 %
004 % CALL:  [Sn,mn4,m0,m2,m4,m1]=wnormspec(S,m0,m2,plotflag);
005 %
006 %        Sn  = the normalized spectrum struct
007 %        mn4 = the fourth spectral moment of Sn.
008 %        mi  = the i'th spectral moment of S.
009 %        S   = the original spectrum struct
010 %   plotflag = 0, do not plot the normalized spectrum (default).
011 %              1, plot the normalized spectrum.
012 %           
013 % Normalization performed such that
014 %    INT S(freq) dfreq = 1       INT freq^2  S(freq) dfreq = 1
015 % where integration limits are given by  freq  and  S(freq)  is the 
016 % spectral density; freq can be frequency or wave number.
017 % The normalization is defined by
018 % A=sqrt(m0/m2); B=1/A/m0; freq'=freq*A; S(freq')=S(freq)*B;
019 %
020 % If S is a directional spectrum then a normalized gravity (.g) is added
021 % to Sn, such that mxx normalizes to 1, as well as m0 and mtt.
022 % (See spec2mom for notation of moments)
023 %
024 % If S is complex-valued cross spectral density which has to be
025 % normalized, then m0, m2 (suitable spectral moments) should be given.
026 %
027 % Example: 
028 %   S = jonswap;
029 %   [Sn,mn4] = wnormspec(S);
030 %   mts = spec2mom(S,2)     % Should be equal to one!
031 % 
032 % See also  spec2mom, spec2spec 
033 
034 % Tested on: Matlab 5.3
035 %
036 % History:
037 % revised ir 31.08.01, introducing normalizing  g for w-spectrum
038 % Revised by jr 20.04.2001
039 % - Condition on nargout related to calculation of mn4; minor change
040 % - Updated information, added example
041 % By es 28.09.1999
042 
043 if  nargin==2 & ~isempty(m0)
044   plotflag=m0;
045 end
046 if (nargin<4 & nargin ~=2)|isempty(plotflag)
047   if nargout==0,
048     plotflag=1;
049   else
050     plotflag=0;
051   end
052 end
053 
054 if strcmp(lower(spectrum.type(end-2:end)),'k2d')
055   intype=spectrum.type;
056   spectrum=spec2spec(spectrum,'dir');
057 end
058   
059 S=spectrum.S;   % size np x nf
060 
061 if isfield(spectrum,'w')
062   f=spectrum.w(:); % length nf
063 elseif isfield(spectrum,'k')
064   f=spectrum.k(:);
065 else
066   f=2*pi*spectrum.f(:);
067   S=S/2/pi;
068 end
069 
070 if (nargin<3)|isempty(m0)|isempty(m2)
071   S1=abs(S);
072   if strcmp(lower(spectrum.type(end-2:end)),'dir')
073     % integrate out theta:
074     S2=trapz(spectrum.theta(:),...
075          S1.*(cos(spectrum.theta(:)*ones(size(f'))).^2),1);
076    
077     m20=trapz(f,f.^4/gravity^2.*S2.',1); % second order moment in x
078     m02=trapz(f,f.^4/gravity^2.*trapz(spectrum.theta(:),S1.*...
079     (sin(spectrum.theta(:)*ones(size(f'))).^2),1).',1); % second order moment in y
080     S1=trapz(spectrum.theta(:),S1,1).'; % integrate out theta
081     mom=trapz(f,[S1  f.*S1 f.^2.*S1 f.^4.*S1] ,1);
082     m0=mom(1); m1=mom(2); m2=mom(3);  m4=mom(4);
083   else
084     S1=S1(:);
085     mom=trapz(f,[S1  f.*S1 f.^2.*S1 f.^4.*S1] ,1);
086     m0=mom(1); m1=mom(2); m2=mom(3);  m4=mom(4);
087     if isfield(spectrum,'w')
088        m02=m4/gravity^2;
089         m20=m02;
090     end    
091   end
092 end
093  
094 SM0=sqrt(m0); SM2=sqrt(m2);
095 A=SM0/SM2;
096 B=SM2/(SM0*m0);
097 
098 f=f*A;
099 S1=S*B;
100 
101 if (nargout >= 2)&nargin<3
102   mn4 = m4*A^5*B;
103 elseif (nargout >= 2)& (nargin>3)
104   mn4=trapz(f,[f.^4.*S1] );
105 end
106 
107 S=spectrum;
108 S.S=S1;
109 if isfield(spectrum,'w')
110   S.w=f;
111 elseif isfield(spectrum,'k')
112   S.k=f;
113 else
114   S.f=f/2/pi;
115   S.S=S.S*2*pi;
116 end
117 if strcmp(lower(spectrum.type(end-2:end)),'dir')
118   S.g=[gravity*sqrt(m0*m20)/m2 gravity*sqrt(m0*m02)/m2];
119   % normalization gravity in x and y, resp.
120 end
121 if isfield(spectrum,'w')
122   S.g=[gravity*sqrt(m0*m20)/m2];
123 end
124 
125 S.norm=1;
126 S.date=datestr(now);
127 if exist('intype')
128   S=spec2spec(S,intype);
129 end  
130 if plotflag
131   wspecplot(S)
132 end
133

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