DAT2SPEC Estimate one-sided spectral density from data. CALL: S = dat2spec(x,L,g,plotflag,p,method,dflag,ftype) S = A structure containing: S = spectral density w = angular frequency tr = transformation g h = water depth (default inf) type = 'freq' note = Memorandum string date = Date and time of creation L = maximum lag size of the window function. CI = lower and upper confidence constant p = confidence level. (Default 0.95). Bw = Bandwidth of the smoothing window which is used in the estimated spectrum. (rad/sec or Hz) x = m column data matrix with sampled times in the first column and values the next columns. L = maximum lag size of the window function. If no value is given the lag size is set to be the lag where the auto correlation is less than 2 standard deviations. (maximum 300) g = the transformation assuming that x is a sample of a transformed Gaussian process. If g is empty then x is a sample of a Gaussian process (Default) plotflag = 1 plots the spectrum, S, 2 plot 10log10(S) and 3 plots both the above plots Method = 'cov' Frequency smoothing using a parzen window function on the estimated autocovariance function. (default) 'psd' Welch's averaged periodogram method with no overlapping batches 'psdo' Welch's averaged periodogram method with overlapping batches 'pmem' Maximum Entropy Method (psd using the Yule-Walker AR method) 'pburg' Burg's method. dflag = specifies a detrending performed on the signal before estimation. 'mean','linear' or 'ma' (= moving average) (default 'mean') ftype = frequency type, 'w' or 'f' (default 'w') Method == 'cov','psd': As L decreases the estimate becomes smoother and Bw increases. If we want to resolve peaks in S which is Bf (Hz or rad/sec) apart then Bw < Bf. Method == 'pmem','pburg': L denotes the order of the AR (AutoRegressive) model. NOTE: The strings method,dflag and ftype may be given anywhere after x and in any order. See also dat2tr, dat2cov
returns the N-point Bingham window in a column vector. | |
Spectrum structure constructor | |
Estimate auto covariance function from data. | |
Estimate transformation, g, from data. | |
Removes a trend from data using a moving average | |
returns the N-point Hanning window in a column vector. | |
returns the N-point Parzen window in a column vector. | |
Transforms process X and up to four derivatives | |
Inverse of the Chi squared distribution function | |
Normalize a spectral density such that m0=m2=1 | |
Plot a spectral density | |
Remove a linear trend from a vector, usually for FFT processing. | |
Discrete Fourier transform. | |
Input argument name. | |
True for character array (string). | |
Convert string to lowercase. | |
Average or mean value. | |
Next higher power of 2. | |
Convert number to string. (Fast version) | |
Power Spectral Density (PSD) estimate via Burg's method. | |
Power Spectral Density estimate. | |
Power Spectral Density (PSD) estimate via Yule-Walker's method. | |
Standard deviation. | |
Compare strings. | |
Display warning message; disable or enable warning messages. |
% CHAPTER1 demonstrates some applications of WAFO | |
% CHAPTER2 Modelling random loads and stochastic waves | |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
Script to computer exercises 1 | |
Script to computer exercises 3 | |
setup all global variables of the RECDEMO | |
Separates the linear component of the Spectrum | |
Joint distribution (pdf) of crest front velocity and wave height: | |
Joint distribution (pdf) of crest front period, Tcf, and crest amplitude, Ac | |
setup all global variables of the WAFODEMO |
0001 function [S,fcut] = dat2spec(xn,varargin) 0002 %DAT2SPEC Estimate one-sided spectral density from data. 0003 % 0004 % CALL: S = dat2spec(x,L,g,plotflag,p,method,dflag,ftype) 0005 % 0006 % S = A structure containing: 0007 % S = spectral density 0008 % w = angular frequency 0009 % tr = transformation g 0010 % h = water depth (default inf) 0011 % type = 'freq' 0012 % note = Memorandum string 0013 % date = Date and time of creation 0014 % L = maximum lag size of the window function. 0015 % CI = lower and upper confidence constant 0016 % p = confidence level. (Default 0.95). 0017 % Bw = Bandwidth of the smoothing window which is used 0018 % in the estimated spectrum. (rad/sec or Hz) 0019 % 0020 % x = m column data matrix with sampled times in the first column 0021 % and values the next columns. 0022 % 0023 % L = maximum lag size of the window function. 0024 % If no value is given the lag size is set to 0025 % be the lag where the auto correlation is less than 0026 % 2 standard deviations. (maximum 300) 0027 % 0028 % g = the transformation assuming that x is a sample of a 0029 % transformed Gaussian process. If g is empty then 0030 % x is a sample of a Gaussian process (Default) 0031 % 0032 % plotflag = 1 plots the spectrum, S, 0033 % 2 plot 10log10(S) and 0034 % 3 plots both the above plots 0035 % 0036 % Method = 'cov' Frequency smoothing using a parzen window function 0037 % on the estimated autocovariance function. (default) 0038 % 'psd' Welch's averaged periodogram method with no overlapping 0039 % batches 0040 % 'psdo' Welch's averaged periodogram method with overlapping 0041 % batches 0042 % 'pmem' Maximum Entropy Method (psd using the Yule-Walker 0043 % AR method) 0044 % 'pburg' Burg's method. 0045 % 0046 % dflag = specifies a detrending performed on the signal before estimation. 0047 % 'mean','linear' or 'ma' (= moving average) (default 'mean') 0048 % ftype = frequency type, 'w' or 'f' (default 'w') 0049 % 0050 % Method == 'cov','psd': 0051 % As L decreases the estimate becomes smoother and Bw increases. If we 0052 % want to resolve peaks in S which is Bf (Hz or rad/sec) apart then Bw < Bf. 0053 % 0054 % Method == 'pmem','pburg': 0055 % L denotes the order of the AR (AutoRegressive) model. 0056 % 0057 % NOTE: The strings method,dflag and ftype may be given anywhere after x 0058 % and in any order. 0059 % 0060 % See also dat2tr, dat2cov 0061 0062 % Secret option: if chopOffHighFreq==1 0063 % Chop off high frequencies in order to 0064 % get the same irregularity factor in the spectrum as 0065 % in x. 0066 0067 0068 % References: 0069 % Georg Lindgren and Holger Rootzen (1986) 0070 % "Stationära stokastiska processer", pp 173--176. 0071 % 0072 % Gareth Janacek and Louise Swift (1993) 0073 % "TIME SERIES forecasting, simulation, applications", 0074 % pp 75--76 and 261--268 0075 % 0076 % Emanuel Parzen (1962), 0077 % "Stochastic Processes", HOLDEN-DAY, 0078 % pp 66--103 0079 0080 % 0081 0082 % Tested on: Matlab 5.3 0083 % history: 0084 % revised pab 03.11.2000 0085 % - changed call from chi2inv to wchi2inv 0086 % revised pab 10.07.00 0087 % - fixed a bug for m>2 and g is given: replaced dat2gaus call with tranproc 0088 % revised pab 12.06.2000 0089 % - added method,dflag,ftype to input arguments 0090 % - removed dT wdef from input arguments 0091 % revised pab 14.02.2000 0092 % - added rate for interpolation in frequency domain 0093 % revised pab 20.01.2000 0094 % - added detrending option linear, ma mean 0095 % - added tapering of data before estimation 0096 % Modified by Per A. Brodtkorb 14.08.98,25.05.98 0097 % - add a nugget effect to ensure that round off errors 0098 % do not result in negative spectral estimates 0099 % Modified by svi 29.09.99 0100 % - The program is not estimating transformation g any more. 0101 % modified pab 22.10.1999 0102 % - fixed so that x in fact can be a m column matrix 0103 % - updated info put into the spectral structure. 0104 % - updated help header 0105 % modified pab 03.11.1999 0106 % - fixed a bug: line 152 wrong array dim. when m>2 0107 0108 0109 % Initialize constants 0110 %~~~~~~~~~~~~~~~~~~~~~ 0111 nugget = 0; %10^-12; 0112 rate = 2; % interpolationrate for frequency 0113 tapery = 0; % taper the data before the analysis 0114 wdef = 1; % 1=parzen window 2=hanning window 0115 0116 % Default values: 0117 %~~~~~~~~~~~~~~~~ 0118 L = []; 0119 g = []; 0120 plotflag = 0; 0121 p = [];%0.95; 0122 chopOffHighFreq=0; % chop off high frequencies in order to get the same 0123 % irregularity factor in the spectrum as in the data 0124 % may not be a good idea => default is 0 0125 0126 method = 'cov'; % cov. other options from signal toolbox: psd = welch's method, pyulear,pmem = 0127 % maximum entropy method 0128 dflag = 'mean'; %'ma','linear' 'mean','none' % detrending option 0129 ftype = 'w' ; %options are 'f' and 'w' 0130 0131 0132 P = varargin; 0133 Np = length(P); 0134 if Np>0 0135 strix = zeros(1,Np); 0136 for ix=1:Np, % finding symbol strings 0137 strix(ix)=ischar(P{ix}); 0138 end 0139 k = find(strix); 0140 Nk = length(k); 0141 if Nk>0 0142 if Nk>3, warning('More than 3 strings are not allowed'), end 0143 for ix = k 0144 switch lower(P{ix}) 0145 case {'f','w'}, ftype = P{ix}; 0146 case {'cov','pmem','mem','psd','psdo','pburg'}, method = P{ix}; 0147 case {'mean','ma','linear','none'}, dflag = P{ix}; 0148 otherwise, warning(['Unknown option:' P{ix}]) 0149 end 0150 end 0151 Np = Np-Nk; 0152 P = {P{find(~strix)}}; % remove strings from input 0153 end 0154 0155 if Np>0 & ~isempty(P{1}), L = P{1};end 0156 if Np>1 & ~isempty(P{2}), g = P{2};end 0157 if Np>2 & ~isempty(P{3}), plotflag = P{3};end 0158 if Np>3 & ~isempty(P{4}), p = P{4};end 0159 if Np>4 & ~isempty(P{5}), chopOffHighFreq = P{5};end 0160 end 0161 0162 0163 %[L,g,plotflag,p,method,dflag,ftype,chopOffHighFreq] = d2schk(varargin); 0164 0165 0166 if (nargout == 0) & (plotflag==0) 0167 plotflag = 1; 0168 end 0169 0170 xx = xn; 0171 [n m] = size(xx); 0172 0173 if min(m,n)==1, 0174 xx = [ (1:n)' xx(:)]; 0175 n = max(m,n); 0176 m = 2; 0177 disp('Warning: The sampling frequency is undetermined and set to 1 Hz.') 0178 end 0179 dT = xx(2,1)-xx(1,1); 0180 0181 0182 if (isempty(g)), 0183 yy = xx; 0184 else 0185 yy = xx; 0186 for ix=2:m 0187 yy(:,ix)= tranproc(xx(:,ix),g); 0188 end 0189 %yy = dat2gaus(xx,g); 0190 end 0191 0192 %display('****1'); 0193 %break; 0194 switch lower(dflag) 0195 case 'mean', 0196 ma = mean(yy(:,2:m)); 0197 yy(:,2:m) = (yy(:,2:m)-ma(ones(n,1),:) ); 0198 case 'linear', 0199 yy(:,2:m) = detrend(yy(:,2:m),1); % signal toolbox detrend 0200 case 'ma', 0201 dL = ceil(1200/2/dT); % approximately 20 min. moving average 0202 yy(:,2:m) = detrendma(yy(:,2:m),dL); 0203 dflag ='mean'; 0204 end 0205 % By using a tapered data window to smooth the data at each 0206 % end of the record has the effect of sharpening 0207 % the spectral window. 0208 % NB! the resulting spectral estimate must be 0209 % normalized in order to correct for the loss of 0210 % amplitude (energy) caused by the data taper. 0211 sa = std(yy(:,2:m)); 0212 if tapery 0213 taper = bingham(n); 0214 yy(:,2:m) = taper(:,ones(1,m-1)).*yy(:,2:m); 0215 end 0216 0217 max_L = 300; % maximum lag if L is undetermined 0218 changed_L = 0; 0219 if isempty(L) 0220 L = min(n-2,ceil(4/3*max_L)); 0221 changed_L = 1; 0222 end 0223 0224 if strcmp(method,'cov') | changed_L, 0225 r=0; 0226 stdev=0; 0227 for ix=2:m 0228 R = dat2cov(yy(:,[1 ix])); 0229 r = r+R.R(:); 0230 stdev = stdev+R.stdev(:); 0231 end 0232 r = r/(m-1); 0233 R.stdev = mean(sa.^2)/r(1)*stdev/(m-1); 0234 r = r*mean(sa.^2)/r(1); 0235 R.R = r; 0236 %covplot(R,150) 0237 if changed_L, 0238 %finding where ACF is less than 2 st. deviations. 0239 L = find(abs(r(1:max_L))>2*R.stdev(1:max_L))+1; % a better L value 0240 L = min(L(end),max_L); 0241 0242 if wdef==1 % modify L so that hanning and Parzen give appr. the same result 0243 L = min(floor(4*L/3),n-2); 0244 end 0245 disp(['The default L is set to ' num2str(L) ]) 0246 end 0247 end 0248 0249 if wdef==1 % Parzen window 0250 v = floor(3.71*n/L); % degrees of freedom used in chi^2 distribution 0251 win = parzen(2*L-1); % Does not give negative estimates of the spectral density 0252 Be = 2*pi*1.33/(L*dT); % bandwidth (rad/sec) 0253 0254 elseif wdef==2 % Hanning window 0255 v = floor(2.67*n/L); % degrees of freedom used in chi^2 distribution 0256 win = hanning(2*L-1); % May give negative estimates of the spectral density 0257 Be = 2*pi/(L*dT); % bandwidth (rad/sec) 0258 end 0259 0260 nf = rate*2^nextpow2(2*L-2); % Interpolate the spectrum with rate 0261 nfft = 2*nf; 0262 0263 S = createspec('freq',ftype); 0264 S.tr = g; 0265 S.note = ['dat2spec(',inputname(1),'), Method = ' method ]; 0266 S.norm = 0; % not normalized 0267 S.L = L; 0268 0269 0270 0271 S.S = zeros(nf+1,m-1); 0272 0273 switch lower(method) 0274 case {'psd','psdo'} % from signal toolbox 0275 noverlap = 0; 0276 if method(end)=='o', noverlap = floor(L/2);end 0277 S.noverlap = noverlap; 0278 [Rper Sc f]=psd(yy(:,2),nfft,1/dT,win,noverlap,p,dflag); 0279 for ix=3:m 0280 Rper = Rper + psd(yy(:,ix),nfft,1/dT,win,noverlap,p,dflag); 0281 end 0282 Rper = Rper/(m-1); 0283 if ( ~isempty(p) ), % Confidence interval constants 0284 [maxS ind] = max(Rper); 0285 S.CI = Sc(ind,:)/Rper(ind); 0286 S.p = p; 0287 end 0288 case {'pyulear','pmem','mem'} % from signal toolbox 0289 [Rper f] = pyulear(yy(:,2),L,nfft,1/dT); 0290 for ix=3:m 0291 Rper = Rper + pyulear(yy(:,ix),L,nfft,1/dT); 0292 end 0293 Rper = Rper/(m-1); 0294 case 'pburg', % from signal toolbox 0295 [Rper f] = pburg(yy(:,2),L,nfft,1/dT); 0296 for ix=3:m 0297 Rper = Rper + pburg(yy(:,ix),L,nfft,1/dT); 0298 end 0299 Rper = Rper/(m-1); 0300 otherwise, % cov method 0301 % add a nugget effect to ensure that round off errors 0302 % do not result in negative spectral estimates 0303 r = r+nugget; 0304 rwin = r(1:L).*win(L:(2*L-1)); 0305 Rper = real(fft([rwin; zeros(nfft-(2*L-1),1); rwin(L:-1:2)],nfft)); 0306 %f = [0:(nf)]'/nf/(2*dT); 0307 if ( ~isempty(p) ), 0308 alpha = (1-p); 0309 % Confidence interval constants 0310 S.CI = [v/wchi2inv( 1-alpha/2 ,v) v/wchi2inv( alpha/2 ,v)]; 0311 S.p = p; 0312 end 0313 end 0314 0315 ind = find(Rper<0); 0316 if any(ind) 0317 Rper(ind) = 0; % set negative values to zero 0318 warning('negative spectral estimates') 0319 end 0320 0321 if strcmp(ftype,'w') 0322 S.w = [0:(nf)]'/nf*pi/dT; % (rad/s) 0323 S.S = real(Rper(1:(nf+1),1))*dT/pi; % (m^2*s/rad)one sided spectrum 0324 S.Bw = Be; 0325 else % ftype == f 0326 S.f = [0:(nf)]'/nf/2/dT; % frequency Hz if dT is in seconds 0327 S.S = 2*real(Rper(1:(nf+1),1))*dT; % (m^2*s) one sided spectrum 0328 S.Bw = Be/(2*pi); % bandwidth in Hz 0329 end 0330 0331 0332 0333 N = floor(nf/10); 0334 % cutting off high frequencies 0335 % in this way may not be a very good idea. 0336 0337 if ((N>3) & chopOffHighFreq), 0338 % The data must be Gaussian in order for this proc to be correct. 0339 [g0 test cmax irr] = dat2tr(xx,'nonlinear'); 0340 ind=(nf-4):(nf+1); 0341 0342 [Sn,m4] = wnormspec(S,0); 0343 0344 while (sqrt(m4) > irr) & (ind(1)>1) 0345 S.S(ind) = 0; 0346 [Sn,m4] = wnormspec(S,0); 0347 ind = ind-5; 0348 end 0349 fcut = S.w(min(ind(end)+5,nf)); % cut off frequency 0350 if ind(1)<1, 0351 disp('DAT2SPEC: Error in cutting off high frequencies, try other L-values') 0352 end 0353 end 0354 0355 0356 0357 %----------------------------------- 0358 % 0359 % Plotting the Spectral Density 0360 % 0361 %----------------------------------- 0362 0363 if plotflag>0 0364 wspecplot(S,plotflag) 0365 end 0366 0367 return 0368 0369 function [L,g,plotflag,p,method,dflag,ftype,chopOffHighFreq] = d2schk(P); 0370 % D2SCHK Helper function for dat2spec. 0371 % 0372 % CALL [L,g,plotflag,p,method,dflag,ftype,chopOffHighFreq]=d2schk(P) 0373 % 0374 % P = the cell array P of input arguments (between 0 and 7 elements) 0375 % xx = must be a two column vector. 0376 % 0377 0378 % Default values: 0379 %~~~~~~~~~~~~~~~~ 0380 L = []; 0381 g = []; 0382 plotflag = 0; 0383 p = 0.95; 0384 chopOffHighFreq=0; % chop off high frequencies in order to get the same 0385 % irregularity factor in the spectrum as in the data 0386 % may not be a good idea => default is 0 0387 0388 method = 'cov'; % cov. other options from signal toolbox: psd = welch's method, pyulear,pmem = 0389 % maximum entropy method 0390 dflag = 'mean'; %'ma','linear' 'mean','none' % detrending option 0391 ftype = 'w' ; %options are 'f' and 'w' 0392 0393 0394 0395 Np=length(P); 0396 strix=zeros(1,Np); 0397 for ix=1:Np, % finding symbol strings 0398 strix(ix)=ischar(P{ix}); 0399 end 0400 k = find(strix); 0401 if any(k) 0402 Nk=length(k); 0403 if Nk>3 0404 warning('More than 3 strings are not allowed in ') 0405 end 0406 for ix = k 0407 switch lower(P{ix}) 0408 case {'f','w'}, ftype = P{ix}; 0409 case {'cov','pmem','mem','psd','psdo','pburg'}, method = P{ix}; 0410 case {'mean','ma','linear','none'}, dflag = P{ix}; 0411 otherwise, warning(['Unknown option:' P{ix}]) 0412 end 0413 end 0414 Np=Np-Nk; 0415 P={P{find(~strix)}}; % remove strings from input 0416 end 0417 0418 if Np>0 & ~isempty(P{1}), L = P{1};end 0419 if Np>1 & ~isempty(P{2}), g = P{2};end 0420 if Np>2 & ~isempty(P{3}), plotflag = P{3};end 0421 if Np>3 & ~isempty(P{4}), p = P{4};end 0422 if Np>4 & ~isempty(P{5}), chopOffHighFreq = P{5};end 0423 0424 0425 return 0426
Comments or corrections to the WAFO group