Home > wafo > onedim > dat2cov.m

dat2cov

PURPOSE ^

Estimate auto covariance function from data.

SYNOPSIS ^

R = dat2cov(xn,varargin)

DESCRIPTION ^

 DAT2COV Estimate auto covariance function from data.
 
  CALL:  R = dat2cov(x,L,plotflag,dT,flag)
 
          R = a structure containing:
                  R = ACF vector length L+1  
                  t = time lags  length L+1  
              stdev = estimated large lag standard deviation of the estimate
                      assuming x is a Gaussian process:
                      if R(k)=0 for all lags k>q then an approximation 
                      of the variance for large samples due to Bartlett
                      var(R(k))=1/N*(R(0)^2+2*R(1)^2+2*R(2)^2+ ..+2*R(q)^2)
                      for  k>q and where  N=length(x). Special case is
                      white noise where it equals R(0)^2/N for k>0
                  h = water depth (default inf)
                 tr = [], transformation 
               type = 'none'
               norm = 0 indicating that R is not normalized
 
         x  = a column data vector or
              two column data matrix with sampled times and values.
         L  = the maximum time-lag for which the ACF is estimated.
                (Default L=n-1)        
  plotflag = 1 then the ACF is plotted vs lag
             2 then the ACF is plotted vs lag in seconds
             3 then the ACF is plotted vs lag and vs lag (sec
        dT = time-step between data points (default xn(2,1)-xn(1,1) or 1 Hz).
      flag = 'biased'  : scales the raw cross-correlation by 1/n. (default)
             'unbiased': scales the raw correlation by 1/(n-abs(k)), 
                         where k is the index into the result.
 
  Note:  - flag may be given anywhere after x.
         - x may contain NaN's  (i.e. missing values).
 
  Example: 
    x = load('sea.dat');
    rf = dat2cov(x,150,2)
 
  See also  dat2cor, covplot

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function R = dat2cov(xn,varargin)
002 %DAT2COV Estimate auto covariance function from data.
003 %
004 % CALL:  R = dat2cov(x,L,plotflag,dT,flag)
005 %
006 %         R = a structure containing:
007 %                 R = ACF vector length L+1  
008 %                 t = time lags  length L+1  
009 %             stdev = estimated large lag standard deviation of the estimate
010 %                     assuming x is a Gaussian process:
011 %                     if R(k)=0 for all lags k>q then an approximation 
012 %                     of the variance for large samples due to Bartlett
013 %                     var(R(k))=1/N*(R(0)^2+2*R(1)^2+2*R(2)^2+ ..+2*R(q)^2)
014 %                     for  k>q and where  N=length(x). Special case is
015 %                     white noise where it equals R(0)^2/N for k>0
016 %                 h = water depth (default inf)
017 %                tr = [], transformation 
018 %              type = 'none'
019 %              norm = 0 indicating that R is not normalized
020 %
021 %        x  = a column data vector or
022 %             two column data matrix with sampled times and values.
023 %        L  = the maximum time-lag for which the ACF is estimated.
024 %               (Default L=n-1)        
025 % plotflag = 1 then the ACF is plotted vs lag
026 %            2 then the ACF is plotted vs lag in seconds
027 %            3 then the ACF is plotted vs lag and vs lag (sec
028 %       dT = time-step between data points (default xn(2,1)-xn(1,1) or 1 Hz).
029 %     flag = 'biased'  : scales the raw cross-correlation by 1/n. (default)
030 %            'unbiased': scales the raw correlation by 1/(n-abs(k)), 
031 %                        where k is the index into the result.
032 %
033 % Note:  - flag may be given anywhere after x.
034 %        - x may contain NaN's  (i.e. missing values).
035 %
036 % Example: 
037 %   x = load('sea.dat');
038 %   rf = dat2cov(x,150,2)
039 %
040 % See also  dat2cor, covplot
041 
042 %(If you have got the signal toolbox you may modify the dat2cor file)%--------
043 
044 % Tested on: Matlab 6.0, 5.3, 5.2, 5.1
045 % History:
046 % revised jr 02.04.2001
047 %  - added example, updating help 
048 % revised pab 09.10.2000
049 %  - added dat2corchk
050 %  - updated documentation
051 %  - added flag option
052 % revised pab 12.08.99
053 %  - updated to handle covariance class/structure
054 % modified by Per A. Brodtkorb 28.10.98 
055 %  removed some code and replaced it wih a call to dat2cor
056 
057 
058 [x,L,plotflag,dT,flag,n] = dat2corchk(varargin,xn,nargout);
059 
060  
061 % ind     = find(~isnan(x));
062 %n0      = length(ind);
063 %c0      = (n0-1)/n0*(std(x(ind)).^2);% allowing missing data
064 [R, c0] = dat2cor(x,L,0,dT,flag);
065 R.R     = R.R*c0;
066 R.stdev = R.stdev*c0;
067 R.norm  = 0;
068 
069 if ((plotflag>0) ),
070   covplot(R,L,plotflag)
071 end
072 
073 
074 function [x,L,plotflag,dT,flag,n] = dat2corchk(P,xn,Na)
075 % DAT2CORCHK Helper function for dat2cor.
076 %
077 % CALL  [x2,L,plotflag,dT,flag,n]=dat2corchk(P,x1) 
078 %
079 %   P = the cell array P of input arguments (between 0 and 7 elements)
080 %   x = a vector.
081 
082 
083 x=xn;
084 [n m]= size(x);
085 
086 if n<m
087  b=m;m=n;n=b; 
088  x=x';
089 end
090 
091 if n<2, 
092   error('The vector must have more than 2 elements!')
093 end
094 
095 % Default values
096 %~~~~~~~~~~~~~~~~  
097 flag = 'biased';
098 plotflag = 0;
099 dT=[]; % sampling frequency unknown
100 L=n-1;
101    % if n<400
102 %    L=floor(n/2);
103 %  else
104 %    L=floor(12*sqrt(n));
105 %  end
106 
107 
108 
109 
110 Np=length(P);
111 strix=zeros(1,Np);
112 for ix=1:Np, % finding symbol strings 
113  strix(ix)=ischar(P{ix});
114 end
115 
116 
117 k=find(strix);
118 if any(k)
119   flag = P{k(1)};
120   Np   = Np-1;
121   P={P{find(~strix)}};
122 end
123 
124 if Np>0 & ~isempty(P{1}), L=P{1}; end
125 if Np>1 & ~isempty(P{2}), plotflag=P{2}; end
126 if Np>2 & ~isempty(P{3}), dT=P{3}; end
127 
128 if (Na == 0)&plotflag==0
129   plotflag=3;
130 end
131 if plotflag>3, 
132   error('Invalid option. Plotflag can only be 0,1,2 or 3')
133 end
134  switch m
135  case 1, x=x(:); if isempty(dT),dT=1;end
136  case 2,  dT=x(2,1)-x(1,1);x=x(:,2);
137  otherwise, error('Wrong dimension of input! dim must be 2xN, 1xN, Nx2 or Nx1 ') 
138               
139 end
140 return
141

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