statistics.py 950 Bytes
Newer Older
Dion Haefner's avatar
Dion Haefner committed
1
import numpy as np
Dion Haefner's avatar
Dion Haefner committed
2
from scipy import integrate, signal
Dion Haefner's avatar
Dion Haefner committed
3

Dion Haefner's avatar
Dion Haefner committed
4 5 6 7 8 9 10 11 12
def norm_L2(a,x):
	a_squared = np.squeeze(a**2)
	n = len(np.shape(a_squared))
	for i in range(n):
		a_squared = integrate.simps(a_squared,x=x[i],axis=0)
	return np.sqrt(a_squared)


def autocorrelation(H, pad=None):
Dion Haefner's avatar
Dion Haefner committed
13 14 15 16 17 18
	"""
	Calculates the autocorrelation of an array via ffts.

	.. math::

		\\text{autocovariance}(H,x) = \\text{fft}^{-1}(|\\text{fft}(H - \\overline{H},k)|^2,x)
Dion Haefner's avatar
Dion Haefner committed
19

Dion Haefner's avatar
Dion Haefner committed
20 21 22
		\\text{autocorrelation}(H,x) = \\text{autocovariance}(x) / \\text{autocovariance}(0)

	:param H: input array
Dion Haefner's avatar
Dion Haefner committed
23
	:param pad: Pads each axis of the input with pad[i] zeros before and after each FFT.
Dion Haefner's avatar
Dion Haefner committed
24 25

	"""
Dion Haefner's avatar
Dion Haefner committed
26 27 28 29
	#pad = np.array(H.shape) if pad is None else pad
	#H = H - np.mean(H)
	#ft_signal = np.fft.fftn(H,pad)
	#B = np.fft.ifftn(np.abs(ft_signal)**2,pad)
30
	reverse = [slice(n,0,-1) for n in H.shape]
Dion Haefner's avatar
Dion Haefner committed
31 32
	autocovariance = signal.fftconvolve(H,H[reverse],mode="full")[reverse]
	return autocovariance / autocovariance.flat[0]