Home > wafo > onedim > dat2cor.m

dat2cor

PURPOSE ^

Estimate auto correlation function from data.

SYNOPSIS ^

[corr,c0] = dat2cor(xn,varargin)

DESCRIPTION ^

 DAT2COR Estimate auto correlation function from data.
 
  CALL:  r = dat2cor(x,L,plotflag,dT,flag)
 
          r = Covariance 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(corr(k))=1/N*(1+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 1/N for k>0
                  h = water depth (default inf)
                 tr = [], transformation 
               type = 'none'
               norm = 1 indicating that R is normalized
 
         x  = a vector or a two column data matrix with 
              sampled times and values. (length n)
         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 correlation by 1/n. (default)
              'unbiased': scales the raw correlation by 1/(n-abs(k)), 
                          where k is the lag.
 
   Note: - flag may be given anywhere after x.
         - x may contain NaN's (i.e. missing values).
 
  Example: 
    x = load('sea.dat');
    rf = dat2cor(x,150,2)
 
  See also  dat2cov, covplot

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function [corr,c0] = dat2cor(xn,varargin)
002 %DAT2COR Estimate auto correlation function from data.
003 %
004 % CALL:  r = dat2cor(x,L,plotflag,dT,flag)
005 %
006 %         r = Covariance 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(corr(k))=1/N*(1+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 1/N for k>0
016 %                 h = water depth (default inf)
017 %                tr = [], transformation 
018 %              type = 'none'
019 %              norm = 1 indicating that R is normalized
020 %
021 %        x  = a vector or a two column data matrix with 
022 %             sampled times and values. (length n)
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 correlation by 1/n. (default)
030 %             'unbiased': scales the raw correlation by 1/(n-abs(k)), 
031 %                         where k is the lag.
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 = dat2cor(x,150,2)
039 %
040 % See also  dat2cov, covplot
041 
042        
043 %(If you do have the signal toolbox you may modify this file)
044 
045 % Tested on: Matlab 6.0, 5.3, 5.2, 5.1
046 % History:
047 % revised jr 02.04.2001
048 %  - added example, updating help 
049 % revised pab 09.10.2000
050 %  - added secret output c0 = variance
051 %  - made computations faster for the case Ncens < n
052 %  - added flag option, dat2corchk
053 %  -fixed a bug: forgot to subtract the mean from x
054 % revised pab 12.08.99
055 % - changed code to handle covariance structure/class
056 % modified by Per A. Brodtkorb 28.10.98 
057 % enabled to handle missing data 
058 
059 
060 [x,L,plotflag,dT,flag,n] = dat2corchk(varargin,xn,nargout);
061 
062 corr       = createcov(0,'t');
063 corr.R     = zeros(L+1,1);
064 corr.t     = zeros(L+1,1);
065 corr.stdev = zeros(L+1,1);
066 corr.h     = inf;
067 corr.norm  = 1; % normalized
068 
069 %----------------------------------------------
070 %  
071 %  Signal toolbox users may modify below
072 %
073 %----------------------------------------------
074 
075 indnan = isnan(x);
076 if any(indnan)
077    x     = x-mean(x(~indnan)); % remove the mean pab 09.10.2000
078    indnan = find(indnan);
079    Ncens = n-length(indnan);
080 else
081    indnan = [];
082    Ncens  = n;
083    x      = x-mean(x);
084 end
085 
086 if 0 %Ncens<n, %Old call for missing data (slow)   
087    r=zeros(L+1,1);
088    % unbiased result, i.e. divide by n-abs(lag)
089    for tau=0:L,
090       Ntau=n-tau;
091       ix=1:Ntau; % constructing a vector if indices
092       %vector multiplication is faster than a for loop in MATLAB
093       tmp=x(ix).*x(ix+tau);
094       if strcmpi(flag,'unbiased'),
095          r(tau+1)=mean(tmp(~isnan(tmp))); % unbiased estimate
096       else
097          r(tau+1)=sum(tmp(~isnan(tmp)))/Ncens; % biased estimate
098       end
099    end
100 else 
101    if any(indnan)
102       x(indnan)=0; % pab 09.10.2000 much faster for censored samples
103    end
104    if 1, % If you haven't got the signal toolbox use this
105      nfft=2^nextpow2(n) ;
106      Rper=abs(fft(x,nfft)).^2/Ncens; % Raw periodogram
107      %plot(Rper)%,size(Rper),n,nfft
108           
109      r=real(fft(Rper))/nfft ; %ifft=fft/nfft since Rper is real!
110       
111      if strcmpi(flag,'unbiased'),
112            % unbiased result, i.e. divide by n-abs(lag)
113          r=r(1:L+1)*Ncens./ (Ncens:-1:Ncens-L)';
114      %else  % biased result, i.e. divide by n
115      %   r=r(1:L+1)/Ncens;
116      end
117 
118   else % If you have got the signal toolbox use this 
119     r=xcorr(x,L,flag);
120     r=(r((L+1):end))'; 
121   end  
122 end
123 
124 c0     = r(1);
125 corr.R = r(1:(L+1))/c0; %normalizing
126 corr.t = (0:L)'*dT;
127 
128 corr.stdev=sqrt(1/Ncens*[ 0; 1 ;(1+2*cumsum(corr.R(2:L).^2))]);
129 
130 if ((plotflag>0) ),
131   covplot(corr,L,plotflag)
132 end
133 
134 return
135 
136 function [x,L,plotflag,dT,flag,n] = dat2corchk(P,xn,Na)
137 % DAT2CORCHK Helper function for dat2cor.
138 %
139 % CALL  [x2,L,plotflag,dT,flag,n]=dat2corchk(P,x1) 
140 %
141 %   P = the cell array P of input arguments (between 0 and 7 elements)
142 %   x = a vector.
143 
144 
145 x=xn;
146 [n m]= size(x);
147 
148 if n<m
149  b=m;m=n;n=b; 
150  x=x';
151 end
152 
153 if n<2, 
154   error('The vector must have more than 2 elements!')
155 end
156 
157 % Default values
158 %~~~~~~~~~~~~~~~~  
159 flag = 'biased';
160 plotflag = 0;
161 dT=[]; % sampling frequency unknown
162 L=n-1;
163    % if n<400
164 %    L=floor(n/2);
165 %  else
166 %    L=floor(12*sqrt(n));
167 %  end
168 
169 
170 
171 
172 Np=length(P);
173 strix=zeros(1,Np);
174 for ix=1:Np, % finding symbol strings 
175  strix(ix)=ischar(P{ix});
176 end
177 
178 
179 k=find(strix);
180 if any(k)
181   flag = P{k(1)};
182   Np   = Np-1;
183   P={P{find(~strix)}};
184 end
185 
186 
187 if Np>0 & ~isempty(P{1}), L=P{1}; end
188 if Np>1 & ~isempty(P{2}), plotflag=P{2}; end
189 if Np>2 & ~isempty(P{3}), dT=P{3}; end
190 
191 if (Na == 0)&plotflag==0
192   plotflag=3;
193 end
194 if plotflag>3, 
195   error('Invalid option. Plotflag can only be 0,1,2 or 3')
196 end
197  switch m
198  case 1, x=x(:); if isempty(dT),dT=1;end
199  case 2,  dT=x(2,1)-x(1,1);x=x(:,2);
200  otherwise, error('Wrong dimension of input! dim must be 2xN, 1xN, Nx2 or Nx1 ') 
201               
202 end
203 return
204

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