DAT2SPEC2 Estimate one-sided spectral density, version 2. using a Parzen window function on the estimated autocovariance function. CALL: S = dat2spec2(x,L,g,plotflag,wdef,dT) 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. 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 wdef = 'parzen' then a Parzen window is used (default). 'hanning' then a Hanning window is used. dT = sampling interval. Default is dT=x(2,1)-x(1,1) or 1Hz. 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. See also dat2tr, dat2cov, dat2spec
Spectrum structure constructor | |
Estimate auto covariance function from data. | |
Transforms x using the transformation g. | |
Estimate transformation, g, from data. | |
returns the N-point Hanning window in a column vector. | |
returns the N-point Parzen window in a column vector. | |
Normalize a spectral density such that m0=m2=1 | |
Plot a spectral density | |
Display message and abort function. | |
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) | |
Compare strings. |
% CHAPTER2 Modelling random loads and stochastic waves | |
% CHAPTER3 Demonstrates distributions of wave characteristics |
001 function S = dat2spec2(xn,L,g,plotflag,wdef,dT,chopOffHighFreq) 002 %DAT2SPEC2 Estimate one-sided spectral density, version 2. 003 % using a Parzen window function on 004 % the estimated autocovariance function. 005 % 006 % CALL: S = dat2spec2(x,L,g,plotflag,wdef,dT) 007 % 008 % S = A structure containing: 009 % S = spectral density 010 % w = angular frequency 011 % tr = transformation g 012 % h = water depth (default inf) 013 % type = 'freq' 014 % note = Memorandum string 015 % date = Date and time of creation 016 % L = maximum lag size of the window function. 017 % Bw = Bandwidth of the smoothing window which is used 018 % in the estimated spectrum. (rad/sec or Hz) 019 % 020 % 021 % x = m column data matrix with sampled times in the first column 022 % and values the next columns. 023 % 024 % L = maximum lag size of the window function. 025 % If no value is given, the lag size is set to 026 % be the lag where the auto correlation is less than 027 % 2 standard deviations. (maximum 300) 028 % 029 % g = the transformation assuming that x is a sample of a 030 % transformed Gaussian process. If g is empty then 031 % x is a sample of a Gaussian process (Default) 032 % 033 % plotflag = 1 plots the spectrum, S, 034 % 2 plot 10log10(S) and 035 % 3 plots both the above plots 036 % 037 % wdef = 'parzen' then a Parzen window is used (default). 038 % 'hanning' then a Hanning window is used. 039 % 040 % dT = sampling interval. Default is dT=x(2,1)-x(1,1) or 1Hz. 041 % 042 % As L decreases the estimate becomes smoother and Bw increases. 043 % If we want to resolve peaks in S which is Bf (Hz or rad/sec) apart 044 % then Bw < Bf. 045 % 046 % See also dat2tr, dat2cov, dat2spec 047 048 % May chop off high frequencies in order to 049 % get the same irregularity factor in the spectrum as 050 % in x. 051 052 053 % References: 054 % Georg Lindgren and Holger Rootzen (1986) 055 % "Stationära stokastiska processer", pp 173--176. 056 % 057 % Gareth Janacek and Louise Swift (1993) 058 % "TIME SERIES forecasting, simulation, applications", 059 % pp 75--76 and 261--268 060 % 061 % Emanuel Parzen (1962), 062 % "Stochastic Processes", HOLDEN-DAY, 063 % pp 66--103 064 065 066 % Tested on: Matlab 5.3 067 % history: 068 % Modified by Per A. Brodtkorb 14.08.98,25.05.98 069 % - add a nugget effect to ensure that round off errors 070 % do not result in negative spectral estimates 071 % Modified by svi 29.09.99 072 % - The program is not estimating transformation g any more. 073 % modified pab 22.10.1999 074 % - fixed so that x in fact can be a m column matrix 075 % - updated info put into the spectral structure. 076 % - updated help header 077 % modified pab 03.11.1999 078 % - fixed a bug: line 152 wrong array dim. when m>2 079 % modified jr 29.01.2000 080 % - modified from dat2spec, new name: dat2spec2 081 % - lines related to confidence bands removed 082 % - call to psd removed 083 % Revised pab Dec 2003 084 % -fixed a bug for wdef and updated helpheader accordingly. 085 % 086 087 nugget=10^-12; 088 089 freqtype='w'; %options are 'f' and 'w' 090 091 xx=xn; 092 istime=1; 093 [n m]= size(xx); 094 095 if min(m,n)==1, 096 xx=[ (1:n)' xx(:)];istime=0; 097 n=max(m,n); 098 m=2; 099 end 100 101 if (nargin < 7) | isempty(chopOffHighFreq), 102 chopOffHighFreq=0; 103 % chop off high frequencies in order to get the same 104 % irregularity factor in the spectrum as in the data 105 % not a good idea 106 end 107 108 if (nargin < 6) | isempty(dT), 109 dT=xx(2,1)-xx(1,1); 110 if ~istime 111 disp(' Warning: The sampling frequency is undetermined and set to 1 Hz.') 112 end 113 end 114 115 if (nargin < 5) | isempty(wdef) 116 wdef=1; % 1=parzen window 2=hanning window 117 elseif ischar(wdef) 118 switch lower(wdef(1)) 119 case 'p', 120 wdef = 1; % parzen 121 case 'h', 122 wdef = 2; % hanning 123 otherwise 124 error('unknown option') 125 end 126 end 127 128 if (nargin < 4) | isempty(plotflag) 129 plotflag = 0; 130 end 131 if (nargout == 0) & (plotflag==0) 132 plotflag = 1; 133 end 134 135 if (nargin<3) 136 g=[]; 137 end 138 139 if (isempty(g)), 140 yy=xx; 141 else 142 yy = dat2gaus(xx,g); 143 end 144 145 % By using a tapered data window to smooth the data at each 146 % end of the record has the effect of sharpening 147 % the spectral window. 148 % NB! the resulting spectral estimate must be 149 % normalized in order to correct for the loss of 150 % amplitude (energy) caused by the data taper. 151 %taper=bingham(n); 152 %yy(:,2)=bingham(n).*yy(:,2); 153 %display('****1'); 154 %break; 155 yy(:,2:m)=(yy(:,2:m)-ones(n,1)*mean(yy(:,2:m)) );%/std(yy(:,2)); 156 157 max_L=300; 158 159 changed_L =0; 160 if (nargin < 2 | isempty(L)) 161 L=min(n-2,ceil(4/3*max_L)); 162 %else 163 changed_L =1; 164 end 165 166 %if ~strcmp(method,'psd') | changed_L, 167 r=0; 168 stdev=0; 169 for ix=2:m 170 R=dat2cov(yy(:,[1 ix])); 171 r=r+R.R(:); 172 stdev=stdev+R.stdev(:); 173 end 174 R.stdev=stdev/(m-1); 175 r=r/(m-1); 176 R.R=r; 177 if changed_L, 178 % finding where ACF is less than 2 st. deviations. 179 L=find(abs(r(1:max_L))>2*R.stdev(1:max_L))+1; % a better L value 180 L=min(L(end),max_L); 181 end 182 %end 183 184 if wdef==1 % Parzen window 185 if changed_L, % modify L so that hanning and Parzen give appr. the same 186 % result 187 L=min(floor(4*L/3),n-2);disp(['The default L is set to ' num2str(L) ]) 188 end 189 v=floor(3.71*n/L); % degrees of freedom used in chi^2 distribution 190 win=parzen(2*L-1); % Does not give negative estimates of the spectral density 191 Be=2*pi*1.33/(L*dT); % bandwidth (rad/sec) 192 193 elseif wdef==2 % Hanning window 194 if changed_L, 195 disp([' The default L is set to ' num2str(L) ]) 196 end 197 v=floor(2.67*n/L); % degrees of freedom used in chi^2 distribution 198 win=hanning(2*L-1); % May give negative estimates of the spectral density 199 Be =2*pi/(L*dT); % bandwidth (rad/sec) 200 end 201 202 nf=2^nextpow2(2*L-2); % Interpolate the spectrum with rate = 2 203 nfft=2*nf; 204 205 S=createspec('freq',freqtype); 206 S.tr=g; 207 S.note=['dat2spec(',inputname(1),')']; 208 S.L=L; 209 S.S=zeros(nf+1,m-1); 210 211 % add a nugget effect to ensure that round off errors 212 % do not result in negative spectral estimates 213 214 r=r+nugget; 215 rwin=r(1:L).*win(L:(2*L-1)); 216 Rper=real(fft([rwin; zeros(nfft-(2*L-1),1); rwin(L:-1:2)],nfft)); 217 218 rpmin=min(Rper); 219 if rpmin<0 220 disp(' Warning: negative spectral estimates') 221 end 222 223 if strcmp(freqtype,'w') 224 S.w=[0:(nf)]'/nf*pi/dT; % (rad/s) 225 S.S=abs(Rper(1:(nf+1),1))*dT/pi; % (m^2*s/rad) one-sided spectrum 226 S.Bw=Be; 227 else % freqtype == f 228 S.f=[0:(nf)]'/nf/2/dT; % frequency Hz if dT is in seconds 229 S.S=2*abs(Rper(1:(nf+1),1))*dT; % (m^2*s) one-sided spectrum 230 S.Bw=Be/(2*pi); % bandwidth in Hz 231 end 232 233 234 235 236 N=floor(nf/10); 237 % cutting off high frequencies 238 % in this way may not be a very good idea. 239 240 if (N>3)&chopOffHighFreq, 241 242 % The data must be Gaussian in order for this proc to be correct. 243 [g0 test cmax irr] = dat2tr(xx,'nonlinear'); 244 ind=(nf-4):(nf+1); 245 [S,m4,m2]=wnormspec(S,0); 246 247 while (sqrt(m4/m2^2) > irr) & (ind(1)>1) 248 S.S(ind)=0; 249 [S,m4,m2]=wnormspec(S,0); 250 ind=ind-5; 251 end 252 fcut=S.w(min(ind(end)+5,nf)); % cut off frequency 253 if ind(1)<1, 254 disp('DAT2SPEC: Error in cutting off high frequencies, try other L-values') 255 end 256 end 257 258 259 260 %----------------------------------- 261 % 262 % Plotting the Spectral Density 263 % 264 %----------------------------------- 265 266 if plotflag>0 267 wspecplot(S,plotflag) 268 end 269 270
Comments or corrections to the WAFO group