@@ -0,0 +1,2322 @@
+#!/usr/bin/env python%0A%0Afrom __future__ import division%0Aimport numpy as np%0Afrom scipy.linalg import hankel%0Aimport scipy.io as sio%0Aimport matplotlib.pyplot as plt%0A%0Afrom tools import *%0Afrom cum2est import *%0Afrom cum3est import *%0Afrom cum4est import *%0A%0A%0Adef cum2x(x, y, maxlag=0, nsamp=0, overlap=0, flag='biased'):%0A  %22%22%22%0A  Cross-covariance%0A  Parameters:%0A      x,y    - data vectors/matrices with identical dimensions%0A               if x,y are matrices, rather than vectors, columns are%0A               assumed to correspond to independent realizations,%0A               overlap is set to 0, and samp_seg to the row dimension.%0A      maxlag - maximum lag to be computed    %5Bdefault = 0%5D%0A    samp_seg - samples per segment  %5Bdefault = data_length%5D%0A     overlap - percentage overlap of segments %5Bdefault = 0%5D%0A               overlap is clipped to the allowed range of %5B0,99%5D.%0A       flag  - 'biased', biased estimates are computed  %5Bdefault%5D%0A              'unbiased', unbiased estimates are computed.%0A%0A  Output:%0A       y_cum - estimated cross-covariance%0A               E x%5E*(n)y(n+m),   -maxlag %3C= m %3C= maxlag%0A  %22%22%22%0A%0A  (lx, nrecs) = x.shape%0A  if (lx, nrecs) != y.shape:%0A    raise ValueError('x,y should have identical dimensions')%0A%0A  if lx == 1:%0A    lx = nrecs%0A    nrecs = 1%0A%0A  if maxlag %3C 0: raise ValueError('maxlag must be non-negative')%0A  if nrecs %3E 1: nsamp = lx%0A  if nsamp %3C= 0 or nsamp %3E lx: nsamp = lx%0A  if nrecs %3E 1: overlap = 0%0A  overlap = max(0,min(overlap,99))%0A%0A  overlap = np.fix(overlap/100 * nsamp)%0A  nadvance = nsamp - overlap%0A  if nrecs == 1:%0A    nrecs = np.fix((lx - overlap)/nadvance)%0A%0A  nlags = 2*maxlag+1%0A  zlag  = maxlag + 1%0A  y_cum = zeros(nlags,1)%0A%0A  if flag == 'biased':%0A    scale = np.ones(%5Bnlags, 1%5D)/nsamp%0A  else:%0A    scale = make_arr((range(lx-maxlag, lx-1, -1), range(lx-1, lx-maxlag, -1)), axis=1).T%0A    scale = np.ones(%5B2*maxlag+1, 1%5D) / scale%0A%0A%0A  ind = np.arange(nsamp).T%0A  for k in xrange(nrecs):%0A    xs = x%5Bind%5D%0A    xs = xs.reshape(1,-1) - np.mean(xs)%0A    ys = y%5Bind%5D%0A    ys = ys.reshape(1,-1) - np.mean(ys)%0A%0A    y_cum%5Bzlag%5D = y_cum%5Bzlag%5D + xs.T * ys%0A%0A    for m in xrange(maxlag):%0A      y_cum%5Bzlag-m%5D = y_cum%5Bzlag-m%5D + xs%5Bm:nsamp%5D.T * ys%5B0:nsamp-m%5D%0A      y_cum%5Bzlag+m%5D = y_cum%5Bzlag+m%5D + xs%5B0:nsamp-m%5D.T * ys%5Bm:nsamp%5D%0A%0A    ind = ind + int(nadvance)%0A%0A  y_cum = y_cum * scale / nrecs%0A%0A  return y_cum%0A%0A

Add Cross-covariance