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:
rf = dat2cor(x,150,2)

## CROSS-REFERENCE INFORMATION

This function calls:
 covplot Plots the auto covariance function (ACF) 1D or 2D. createcov Covariance class constructor error Display message and abort function. fft Discrete Fourier transform. ischar True for character array (string). mean Average or mean value. nextpow2 Next higher power of 2. strcmpi Compare strings ignoring case. xcorr Cross-correlation function estimates.
This function is called by:
 dat2cov Estimate auto covariance function from data.

## 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 %
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