Sample autocorrelation at specified lag.
    The autocorrelation is the correlation of x_i with x_{i+lag}.
336226
def autocorr(x, lag=1):
    

    if not lag:
        return 1
    if lag < 0:
        return
    # x = np.squeeze(asarray(x))
    #     mu = x.mean()
    #     v = x.var()
    #     return ((x[:-lag]-mu)*(x[lag:]-mu)).sum()/v/(len(x) - lag)
    S = autocov(x, lag)
    return S[0, 1] / sqrt(prod(diag(S)))
910968
function autocorrelation(signal) {
        
        var n=signal.length;
        
        var i;
        
        var n2=nextpow2(n);

	if (n===0){
	    n2=0;
	}
        
        var paddedSignal=[];
        
        if (n==n2){
            paddedSignal=signal;
        }
        else{
            for (i=0;i<n2;i++){
                if (i<n){
                    paddedSignal[i]=signal[i];
                }
                else{
                    paddedSignal[i]=0;
                }
            }
        }

        var acv=autocovariance(paddedSignal);
        
        var variance=acv[0];
        
        var acf=[];
        
        var i;
        
        if (variance>0){
            for (i=0; i<n; i++){
                acf[i]=acv[i]/variance;    
            }
        }
        else{
            for (i=0; i<n; i++){
                acf[i]=0;    
            }            
        }
        return acf;
    
        function autocovariance(signal) {
            
            var n=signal.length;
            
            if (n===0)
                return [];
            
            var mean=0;
            var i;
            for (i=0; i<n; i++){
                mean+=signal[i];
            }
            mean/=n;
                      
            var padded_signal=[];
            
            for (i=0; i<n; i++){
                padded_signal[i]=signal[i]-mean;
                padded_signal[n+i]=0;
            }
                        
            var ft_signal = fft(padded_signal);
                       
            var pseudo_powerSpectralDensity=[];
            
            for (i=0; i<ft_signal.length; i++){
                pseudo_powerSpectralDensity[i]=complexMultiply(ft_signal[i], complexConjugate(ft_signal[i]));
            }
           
            var pseudo_autocovariance = ifft(pseudo_powerSpectralDensity);
    
            var mask=[];
            
            for (i=0; i<n; i++){
                mask[i]=1;
                mask[i+n]=0;
            }
            
            var ft_mask = fft(mask);
    
            var mask_powerSpectralDensity=[];
            
            for (i=0; i<ft_mask.length; i++){
                mask_powerSpectralDensity[i]=complexMultiply(ft_mask[i], complexConjugate(ft_mask[i]));
            }
            
            var mask_correction_factors = ifft(mask_powerSpectralDensity);
    
            var acv=[];
            
            for (i=0; i<n; i++){
                acv[i]=complexDivideRealPart(pseudo_autocovariance[i],mask_correction_factors[i]);
            }    
            
            return acv;
            /*
             * ComplexMultiply
             * 
             * From complex.js in fft-js in https://github.com/vail-systems/node-fft
             *
             * (c) Vail Systems. Joshua Jung and Ben Bryan. 2015
             */
            function complexMultiply(a, b){
                return [(a[0] * b[0] - a[1] * b[1]), 
                        (a[0] * b[1] + a[1] * b[0])];
            }
            
            function complexConjugate(a){
                return [a[0], -a[1]];
            }
            /*function complexDivide(a,b){
                return [ (a[0]*b[0]+a[1]*b[1])/(b[0]*b[0]+b[1]*b[1]),
                         (a[1]*b[0]-a[0]*b[1])/(b[0]*b[0]+b[1]*b[1])];
            }*/
            //We only need the real part here
            function complexDivideRealPart(a,b){
                return (a[0]*b[0]+a[1]*b[1])/(b[0]*b[0]+b[1]*b[1]);
            }           
            
        }
    }
689331
def covariance(left, right, where=None, how='sample'):
    
    expr = ops.Covariance(left, right, how, where).to_expr()
    return expr
807759
def autocov(x):
    

    acorr = autocorr(x)
    varx = np.var(x, ddof=1) * (len(x) - 1) / len(x)
    acov = acorr * varx
    return acov
764701
def correlate(self, signal):
        
        s = asarray(signal)

        if s.ndim == 1:
            if size(s) != self.shape[-1]:
                raise ValueError("Length of signal '%g' does not match record length '%g'"
                                 % (size(s), self.shape[-1]))

            return self.map(lambda x: corrcoef(x, s)[0, 1], index=[1])

        elif s.ndim == 2:
            if s.shape[1] != self.shape[-1]:
                raise ValueError("Length of signal '%g' does not match record length '%g'"
                                 % (s.shape[1], self.shape[-1]))
            newindex = arange(0, s.shape[0])
            return self.map(lambda x: array([corrcoef(x, y)[0, 1] for y in s]), index=newindex)

        else:
            raise Exception('Signal to correlate with must have 1 or 2 dimensions')
493865
def at_time(self, time, nearest_sample=False):
        
        if nearest_sample:
            time += self.delta_t / 2.0
        return self[int((time-self.start_time)*self.sample_rate)]
1102124
def util_pattern_space(time_series, lag, dim):
    
    n = len(time_series)

    if lag * dim > n:
        raise Exception('Result matrix exceeded size limit, try to change lag or dim.')
    elif lag < 1:
        raise Exception('Lag should be greater or equal to 1.')

    pattern_space = np.empty((n - lag * (dim - 1), dim))
    for i in range(n - lag * (dim - 1)):
        for j in range(dim):
            pattern_space[i][j] = time_series[i + j * lag]

    return pattern_space
360917
def compamp_to_ac(compamp, window=np.hanning):  # convert single or multi-subband compamps into autocorrelation waterfall

  '''
  
  
  '''
  header = compamp.header()
 
  cdata = compamp.complex_data()

  #Apply Windowing and Padding
  cdata = np.multiply(cdata, window(cdata.shape[2]))  # window for smoothing sharp time series start/end in freq. dom.
  cdata_normal = cdata - cdata.mean(axis=2)[:, :, np.newaxis]  # zero mean, does influence a minority of lines in some plots

  cdata = np.zeros((cdata.shape[0], cdata.shape[1], 2 * cdata.shape[2]), complex)
  cdata[:, :, cdata.shape[2]/2:cdata.shape[2] + cdata.shape[2]/2] = cdata_normal  # zero-pad to 2N

  #Perform Autocorrelation
  cdata = np.fft.fftshift(np.fft.fft(cdata), 2)  # FFT all blocks separately and arrange correctly
  cdata = cdata.real**2 + cdata.imag**2  # FFT(AC(x)) = FFT(x)FFT*(x) = abs(x)^2
  cdata = np.fft.ifftshift(np.fft.ifft(cdata), 2)  # AC(x) = iFFT(abs(x)^2) and arrange correctly
  cdata = np.abs(cdata)  # magnitude of AC

  # normalize each row to sqrt of AC triangle
  cdata = np.divide(cdata, np.sqrt(np.sum(cdata, axis=2))[:, :, np.newaxis])  

  return cdata
289240
def integratedAutocorrelationTimeMultiple(A_kn, fast=False):
    

    g = statisticalInefficiencyMultiple(A_kn, fast, False)
    tau = (g - 1.0) / 2.0
    return tau
1136996
function sample_correlation(x, y) {
        var cov = sample_covariance(x, y),
            xstd = sample_standard_deviation(x),
            ystd = sample_standard_deviation(y);

        if (cov === null || xstd === null || ystd === null) {
            return null;
        }

        return cov / xstd / ystd;
    }
31657
def correlation(args):
    
    p = OptionParser(correlation.__doc__)
    opts, args, iopts = p.set_image_options(args, figsize="12x8")

    if len(args) != 1:
        sys.exit(not p.print_help())

    tsvfile, = args
    df = pd.read_csv(tsvfile, sep="\t")
    composite_correlation(df, size=(iopts.w, iopts.h))
    outfile = tsvfile.rsplit(".", 1)[0] + ".correlation.pdf"
    savefig(outfile)
61141
func AutoCorrelation(data Float64Data, lags int) (float64, error) {
	if len(data) < 1 {
		return 0, EmptyInputErr
	}

	mean, _ := Mean(data)

	var result, q float64

	for i := 0; i < lags; i++ {
		v := (data[0] - mean) * (data[0] - mean)
		for i := 1; i < len(data); i++ {
			delta0 := data[i-1] - mean
			delta1 := data[i] - mean
			q += (delta0*delta1 - q) / float64(i+1)
			v += (delta1*delta1 - v) / float64(i+1)
		}

		result = q / v
	}

	return result, nil
}
395829
async function measure(test, bench) {
  let {
    delay,
    minTime,
    minSamples,
    minWarmups,
    onCycle,
    onSample,
  } = bench.config
  let { name, run } = test

  delay *= 1e3
  minTime *= 1e3

  let samples = []
  let cycle = {
    name, // test name
    hz: null, // samples per second
    size: null, // number of samples
    time: null, // elapsed time (including delays)
    stats: null,
  }

  let n = 0
  let t = process.hrtime()
  let warmups = -1

  // synchronous test
  if (run.length == 0) {
    let start
    while (true) {
      bench.before.forEach(call)
      if (test.before) test.before()

      start = process.hrtime()
      run()
      let sample = clock(start)

      if (test.after) test.after()
      bench.after.forEach(call)

      if (warmups == -1) {
        warmups = Math.max(minWarmups, sample > 100 ? 1 : sample > 10 ? 5 : 50)
      }
      if (warmups > 0) {
        warmups--
      } else {
        samples[n++] = sample
        onSample(samples[n - 1], cycle)

        if (minTime <= clock(t) && minSamples <= n) {
          break // all done!
        }
      }

      // wait then repeat
      await wait(delay)
    }
  }

  // asynchronous test
  else {
    let start
    const next = function() {
      bench.before.forEach(call)
      if (test.before) test.before()

      start = process.hrtime()
      run(done)
    }

    // called by the test function for every sample
    const done = function() {
      let sample = clock(start)

      if (test.after) test.after()
      bench.after.forEach(call)

      if (warmups == -1) {
        warmups = Math.max(minWarmups, sample > 100 ? 1 : sample > 10 ? 5 : 50)
      }
      if (warmups > 0) {
        warmups--
      } else {
        samples[n++] = sample
        onSample(samples[n - 1], cycle)

        if (minTime <= clock(t) && minSamples <= n) {
          return cycle.done() // all done!
        }
      }

      // wait then repeat
      wait(delay).then(next)
    }

    defer(cycle) // wrap cycle with a promise
    next() // get the first sample

    // wait for samples
    await cycle.promise
  }

  const time = clock(t)
  const stats = analyze(samples)
  cycle.hz = 1000 / stats.mean
  cycle.size = n
  cycle.time = time
  cycle.stats = stats

  onCycle(cycle, bench.config)
}
1041228
func (f Float64Data) AutoCorrelation(lags int) (float64, error) {
	return AutoCorrelation(f, lags)
}
121980
public static double computeCorrelation(final List<Double> list1, final List<Double> list2)
    {
        final double[] doubleArray1 = new double[list1.size()];
        final double[] doubleArray2 = new double[list2.size()];

        int off1 = 0;
        for (final double item : list1) {
            doubleArray1[off1] = item;
            off1++;
        }

        int off2 = 0;
        for (final double item : list2) {
            doubleArray2[off2] = item;
            off2++;
        }

        final double correlation = computeCorrelation(doubleArray1, doubleArray2);

        return correlation;
    }
19690
def _find_max_lag(x, rho_limit=0.05, maxmaxlag=20000, verbose=0):
    """"""

    # Fetch autocovariance matrix
    acv = autocov(x)
    # Calculate rho
    rho = acv[0, 1] / acv[0, 0]

    lam = -1. / np.log(abs(rho))

    # Initial guess at 1.5 times lambda (i.e. 3 times mean life)
    maxlag = int(np.floor(3. * lam)) + 1

    # Jump forward 1% of lambda to look for rholimit threshold
    jump = int(np.ceil(0.01 * lam)) + 1

    T = len(x)

    while ((abs(rho) > rho_limit) & (maxlag < min(T / 2, maxmaxlag))):

        acv = autocov(x, maxlag)
        rho = acv[0, 1] / acv[0, 0]
        maxlag += jump

    # Add 30% for good measure
    maxlag = int(np.floor(1.3 * maxlag))

    if maxlag >= min(T / 2, maxmaxlag):
        maxlag = min(min(T / 2, maxlag), maxmaxlag)
        "maxlag fixed to %d" % maxlag
        return maxlag

    if maxlag <= 1:
        print_("maxlag = %d, fixing value to 10" % maxlag)
        return 10

    if verbose:
        print_("maxlag = %d" % maxlag)
    return maxlag
1261454
def crosscorr(self, signal, lag=0):
        
        from scipy.linalg import norm

        s = asarray(signal)
        s = s - mean(s)
        s = s / norm(s)

        if size(s) != size(self.index):
            raise Exception('Size of signal to cross correlate with, %g, '
                            'does not match size of series' % size(s))

        # created a matrix with lagged signals
        if lag is not 0:
            shifts = range(-lag, lag+1)
            d = len(s)
            m = len(shifts)
            sshifted = zeros((m, d))
            for i in range(0, len(shifts)):
                tmp = roll(s, shifts[i])
                if shifts[i] < 0:
                    tmp[(d+shifts[i]):] = 0
                if shifts[i] > 0:
                    tmp[:shifts[i]] = 0
                sshifted[i, :] = tmp
            s = sshifted
        else:
            shifts = [0]

        def get(y, s):
            y = y - mean(y)
            n = norm(y)
            if n == 0:
                b = zeros((s.shape[0],))
            else:
                y /= n
                b = dot(s, y)
            return b

        return self.map(lambda x: get(x, s), index=shifts)