Home > wafo > wsim > spec2linspec.m

spec2linspec

PURPOSE ^

Separates the linear component of the Spectrum

SYNOPSIS ^

[SL,SN]=spec2linspec(S,np,dt,iseed,fnLimit)

DESCRIPTION ^

 SPEC2LINSPEC  Separates the linear component of the Spectrum  
               according to 2nd order wave theory 
             
    CALL: [SL,SN] = spec2linspec(S,[np cases],dt,iseed,fnLimit); 
   
    SL    = spectral density structure with linear components only. 
    SN    = spectral density structure with non-linear components only. 
    S     = spectral density structure with linear and non-linear. 
            components 
    np    = giving np load points.  (default length(S)-1=n-1). 
            If np>n-1 it is assummed that S(k)=0 for all k>n-1 
    cases = number of cases (default=20)  
    dt    = step in grid (default dt is defined by the Nyquist freq) 
    iseed = starting seed number for the random number generator  
           (default none is set) 
  fnLimit = normalized upper frequency limit of spectrum for 2'nd order 
            components. The frequency is normalized with  
            sqrt(gravity*tanh(kbar*waterDepth)/Amax)/(2*pi) 
            (default sqrt(2), i.e., Convergence criterion). 
            Generally this should be the same as used in the final 
            non-linear simulation (see example below). 
  
   SPEC2LINSPEC separates the linear and non-linear component of the Spectrum  
   according to 2nd order wave theory. This is useful when simulating 
   non-linear waves because: 
   If the spectrum does not decay rapidly enough towards zero, the 
   contribution from the 2nd order wave components at the upper tail can 
   be very large and unphysical. 
   Another option to ensure convergence of the perturbation series in the 
   simulation, is to truncate the upper tail of the 
   spectrum at FNLIMIT in the calculation of the 2nd order 
   wave components, i.e., in the calculation of sum and difference 
   frequency effects.  
  
  Example: 
   np = 10000; 
   iseed = 1; 
   pflag = 2; 
   S  = jonswap(10); 
   fnLimit = inf;   
   [SL,SN] = spec2linspec(S,np,[],[],fnLimit); 
   x0 = spec2nlsdat(SL,8*np,[],iseed,[],fnLimit); 
   x1 = spec2nlsdat(S,8*np,[],iseed,[],fnLimit);  
   x2 = spec2nlsdat(S,8*np,[],iseed,[],sqrt(2));   
   Se0 = dat2spec(x0); 
   Se1 = dat2spec(x1); 
   Se2 = dat2spec(x2);  
   clf   
   wspecplot(SL,'r',pflag),  % Linear components 
    hold on 
   wspecplot(S,'b',pflag)    % target spectrum for simulated data 
   wspecplot(Se0,'m',pflag), % approx. same as S  
   wspecplot(Se1,'g',pflag)  % unphysical spectrum 
   wspecplot(Se2,'k',pflag)  % approx. same as S 
   axis([0 10 -80 0]) 
   hold off 
    
  See also  spec2nlsdat

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [SL,SN]=spec2linspec(S,np,dt,iseed,fnLimit) 
002 %SPEC2LINSPEC  Separates the linear component of the Spectrum  
003 %              according to 2nd order wave theory 
004 %            
005 %   CALL: [SL,SN] = spec2linspec(S,[np cases],dt,iseed,fnLimit); 
006 %  
007 %   SL    = spectral density structure with linear components only. 
008 %   SN    = spectral density structure with non-linear components only. 
009 %   S     = spectral density structure with linear and non-linear. 
010 %           components 
011 %   np    = giving np load points.  (default length(S)-1=n-1). 
012 %           If np>n-1 it is assummed that S(k)=0 for all k>n-1 
013 %   cases = number of cases (default=20)  
014 %   dt    = step in grid (default dt is defined by the Nyquist freq) 
015 %   iseed = starting seed number for the random number generator  
016 %          (default none is set) 
017 % fnLimit = normalized upper frequency limit of spectrum for 2'nd order 
018 %           components. The frequency is normalized with  
019 %           sqrt(gravity*tanh(kbar*waterDepth)/Amax)/(2*pi) 
020 %           (default sqrt(2), i.e., Convergence criterion). 
021 %           Generally this should be the same as used in the final 
022 %           non-linear simulation (see example below). 
023 % 
024 %  SPEC2LINSPEC separates the linear and non-linear component of the Spectrum  
025 %  according to 2nd order wave theory. This is useful when simulating 
026 %  non-linear waves because: 
027 %  If the spectrum does not decay rapidly enough towards zero, the 
028 %  contribution from the 2nd order wave components at the upper tail can 
029 %  be very large and unphysical. 
030 %  Another option to ensure convergence of the perturbation series in the 
031 %  simulation, is to truncate the upper tail of the 
032 %  spectrum at FNLIMIT in the calculation of the 2nd order 
033 %  wave components, i.e., in the calculation of sum and difference 
034 %  frequency effects.  
035 % 
036 % Example: 
037 %  np = 10000; 
038 %  iseed = 1; 
039 %  pflag = 2; 
040 %  S  = jonswap(10); 
041 %  fnLimit = inf;   
042 %  [SL,SN] = spec2linspec(S,np,[],[],fnLimit); 
043 %  x0 = spec2nlsdat(SL,8*np,[],iseed,[],fnLimit); 
044 %  x1 = spec2nlsdat(S,8*np,[],iseed,[],fnLimit);  
045 %  x2 = spec2nlsdat(S,8*np,[],iseed,[],sqrt(2));   
046 %  Se0 = dat2spec(x0); 
047 %  Se1 = dat2spec(x1); 
048 %  Se2 = dat2spec(x2);  
049 %  clf   
050 %  wspecplot(SL,'r',pflag),  % Linear components 
051 %   hold on 
052 %  wspecplot(S,'b',pflag)    % target spectrum for simulated data 
053 %  wspecplot(Se0,'m',pflag), % approx. same as S  
054 %  wspecplot(Se1,'g',pflag)  % unphysical spectrum 
055 %  wspecplot(Se2,'k',pflag)  % approx. same as S 
056 %  axis([0 10 -80 0]) 
057 %  hold off 
058 %   
059 % See also  spec2nlsdat 
060  
061 % Reference  
062 % P. A. Brodtkorb (2004),  
063 % The probability of Occurrence of dangerous Wave Situations at Sea. 
064 %  Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
065 %  Trondheim, Norway. 
066 %   
067 % Nestegaard, A  and Stokka T (1995) 
068 % A Third Order Random Wave model. 
069 % In proc.ISOPE conf., Vol III, pp 136-142. 
070 % 
071 % R. S Langley (1987) 
072 % A statistical analysis of non-linear random waves. 
073 % Ocean Engng, Vol 14, pp 389-407 
074 % 
075 % Marthinsen, T. and Winterstein, S.R (1992) 
076 % 'On the skewness of random surface waves' 
077 % In proc. ISOPE Conf., San Francisco, 14-19 june. 
078  
079 % tested on: Matlab 5.3 
080 % History: 
081 % Revised pab Nov 2004   
082 % changed the default constant controlling its 
083 %                 performance. Can be improved further 
084 % by pab 13.08.2002 
085  
086 % TODO % Replace inputs with options structure 
087 % TODO % Can be improved further. 
088 error(nargchk(1,5,nargin)) 
089  
090 % Define some constants 
091 %fnLimit = sqrt(inf) 
092 method = 'apstochastic'; 
093 trace     = 1; % trace the convergence 
094 maxSim    = 30; 
095 tolerance = 5e-5; 
096 cases     = 20; 
097 L         = 200; %maximum lag size of the window function used in 
098                  %spectral estimate 
099 ftype = freqtype(S); %options are 'f' and 'w' and 'k' 
100 n     = length(getfield(S,ftype)); 
101 switch ftype 
102  case 'f',  
103   ftype = 'w'; 
104   S = ttspec(S,ftype); 
105 end 
106 Hm0  = spec2char(S,'Hm0'); 
107 Tm02 = spec2char(S,'Tm02'); 
108  
109 if (nargin<5 | isempty(fnLimit)) 
110   fnLimit = sqrt(2); 
111 end 
112 if (nargin>3 & ~isempty(iseed)),  
113   randn('state',iseed); % set the the seed  
114 else 
115   iseed = 0; 
116 end    
117 if (nargin<2 | isempty(np)),     np = max(n-1,5000);  end 
118 if (nargin>2 & ~isempty(dt)),    S = specinterp(S,dt);end  % interpolate spectrum   
119                                
120 switch  length(np)  
121   case 1,  
122   case 2, cases=np(2); np=np(1); 
123   otherwise, error('Wrong input. Too many arguments') 
124 end 
125 np = np + mod(np,2); % make sure np is even     
126  
127  
128  
129 waterDepth = abs(S.h); 
130 kbar = w2k(2*pi/Tm02,0,waterDepth); 
131  
132 % Expected maximum amplitude for 1000 waves seastate 
133 numWaves = 10000;   
134 Amax = sqrt(2*log(numWaves))*Hm0/4;  
135    
136 fLimitLo = sqrt(gravity*tanh(kbar*waterDepth)*Amax/waterDepth^3); 
137  
138  
139 freq = getfield(S,ftype); 
140 freq(end) = freq(end)-sqrt(eps); 
141 Hw2  = 0; 
142  
143 SL = S; 
144  
145 indZero = find(freq<fLimitLo); 
146 if any(indZero) 
147   SL.S(indZero) = 0; 
148 end 
149 maxS = max(S.S); 
150 Fs = 2*freq(end)+eps; % sampling frequency 
151  
152 for ix=1:maxSim 
153   [x2,x1] = spec2nlsdat(SL,[np,cases],[],iseed,method,fnLimit); 
154   %x2(:,2:end) = x2(:,2:end) -x1(:,2:end); 
155   S2 = dat2spec(x2,L);    
156   S1 = dat2spec(x1,L); 
157   %[tf21,fi] = tfe(x2(:,2),x1(:,2),1024,Fs,[],512); 
158   %Hw11 = interp1q(fi,tf21.*conj(tf21),freq); 
159    
160   Hw1  = interp1q( S2.w,abs(S1.S./S2.S),freq); 
161   %Hw1  = (interp1q( S2.w,abs(S1.S./S2.S),freq)+Hw2)/2; 
162   %plot(freq, abs(Hw11-Hw1),'g') 
163   %title('diff') 
164   %pause 
165   %clf 
166    
167   %d1 = interp1q( S2.w,S2.S,freq);; 
168    
169   SL.S = (Hw1.*S.S); 
170    
171   if any(indZero) 
172     SL.S(indZero) = 0; 
173   end 
174   k = find(SL.S< 0); 
175   if any(k), % Make sure that the current guess is larger than zero 
176     k 
177     Hw1(k) 
178     Hw1(k)  = min((S1.S(k)*0.9),(S.S(k))); 
179     SL.S(k) = max(Hw1(k).*S.S(k),eps); 
180   end 
181   Hw12 = Hw1-Hw2; 
182   maxHw12 = max(abs(Hw12)); 
183   if trace==1, 
184     figure(1), plot(freq,Hw1,'r'), hold on, title('Hw') 
185     set(gca,'yscale','log') 
186     figure(2), plot(freq,abs(Hw12),'r'), hold on, title('Hw-HwOld') 
187     set(gca,'yscale','log') 
188     drawnow 
189     pause(3) 
190     figure(1), plot(freq,Hw1,'b'), hold on, title('Hw') 
191     figure(2), plot(freq,abs(Hw12),'b'), hold on, title('Hw-HwOld') 
192     tilefigs 
193   end 
194     
195    disp(['Iteration ',num2str(ix),... 
196         ' Hw12 :  ' num2str(maxHw12), ... 
197     ' Hw12/maxS : ' num2str(maxHw12/maxS)]), 
198   if (maxHw12<maxS*tolerance & Hw1(end)<Hw2(end) ) 
199     break 
200   end 
201   Hw2 = Hw1; 
202 end 
203  
204 %Hw1(end) 
205 %maxS*1e-3 
206 %if Hw1(end)*S.>maxS*1e-3, 
207 %  warning('The Nyquist frequency of the spectrum may be too low') 
208 %end 
209  
210 SL.date = datestr(now); 
211 if nargout>1 
212   SN   = SL; 
213   SN.S = S.S-SL.S; 
214   SN.note = cat(2,SN.note,' non-linear component (spec2linspec)'); 
215 end 
216 SL.note = cat(2,SL.note,' linear component (spec2linspec)'); 
217  
218 return 
219  
220  
221  
222

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