% calculate cross correlations % function [lag, R_xy, P_xy, n,Sig_xy] = cross_corr(x, y, nlags, bad_flag) % % input: % x,y the data sets % nlags are number of lags to process % bad_flag is data value that indicates missing or bad data % % output: % lag is lag from -nlag to +nlag % R_xy are covariances. % P_xy are correlations % n is array containing the number of data points used to % calculate each R % put in column format x=x(:); y=y(:); % this many data points N=length(x); % initialize output R_xy = zeros(2*nlags+1,1); P_xy = zeros(2*nlags+1,1); Sig_xy = zeros(2*nlags+1,1); n = zeros(2*nlags+1,1); % check that y is same length if (length(y)~=N) fprintf('x and y different lengths\n') return; end % find means if(isnan(bad_flag)) bad_flag=1e35; id = find(isnan(x)); x(id) = bad_flag+0*id; id = find(isnan(y)); y(id) = bad_flag+0*id; end good_id = find(x~=bad_flag); if(length(good_id)>0) mean_x = mean(x(good_id)); else fprintf('no data found\n') return; end good_id = find(y~=bad_flag); if(length(good_id)>0) mean_y = mean(y(good_id)); else fprintf('no data found\n') return; end % do the lags cnt = 0; for l=-nlags:1:nlags, cnt = cnt + 1; % check for neg./pos lag if (l<0) k=(-1)*l; lag2_id = [1:1:(N-k)]'; lag1_id = lag2_id+k; else k=l; lag1_id = [1:1:(N-k)]'; lag2_id = lag1_id+k; end % find good data in x series good_id = find( (x(lag1_id)~=bad_flag) ); Ngoodx = length(good_id); % continue with this lag if there are data if (Ngoodx>0) lag1_id = lag1_id(good_id); lag2_id = lag2_id(good_id); % find good data in y-series where x series was good good_id = find( (y(lag2_id)~=bad_flag) ); Ngood = length(good_id); % continue only if there are data if (Ngood>0) n(cnt) = Ngood; lag1_id = lag1_id(good_id); lag2_id = lag2_id(good_id); % calculate statistics % dudley's method if (1) mean_1 = mean(x(lag1_id)); mean_2 = mean(y(lag2_id)); z1=x(lag1_id)-mean_1; z2=y(lag2_id)-mean_2; end % nathaniel's method if (0) z1=x(lag1_id)-mean_x; z2=y(lag2_id)-mean_y; end % get the normalizing variances std_1 = sqrt(z1'*z1/Ngood ); std_2 = sqrt(z2'*z2/Ngood); % estimate cov. and corr. R_xy(cnt) = z1'*z2/Ngood; % Equation 6.2.78 P_xy(cnt) = R_xy(cnt)/(std_1*std_2);% Equation 6.2.82 % use matlab function to get P values [R,P]= corrcoef(z1,z2); Sig_xy(cnt)=P(1,2); end % check for good x -data end % check for good y -data end % all lags % reproduce lags lag = [-nlags:1:nlags]';