ricampionamento, matrice interpolante

Question

Sto provando a interpolare alcuni dati allo scopo di tracciare. Ad esempio, dati N punti di dati, mi piacerebbe essere in grado di generare un grafico "liscio", composto da 10 * N o punti di dati interpolati.

Il mio approccio è quello di generare una matrice N-by-10 * N e calcolare il prodotto interno del vettore originale e della matrice I generata, ottenendo un vettore 1-by-10 * N. Ho già elaborato i calcoli che mi piacerebbe utilizzare per l'interpolazione, ma il mio codice è piuttosto lento. Sono abbastanza nuovo in Python, quindi spero che alcuni degli esperti qui possano darmi qualche idea su come posso provare ad accelerare il mio codice.

Credo parte del problema è che la generazione della matrice richiede 10 * N^2 chiamate alla seguente funzione:.

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

(Questo comes from sampling theory Essenzialmente, sto cercando di ricreare un segnale dal suo . campioni, e upsample ad una frequenza superiore)

la matrice viene generato dal seguente:

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) 

sto considerando rompere il compito in a pezzi più piccoli perché non mi piace l'idea di una matrice N^2 seduta nella memoria. Probabilmente potrei fare "resampleMatrix" in una funzione di generatore e fare il prodotto interno riga per riga, ma non credo che questo acceleri il mio codice molto prima di iniziare a fare il paging dentro e fuori dalla memoria.

Grazie in anticipo per i vostri suggerimenti!

Answer 1

Questo è il sovracampionamento. Vedere Help with resampling/upsampling per alcune soluzioni di esempio.

Un modo rapido per eseguire questa operazione (per i dati offline, come l'applicazione di stampa) consiste nell'utilizzare FFT. Questo è ciò che fa nativo di SciPy resample() function. Presuppone tuttavia un segnale periodico, so it's not exactly the same. Vedi this reference:

Ecco il secondo problema riguardante l'interpolazione del segnale reale nel dominio del tempo, ed è davvero un grosso problema. Questo algoritmo di interpolazione esatto fornisce risultati corretti solo se la sequenza x (n) originale è periodica entro il suo intervallo di tempo completo.

La tua funzione presuppone che i campioni del segnale siano tutti 0 all'esterno dell'intervallo definito, quindi i due metodi divergeranno dal punto centrale. Se si riempie il segnale con un sacco di zeri prima, produrrà un risultato molto vicino. Ci sono molti più zeri oltre il bordo della trama non illustrati:

interpolazione cubica non sarà corretto per scopi di ricampionamento. Questo esempio è un caso estremo (vicino alla frequenza di campionamento), ma come potete vedere, l'interpolazione cubica non è nemmeno vicina. Per le frequenze più basse dovrebbe essere abbastanza preciso.

Answer 2

Piccolo miglioramento. Utilizzare la funzione numpy.sinc (x) integrata che viene eseguita nel codice C compilato.

Possibile miglioramento maggiore: è possibile eseguire l'interpolazione al volo (come avviene il disegno)? O sei legato a una libreria di stampa che accetta solo una matrice?

Answer 3

La tua domanda non è del tutto chiara; stai cercando di ottimizzare il codice che hai postato, giusto?

La riscrittura sinc come questa dovrebbe accelerare notevolmente. Questa implementazione consente di evitare la verifica che il modulo per la matematica sia importato a ogni chiamata, non fa attributo di accesso per tre volte, e sostituisce gestione delle eccezioni con un'espressione condizionale:

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

Si potrebbe anche cercare di evitare la creazione di matrice due volte (e tenendolo per due volte in parallelo in memoria) creando direttamente un numpy.array (non da una lista di liste):

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 

(sostituire xrange con portate su Python 3.0 e superiori)

Infine, è possibile crea righe con numpy.arange nonché chiamando numpy.sinc su ogni riga o anche su l'intera matrice:

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) 

Questo dovrebbe essere significativamente più veloce rispetto l'implementazione originale. Prova diverse combinazioni di queste idee e prova le loro prestazioni, guarda quale funziona meglio!

Answer 4

Se si desidera interpolare i dati in modo abbastanza generale e veloce, le spline o i polinomi sono molto utili. Scipy ha il modulo scipy.interpolate, che è molto utile. Puoi trovare many examples nelle pagine ufficiali.

Answer 5

Ecco un esempio minimale di interpolazione 1d con scipy: non è tanto divertente quanto reinventare, ma.
La trama assomiglia a sinc, che non è un caso: prova google spline resample "approssimativo sinc".
(presumibilmente meno locali/più rubinetti ⇒ migliore approssimazione, ma non ho idea di come UnivariateSplines locali sono.)

""" 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() 

Aggiunto feb 2010: si veda anche basic-spline-interpolation-in-a-few-lines-of-numpy

Answer 6

io non sono abbastanza sicuro che cosa stai cercando di fare, ma ci sono alcuni aumenti che puoi fare per creare la matrice. Braincore's suggestion per usare numpy.sinc è un primo passo, ma il secondo è realizzare che le funzioni di numpy vogliono lavorare su array di numpy, dove possono fare loop in C speen e possono farlo più velocemente rispetto ai singoli elementi.

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

Il trucco è che indicizzando i aranges con numpy.newaxis, numpy converte la matrice con forma i per uno con forma i x 1, e l'array di forma j, per sagomare 1 x j. Alla fase di sottrazione, numpy "trasmetterà" ciascun input per agire come un array a forma di i x e fare la sottrazione. ("Broadcast" è il termine di Numpy, che riflette il fatto che non viene eseguita alcuna copia aggiuntiva per estendere ix 1 a ix j.)

Ora numpy.sinc può scorrere su tutti gli elementi nel codice compilato, molto più rapidamente di qualsiasi altro per -solo potresti scrivere

(C'è uno speed-up supplementare disponibile se si fa la divisione prima della sottrazione, soprattutto perché inthe quest'ultima la divisione annulla la moltiplicazione.)

L'unico inconveniente è che ora paga per un extra di Nx10 * N array per contenere la differenza. Questo potrebbe essere un dealbreaker se N è grande e la memoria è un problema.

Altrimenti, dovresti essere in grado di scrivere utilizzando numpy.convolve. Da quel poco che ho imparato sull'interpolazione del sinc, direi che vuoi qualcosa come numpy.convolve(orig,numpy.sinc(numpy.arange(j)),mode="same"). Ma probabilmente ho torto riguardo le specifiche.

Answer 7

Se il vostro unico interesse è quello di 'generare una trama 'liscio'' Vorrei solo andare con un semplice polinomiale curva spline in forma:

Per ogni coppia di punti dati adiacenti i coefficienti di un terzo grado funzione polinomiale può essere calcolato dalle coordinate di quei punti dati e dai due punti addizionali alla loro sinistra e destra (trascurando i punti di confine). Questo genererà punti su una bella curva liscia con un primo dirivibile continuo. C'è una formula semplice per convertire 4 coordinate in 4 coefficienti polinomiali, ma non voglio privarti del divertimento di cercarlo; o).

Answer 8

Si consiglia di controllare l'algoritmo, in quanto si tratta di un problema non banale. Nello specifico, suggerisco di ottenere l'articolo "Funzione di tracciamento usando spline coniche" (IEEE Computer Graphics and Applications) di Hu and Pavlidis (1991). L'implementazione dell'algoritmo consente il campionamento adattativo della funzione, in modo tale che il tempo di rendering sia inferiore rispetto agli approcci con intervalli regolari.

L'estratto segue:

Procedimento viene presentato quale, data una descrizione matematica di una funzione , una spline conica approssima viene prodotto il grafico della funzione. archi coniche sono stati selezionati come curve primitive perché ci sono semplici algoritmi di tracciato incrementali per coniche già incluse in alcuni driver di periferica , e ci sono semplici algoritmi per approssimazioni locali da coniche. Un algoritmo split-and-merge per la scelta adattiva dei nodi, in base all'analisi della forma della funzione originale basata sulle sue derivate del primo ordine , è introdotte.