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

## CROSS-REFERENCE INFORMATION

This function calls:
 dat2spec Estimate one-sided spectral density from data. freqtype returns the frequency type of a Spectral density struct. gravity returns the constant acceleration of gravity spec2char Evaluates spectral characteristics and their covariance spec2nlsdat Simulates a Randomized 2nd order non-linear wave X(t) specinterp Interpolation and zero-padding of spectrum ttspec Toggle Transform between angular frequency and frequency spectrum w2k Translates from frequency to wave number cat Concatenate arrays. datestr String representation of date. drawnow Flush pending graphics events. error Display message and abort function. figure Create figure window. gca Get handle to current axis. getfield Get structure field contents. hold Hold current graph. interp1q Quick 1-D linear interpolation. now Current date and time as date number. num2str Convert number to string. (Fast version) pause Wait for user response. plot Linear plot. set Set object properties. tilefigs Tile figure windows. title Graph title.
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 %
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