# -*- coding: utf-8 -*-
# 2009-12-07 Moritz

import scipy.signal as sg
import matplotlib.pyplot as plt
import numpy as np
plt.close('all')

def xcorr(tr1, tr2):
    """
    Cross Correlation, normalized.

    Do the convolution for more than 5000 points in the frequency domain,
    else in the time domain.
    """
    denom = np.sqrt( np.dot(tr1,tr1)*np.dot(tr2,tr2) )
    if len(tr1)+len(tr2) < 5000:
        return np.convolve(tr1[::-1], tr2)/denom
    else:
        return sg.fftconvolve(tr1[::-1], tr2)/denom


# Winlen
N = 1000


# Boxcar
data = np.zeros(N)
data[300:400] = 1.0

# Uncomment for using random data
#data = np.random.rand(N) - .5

# Do the auto correlation
autocorr = xcorr(data, data)


#
# The plotting part
#
plt.figure(1)
plt.subplot(211)
plt.plot(boxcar, 'k', label='Original Data')
plt.ylim([plt.ylim()[0]-0.1, plt.ylim()[1]+0.1])
plt.legend()
plt.subplot(212)
plt.plot(autocorr, 'b', label='Max Xcorr Value')
plt.legend()
plt.xlabel("Time [s]")

plt.show()