9
votes

Rééchantillonnage, matrice d'interpolation

J'essaie d'interpoler certaines données dans le but de tracer. Par exemple, étant donné n points de données, j'aimerais pouvoir générer une parcelle "lisse", composée de points de données de 10 * N ou ainsi interpolés.

Mon approche consiste à générer un N-by-10 * N matrice et calculez le produit intérieur Le vecteur d'origine et la matrice que j'ai générée, donnant un vecteur de 1 by-10 * N. J'ai déjà travaillé sur les mathématiques que j'aimerais utiliser pour l'interpolation, mais mon code est assez lent. Je suis assez nouveau à Python, alors j'espère que certains des experts ici peuvent me donner quelques idées de façons d'essayer de accélérer mon code. P>

Je pense qu'une partie du problème est que La génération de la matrice nécessite des appels 10 * N ^ 2 à la fonction suivante: p> xxx pré>

(ceci provient de la théorie de l'échantillonnage . Essentiellement, je tente de recréer un signal de ses échantillons et de la remettre à une fréquence supérieure.) P> La matrice est générée par les éléments suivants: P>

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)


5 commentaires

Complètement mis à part ce que vous essayez de faire avec votre code, l'idée que vous pouvez simplement interpoler des points supplémentaires sans modèle génératif des données n'est fausse. Si vous voulez le faire de manière statistique de manière statistique, vous devez effectuer une sorte de régression. Voir en.wikipedia.org/wiki/generative_model


On dirait que Phil ne veut que utiliser l'interpolation pour le traçage. Tant que les points interpolés ne sont pas utilisés à un autre but, je ne vois pas pourquoi on aurait besoin d'un modèle génératif


@Phil: Une raison particulière pour laquelle vous souhaitez utiliser une interpolation de Sinc, étant donné que c'est un algorithme O (n ^ 2) et d'autres méthodes telles que la spline cubique ne sont que O (n)?


@ Twole18: Le modèle des données est qu'il a été échantillonné selon le FR. wikipedia.org/wiki/nyquist%E2%80%93Shannon_sampling_Theor em . Vous pouvez récupérer l'original exactement en utilisant des fonctions Sinc.


NUMPY a déjà un SINC () Fonction, à la manière. docs.cipy.org/doc/numpy/reference/generated/ numpy.sinc.html


8 Réponses :


0
votes

Petite amélioration. Utilisez la fonction intégrée numpy.sinc (x) qui fonctionne dans le code C compilé.

Amélioration possible: pouvez-vous faire l'interpolation à la volée (comme le tracé survient)? Ou êtes-vous attaché à une bibliothèque de tracé qui accepte uniquement une matrice?


3 commentaires

Merci pour le commentaire. Étrangement, le code exécute environ 10 fois plus lentement lorsque j'ai utilisé Numpy.SIgn (x). Je suis surpris!


Le morceau de traçage de la description était à des fins d'illustration. Je ne suis pas vraiment inquiet de dessiner l'intrigue, de faire le calcul réel plus rapidement. Finalement, cela sera plus d'une tâche «à la volée», car je vais traiter des tranches d'un ensemble de données volumineux. Cependant, comme cela se tient maintenant, courir à travers ce que je considérerais la plus petite tranche de données utile nécessite plus de temps qu'il ne prend que le prochain jeu de données pour arriver ...


TSO = temps d'échantillon initial, TSF = temps d'échantillon final. Donc, si je commence par un signal échantillonné à 1 kHz et que je souhaite générer 10 points interpolés pour chaque échantillon (le nouveau taux d'échantillon sera de 10 kHz), TSO = 0,001, TSF = 0,0001.



1
votes

Votre question n'est pas tout à fait claire; Vous essayez d'optimiser le code que vous avez posté, non?

La ré-écrit Sinc, comme cela devrait accélérer considérablement. Cette implémentation évite de vérifier que le module de mathématiques est importé sur chaque appel, n'accepte pas l'accès à trois reprises et remplace la manipulation des exceptions avec une expression conditionnelle: P>

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)


4 commentaires

"Remplace la manipulation des exceptions avec une expression conditionnelle" mais des exceptions sont plus rapides que les conditionnels de Python. De plus, il serait plus rapide de faire pi * x une fois et l'utiliser deux fois, non?


@endolith Il n'est pas vrai que «les exceptions sont plus rapides que les conditionnels de Python», cela dépend vraiment de la fréquence à laquelle l'état exceptionnel se produit. Quoi qu'il en soit, cela devrait être assez insignifiant ici comparé à l'évitement de la recherche d'import et d'attribut sur chaque appel de fonction. Ne pas utiliser essayer / sauf ici est une question de style et de clarté de code.


@endolith comme pour pi * x , je ne suis pas sûr de créer une nouvelle variable locale pour éviter qu'une multiplication à flotteur unique serait bénéfique. C'est l'une de ces choses que vous devez simplement tester. Encore une fois, cependant, c'est vraiment insignifiant par rapport aux autres changements que j'ai suggéré, ce qui aurait un impact important.


Oui, les exceptions sont plus rapides que les conditionnels, de sorte que s'ils se produisent rarement, les codes les utilisant seront également plus rapides. Dans ce cas, le conditionnel n'arrivera que si l'entrée est exactement 0, ce qui est très rare, il est donc plus rapide d'utiliser l'exception. Dans un test rapide, la version d'exception est d'environ 30% plus rapide pour une entrée aléatoire et utilisant pix = pi * x augmente d'environ 40% environ.



3
votes

Si vous souhaitez interpoler des données de manière assez générale et rapide, des splines ou des polynômes sont très utiles. Scipey a le module sciped.interpoler, qui est très utile. Vous pouvez trouver nombreux exemples dans les pages officielles.


0 commentaires

1
votes

Voici un exemple minimal d'une interpolation 1D avec Scipy - pas autant de plaisir que réinventant, mais.
La parcelle ressemble à sinc code>, qui n'est pas une coïncidence: Essayez Google Spline Rééchantillonner "Sinc approximatif".
(Probablement moins local / plus de robinets ⇒ meilleure approximation, Mais je n'ai aucune idée de la manière dont les plans d'univariats locaux sont.)

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


0 commentaires

1
votes

Je ne suis pas tout à fait sûr de ce que vous essayez de faire, mais il y a des effectifs que vous pouvez faire pour créer la matrice. Suggestion de Braincore à utiliser numpy.sinc est une première étape , mais la seconde consiste à réaliser que les fonctions numpy veulent travailler sur des tableaux numpopiques, où ils peuvent faire des boucles à C Sween, et peut le faire plus rapidement que sur des éléments individuels. XXX

Le truc est qu'en indexant les aranges avec le Numpy.Newaxis, NUMPY convertit le tableau avec une forme I à une avec la forme IX 1 et la matrice avec forme J, pour former 1 x j. À l'étape de soustraction, NUMPY "diffusera" l'entrée de chaque entrée pour agir comme une matrice I x J X et la soustraction. ("Diffusion" est à terme de Numpy, reflétant le fait que aucune copie supplémentaire n'est faite pour étirer l'IX 1 à ix j.)

MAINTENANT LE NUMPY.SIgn peut itérair sur tous les éléments du code compilé, beaucoup plus rapide que tout pour la boucle que vous puissiez écrire.

(il y a une vitesse supplémentaire disponible si vous faites la division avant la soustraction, d'autant plus que la dernière division annule la multiplication.)

Le seul inconvénient est que vous payez maintenant pour un réseau NX10 * N supplémentaire pour tenir la différence. Cela pourrait être un prefreaker si n est grand et la mémoire est un problème.

Sinon, vous devriez pouvoir écrire ceci en utilisant numpy.convolve . D'après ce que je viens d'apprendre à Sinc-Interpolation, je dirais que vous voulez quelque chose comme numpy.convolve (origy.SIgn (numpy.arange (j)), mode = "même") . Mais je me trompe probablement sur les détails.


1 commentaires

Je tente une convolution, alors je pense que Numpy.Convolve pourrait être la bonne direction à prendre.



1
votes

Si votre seul intérêt est de "générer une" parcelle "lisse", j'irais simplement avec une simple courbe de courbe de spline polynomiale:

Pour tous deux points de données adjacents, les coefficients d'une fonction polynomiale de tiers degrés peuvent être calculés à partir des coordonnées de ces points de données et les deux points supplémentaires à gauche et à droite (ignorant les points limites.) Cela générera des points sur une belle courbe lisse avec un premier dirivi continu. Il y a une formule directe pour convertir 4 coordonnées à 4 coefficients polynomiaux, mais je ne veux pas vous priver du plaisir de la regarder; o).


0 commentaires

0
votes

Je vous recommande de vérifier votre algorithme, car c'est un problème non trivial. Spécifiquement, je vous suggère d'accéder à l'article "Fonction tracé à l'aide de Splines coniques" (ICEE Informatique et applications) de Hu et de Pavlidis (1991). Leur implémentation d'algorithme permet un échantillonnage adaptatif de la fonction, de sorte que le temps de rendu est inférieur à celui des approches régulièrement espacées.

L'abstrait suit:

Une méthode est présentée par lequel, donné un Description mathématique d'un fonction, une spline conique approximative L'intrigue de la fonction est produite. Les arcs coniques ont été sélectionnés comme le courbes primitives parce qu'il y a Algorithmes de tracé incrémentiels simples pour les coniques déjà inclus dans certains pilotes de périphérique, et il y a simple algorithmes pour approximations locales par conique. Un algorithme scindé et fusion Pour choisir les nœuds de manière adaptative, Selon l'analyse de forme de la fonction originale basée sur son dérivés de premier ordre, est introduit.


2 commentaires

Mon algorithme provient de la théorie de l'échantillonnage. Essentiellement, je tente de recréer un signal de ses échantillons et de la rétablimer à une fréquence plus élevée. Aux fins de la traçage, je suis sûr que ma solution n'est pas la meilleure méthode ...


@Phil: Vous auriez dû dire que dans la question



9
votes