signal iir help for fir filtros filtro and python numpy signal-processing

python - iir - scipy help



Estimar la autocorrelación usando Python (5)

Me gustaría realizar la autocorrelación en la señal que se muestra a continuación. El tiempo entre dos puntos consecutivos es 2.5ms (o una tasa de repetición de 400Hz).

Esta es la ecuación para estimar la autoacrrelación que me gustaría usar (tomada de http://en.wikipedia.org/wiki/Autocorrelation , sección Estimación)

¿Cuál es el método más simple para encontrar la autocorrelación estimada de mis datos en python? ¿Hay algo similar a numpy.correlate que pueda usar?

¿O debería simplemente calcular la media y la varianza?

Editar:

Con la ayuda de unutbu , he escrito:

from numpy import * import numpy as N import pylab as P fn = ''data.txt'' x = loadtxt(fn,unpack=True,usecols=[1]) time = loadtxt(fn,unpack=True,usecols=[0]) def estimated_autocorrelation(x): n = len(x) variance = x.var() x = x-x.mean() r = N.correlate(x, x, mode = ''full'')[-n:] #assert N.allclose(r, N.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)])) result = r/(variance*(N.arange(n, 0, -1))) return result P.plot(time,estimated_autocorrelation(x)) P.xlabel(''time (s)'') P.ylabel(''autocorrelation'') P.show()

El método que escribí a partir de mi última edición ahora es más rápido que incluso scipy.statstools.acf con fft=True hasta que el tamaño de la muestra es muy grande.

Análisis de errores Si desea ajustar los sesgos y obtener estimaciones de error altamente precisas: vea mi código here que implementa este documento de Ulli Wolff ( u original de UW en Matlab )

Funciones probadas

  • a = correlatedData(n=10000) es de una rutina que se encuentra here
  • gamma() es del mismo lugar que correlated_data()
  • acorr() es mi función a continuación
  • estimated_autocorrelation se encuentra en otra respuesta
  • acf() es from statsmodels.tsa.stattools import acf

Tiempos

%timeit a0, junk, junk = gamma(a, f=0) # puwr.py %timeit a1 = [acorr(a, m, i) for i in range(l)] # my own %timeit a2 = acf(a) # statstools %timeit a3 = estimated_autocorrelation(a) # numpy %timeit a4 = acf(a, fft=True) # stats FFT ## -- End pasted text -- 100 loops, best of 3: 7.18 ms per loop 100 loops, best of 3: 2.15 ms per loop 10 loops, best of 3: 88.3 ms per loop 10 loops, best of 3: 87.6 ms per loop 100 loops, best of 3: 3.33 ms per loop

Editar ... Verifiqué nuevamente manteniendo l=40 y cambiando n=10000 a n=200000 muestras. Los métodos de FFT comienzan a obtener un poco de tracción y las implementaciones de statsmodels simplemente lo limitan ... (el orden es el mismo)

## -- End pasted text -- 10 loops, best of 3: 86.2 ms per loop 10 loops, best of 3: 69.5 ms per loop 1 loops, best of 3: 16.2 s per loop 1 loops, best of 3: 16.3 s per loop 10 loops, best of 3: 52.3 ms per loop

Edición 2: cambié mi rutina y volví a probar frente a la FFT para n=10000 n=20000

a = correlatedData(n=200000); b=correlatedData(n=10000) m = a.mean(); rng = np.arange(40); mb = b.mean() %timeit a1 = map(lambda t:acorr(a, m, t), rng) %timeit a1 = map(lambda t:acorr.acorr(b, mb, t), rng) %timeit a4 = acf(a, fft=True) %timeit a4 = acf(b, fft=True) 10 loops, best of 3: 73.3 ms per loop # acorr below 100 loops, best of 3: 2.37 ms per loop # acorr below 10 loops, best of 3: 79.2 ms per loop # statstools with FFT 100 loops, best of 3: 2.69 ms per loop # statstools with FFT

Implementación

def acorr(op_samples, mean, separation, norm = 1): """autocorrelation of a measured operator with optional normalisation the autocorrelation is measured over the 0th axis Required Inputs op_samples :: np.ndarray :: the operator samples mean :: float :: the mean of the operator separation :: int :: the separation between HMC steps norm :: float :: the autocorrelation with separation=0 """ return ((op_samples[:op_samples.size-separation] - mean)*(op_samples[separation:]- mean)).ravel().mean() / norm

4x aceleración se puede lograr a continuación. Debe tener cuidado de solo pasar op_samples=a.copy() ya que modificará la matriz a por a-=mean contrario:

op_samples -= mean return (op_samples[:op_samples.size-separation]*op_samples[separation:]).ravel().mean() / norm

Prueba de cordura

Ejemplo de análisis de error

Esto está un poco fuera de alcance, pero no puedo molestarme en rehacer la figura sin el cálculo del tiempo de autocorrelación integrado o de la ventana de integración. Las autocorrelaciones con errores son claras en la gráfica inferior.

Encontré que esto obtuvo los resultados esperados con solo un ligero cambio:

def estimated_autocorrelation(x): n = len(x) variance = x.var() x = x-x.mean() r = N.correlate(x, x, mode = ''full'') result = r/(variance*n) return result

Pruebas contra los resultados de autocorrelación de Excel.

No creo que haya una función NumPy para este cálculo en particular. Así es como lo escribiría:

def estimated_autocorrelation(x): """ http://.com/q/14297012/190597 http://en.wikipedia.org/wiki/Autocorrelation#Estimation """ n = len(x) variance = x.var() x = x-x.mean() r = np.correlate(x, x, mode = ''full'')[-n:] assert np.allclose(r, np.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)])) result = r/(variance*(np.arange(n, 0, -1))) return result

La declaración de afirmación está ahí para verificar el cálculo y documentar su intención.

Cuando esté seguro de que esta función se comporta como se espera, puede comentar la declaración de afirmación o ejecutar su script con python -O . (La -O le dice a Python que ignore las aseveraciones).

Tomé parte del código de la función pandas autocorrelation_plot (). Comprobé las respuestas con R y los valores coinciden exactamente.

import numpy def acf(series): n = len(series) data = numpy.asarray(series) mean = numpy.mean(data) c0 = numpy.sum((data - mean) ** 2) / float(n) def r(h): acf_lag = ((data[:n - h] - mean) * (data[h:] - mean)).sum() / float(n) / c0 return round(acf_lag, 3) x = numpy.arange(n) # Avoiding lag 0 calculation acf_coeffs = map(r, x) return acf_coeffs