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
Estimate one-sided spectral density from data. | |
returns the frequency type of a Spectral density struct. | |
returns the constant acceleration of gravity | |
Evaluates spectral characteristics and their covariance | |
Simulates a Randomized 2nd order non-linear wave X(t) | |
Interpolation and zero-padding of spectrum | |
Toggle Transform between angular frequency and frequency spectrum | |
Translates from frequency to wave number | |
Concatenate arrays. | |
String representation of date. | |
Flush pending graphics events. | |
Display message and abort function. | |
Create figure window. | |
Get handle to current axis. | |
Get structure field contents. | |
Hold current graph. | |
Quick 1-D linear interpolation. | |
Current date and time as date number. | |
Convert number to string. (Fast version) | |
Wait for user response. | |
Linear plot. | |
Set object properties. | |
Tile figure windows. | |
Graph title. |
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
Comments or corrections to the WAFO group