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
Plots the auto covariance function (ACF) 1D or 2D. | |
Covariance class constructor | |
Display message and abort function. | |
Discrete Fourier transform. | |
True for character array (string). | |
Average or mean value. | |
Next higher power of 2. | |
Compare strings ignoring case. | |
Cross-correlation function estimates. |
Estimate auto covariance function from data. |
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
Comments or corrections to the WAFO group