Home > wafo > onedim > dat2spec.m

dat2spec

PURPOSE ^

Estimate one-sided spectral density from data.

SYNOPSIS ^

[S,fcut] = dat2spec(xn,varargin)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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