remuestreo, matriz de interpolación

Question

Estoy tratando de interpolar algunos datos con el propósito de planear. Por ejemplo, los puntos de datos N dadas, me gustaría ser capaz de generar una trama "suave", compuestas de 10 * N o puntos de datos así interpolados.

Mi enfoque es generar una matriz * N N-por-10 y calcular el producto interno del vector original y la matriz I generado, produciendo un * N vector 1-por-10. Ya he trabajado las matemáticas me gustaría utilizar para la interpolación, pero mi código es bastante lento. Soy bastante nuevo en Python, así que tengo la esperanza de que algunos de los expertos aquí me puede dar algunas ideas de maneras que puede tratar de acelerar mi código.

Creo que parte del problema es que la generación de la matriz requiere 10 * n ^ 2 llamadas a la función siguiente:

def sinc(x):
    import math
    try:
        return math.sin(math.pi * x) / (math.pi * x)
    except ZeroDivisionError:
        return 1.0

( proviene de la teoría del muestreo . En esencia, estoy intentar volver a crear una señal a partir de sus muestras, y aumentará la muestra a una frecuencia más alta.)

La matriz se genera por la siguiente:

def resampleMatrix(Tso, Tsf, o, f):
    from numpy import array as npar
    retval = []

    for i in range(f):
        retval.append([sinc((Tsf*i - Tso*j)/Tso) for j in range(o)])

    return npar(retval)

Estoy pensando en romper la tarea en trozos más pequeños, porque no me gusta la idea de una matriz de N ^ 2 sentados en la memoria. Probablemente podría hacer 'resampleMatrix' en una función de generador y hacer el producto interno fila por fila, pero no creo que va a acelerar mi código mucho hasta que comience la materia paginación dentro y fuera de la memoria.

Gracias de antemano por sus sugerencias!

Answer 1

Este es el muestreo ascendente. Ver Ayuda con remuestreo / muestreo superior para algunos modelos de soluciones.

Una manera rápida de hacer esto (para los datos fuera de línea, al igual que su aplicación trazado) es utilizar FFT. Esto es lo que hace de SciPy nativa resample() función . Supone una señal periódica, sin embargo, así que no es exactamente el mismo . Ver esta referencia :

  

Aquí está la segunda cuestión relativa a la interpolación de dominio de tiempo real de la señal, y es de hecho una gran cosa. Este algoritmo exacto de interpolación proporciona resultados correctos sólo si los x originales (n) secuencia es periódica dentro de su intervalo de tiempo completo.

Su función asume muestras de la señal son todos 0 fuera del rango de tolerancia, por lo que los dos métodos divergen desde el punto central. Si la almohadilla de la señal con una gran cantidad de ceros en primer lugar, se producirá un resultado muy cercano. Hay varios ceros más allá del borde de la trama no se muestra aquí:

interpolación cúbica no será correcta para fines de remuestreo. Este ejemplo es un caso extremo (cerca de la frecuencia de muestreo), pero como se puede ver, la interpolación cúbica no es ni de lejos. Para frecuencias más bajas que debe ser bastante precisa.

Answer 2

Si desea interpolar datos de una manera bastante general y rápido, estrías o polinomios son muy útiles. Scipy tiene el módulo scipy.interpolate, que es muy útil. Puede encontrar muchos ejemplos en las páginas oficiales.

Answer 3

Su pregunta no es del todo claro; usted está tratando de optimizar el código que envió, ¿verdad?

Re-escritura SINC como este debería acelerarlo considerablemente. Esta aplicación evita la comprobación de que el módulo de matemáticas ha sido importada en cada llamada, no hace atribuir el acceso tres veces, y sustituye a la gestión de excepciones con una expresión condicional:

from math import sin, pi
def sinc(x):
    return (sin(pi * x) / (pi * x)) if x != 0 else 1.0

Usted podría también tratar de evitar la creación de la matriz dos veces (y manteniéndolo dos veces en paralelo en memoria) mediante la creación de un numpy.array directamente (no a partir de una lista de listas):

def resampleMatrix(Tso, Tsf, o, f):
    retval = numpy.zeros((f, o))
    for i in xrange(f):
        for j in xrange(o):
            retval[i][j] = sinc((Tsf*i - Tso*j)/Tso)
    return retval

(reemplace xrange con la gama en Python 3.0 y superior)

Por último, puede crear filas con numpy.arange, así como llamar numpy.sinc en cada fila o incluso en toda la matriz:

def resampleMatrix(Tso, Tsf, o, f):
    retval = numpy.zeros((f, o))
    for i in xrange(f):
        retval[i] = numpy.arange(Tsf*i / Tso, Tsf*i / Tso - o, -1.0)
    return numpy.sinc(retval)

Esto debe ser significativamente más rápido que su implementación original. Pruebe diferentes combinaciones de estas ideas y poner a prueba su rendimiento, ver lo que da el mejor!

Answer 4

No estoy muy seguro de lo que estás tratando de hacer, pero hay algunas aceleraciones que puede hacer para crear la matriz. sugerencia de Braincore utilizar numpy.sinc es un primer paso, pero el segundo es darse cuenta de que numpy funciones quieren trabajar en matrices numpy, donde pueden hacer bucles en C speen, y pueden hacerlo más rápido que en los elementos individuales.

def resampleMatrix(Tso, Tsf, o, f):
    retval = numpy.sinc((Tsi*numpy.arange(i)[:,numpy.newaxis]
                         -Tso*numpy.arange(j)[numpy.newaxis,:])/Tso)
    return retval

El truco es que mediante la indexación de los aranges con los numpy.newaxis, numpy convierte la matriz con forma i a uno con forma i x 1, y la matriz con forma de j, para dar forma a 1 x j. En la etapa de sustracción, numpy se "emisión" la cada entrada para actuar como un i x conjunto con forma de j y el hacer la resta. ( "Broadcast" es el término de numpy, lo que refleja el hecho de no se hace ninguna copia adicional para estirar el i x 1 para i x j.)

Ahora el numpy.sinc puede iterar sobre todos los elementos de código compilado, mucho más rápido que cualquier bucle for se podría escribir.

(Hay un adicional de aceleración disponible si lo hace la división antes de la resta, sobre todo porque enel último la división cancela la multiplicación.)

El único inconveniente es que ahora se paga por una serie adicional NX10 * N para mantener la diferencia. Esto podría ser un motivo de ruptura si N es grande y la memoria es un problema.

De lo contrario, usted debe ser capaz de escribir esto utilizando numpy.convolve. Por lo poco que acabo de aprender acerca de sincronismo-interpolación, yo diría que usted quiere algo así como numpy.convolve(orig,numpy.sinc(numpy.arange(j)),mode="same"). Pero probablemente estoy equivocado acerca de los detalles.

Answer 5

Si su único interés es el de 'generar una trama 'suave'' Me gustaría ir con una simple curva spline ajuste polinómico:

Para cualquier par de puntos de datos adyacentes a los coeficientes de una función polinómica de tercer grado se pueden calcular a partir de las coordenadas de los puntos de datos y los dos puntos adicionales a su izquierda y derecha (sin tener en cuenta los puntos limítrofes.) Esto generará puntos en un buen curva suave con una primera dirivitive continua. Hay una fórmula recta hacia adelante para la conversión de coordenadas 4 a 4 coeficientes de los polinomios, pero no quiero privarte de la diversión de mirar hacia arriba; o).

Answer 6

Este es un ejemplo de un mínimo de interpolación 1D con scipy -. No es tan divertido como reinventar, pero
La trama se parece sinc, que no es casualidad: Prueba Google spline volver a muestrear "SINC aproximada".
(Es de suponer menos locales / más grifos ⇒ mejor aproximación, pero no tengo ni idea de cómo son UnivariateSplines locales.)

""" interpolate with scipy.interpolate.UnivariateSpline """
from __future__ import division
import numpy as np
from scipy.interpolate import UnivariateSpline
import pylab as pl

N = 10 
H = 8
x = np.arange(N+1)
xup = np.arange( 0, N, 1/H )
y = np.zeros(N+1);  y[N//2] = 100

interpolator = UnivariateSpline( x, y, k=3, s=0 )  # s=0 interpolates
yup = interpolator( xup )
np.set_printoptions( 1, threshold=100, suppress=True )  # .1f
print "yup:", yup

pl.plot( x, y, "green",  xup, yup, "blue" )
pl.show()

añadido Feb 2010: véase también básica-spline de interpolación-in-a-pocos-líneas-de-numpy

Answer 7

pequeña mejora. Utilice la función incorporada numpy.sinc (x) que se ejecuta en código C compilado.

Posible mejora más grande: ¿Se puede hacer la interpolación sobre la marcha (como se produce el trazado)? ¿O atado a una biblioteca de trazado que sólo acepta una matriz?

Answer 8

Le recomiendo que compruebe su algoritmo, ya que es un problema no trivial. En concreto, le sugiero que accede a la "dibujo de las funciones Utilizando cónicas Splines" artículo (IEEE Computer Graphics and Applications) por Hu y Pavlidis (1991). Su aplicación algoritmo permite para el muestreo adaptativo de la función, de manera que el tiempo de presentación es menor que con los enfoques regularmente espaciados.

El resumen a continuación:

  

Se presenta un método mediante el cual, dada una   Descripción matemática de una   función, un spline cónica se aproxima   la trama de la función se produce.   arcos cónicas fueron seleccionados como la   curvas primitivas porque hay   sencillos algoritmos de trazado incrementales   para las cónicas ya incluido en algunos   controladores de dispositivos, y no son simples   algoritmos para aproximaciones locales por   cónicas. Un algoritmo de división y fusión   para la elección de los nudos de manera adaptativa,   de acuerdo con la forma a análisis de la   función original en función de su   derivados de primer orden, es   introducido.