Calcul du PDF d'une forme d'onde à partir de ses échantillons

27

Il y a quelque temps, j'essayais différentes façons de dessiner des formes d'onde numériques , et l'une des choses que j'ai essayées était, au lieu de la silhouette standard de l'enveloppe d'amplitude, de l'afficher plus comme un oscilloscope. Voici à quoi ressemble une onde sinusoïdale et carrée sur une lunette:

entrez la description de l'image ici

La façon naïve de le faire est:

  1. Divisez le fichier audio en un morceau par pixel horizontal dans l'image de sortie
  2. Calculer l'histogramme des amplitudes d'échantillonnage pour chaque bloc
  3. Tracer l'histogramme par luminosité sous forme de colonne de pixels

Cela produit quelque chose comme ceci: entrez la description de l'image ici

Cela fonctionne bien s'il y a beaucoup d'échantillons par morceau et que la fréquence du signal n'est pas liée à la fréquence d'échantillonnage, mais pas autrement. Si la fréquence du signal est un sous-multiple exact de la fréquence d'échantillonnage, par exemple, les échantillons se produiront toujours exactement aux mêmes amplitudes dans chaque cycle et l'histogramme ne sera que de quelques points, même si le signal reconstruit réel existe entre ces points. Cette impulsion sinusoïdale devrait être aussi lisse que ci-dessus à gauche, mais ce n'est pas parce qu'elle est exactement de 1 kHz et que les échantillons se produisent toujours autour des mêmes points:

entrez la description de l'image ici

J'ai essayé le suréchantillonnage pour augmenter le nombre de points, mais cela ne résout pas le problème, aide simplement à lisser les choses dans certains cas.

Donc ce que j'aimerais vraiment, c'est un moyen de calculer le vrai PDF (probabilité vs amplitude) du signal reconstruit continu à partir de ses échantillons numériques (amplitude vs temps). Je ne sais pas quel algorithme utiliser pour cela. En général, le PDF d'une fonction est la dérivée de sa fonction inverse .

PDF de sin (x):ddxarcsinx=11x2

Mais je ne sais pas comment calculer cela pour des vagues où l'inverse est une fonction à valeurs multiples , ni comment le faire rapidement. Le diviser en branches et calculer l'inverse de chacune, prendre les dérivées et les additionner toutes ensemble? Mais c'est assez compliqué et il y a probablement un moyen plus simple.

Ce "PDF de données interpolées" est également applicable à une tentative que j'ai faite de faire une estimation de la densité du noyau d'une trace GPS. Il aurait dû être en forme d'anneau, mais parce qu'il ne regardait que les échantillons et ne tenait pas compte des points interpolés entre les échantillons, le KDE ressemblait plus à une bosse qu'à un anneau. Si les échantillons sont tout ce que nous savons, alors c'est le mieux que nous puissions faire. Mais les échantillons ne sont pas tout ce que nous savons. Nous savons également qu'il existe un chemin entre les échantillons. Pour le GPS, il n'y a pas de reconstruction parfaite de Nyquist comme pour le son à bande limitée, mais l'idée de base s'applique toujours, avec quelques suppositions dans la fonction d'interpolation.

endolith
la source
Avez-vous un exemple de fonction à plusieurs valeurs qui vous intéresse? Vous devrez probablement l'évaluer le long d'une coupe de branche qui a le plus de sens pour vos données physiques.
Lorem Ipsum,
Êtes-vous plus intéressé par des façons de dessiner ce type de tracé, ou le tracé est-il simplement une motivation pour la question sur le calcul du PDF?
datageist
@yoda: Eh bien, la fonction ci-dessus pour l'onde sinusoïdale est trouvée en prenant seulement un demi-cycle, en inversant et en prenant le dérivé, parce que chaque demi-cycle a le même PDF que le suivant. Mais pour obtenir la valeur d'un signal audio arbitraire entier, vous ne pouvez pas faire cette hypothèse. Je pense que vous auriez besoin de le diviser en "coupes de branche", de prendre le PDF de chacun à son tour, et de les additionner tous ensemble?
endolith
@ datatist: Hmm. Je m'intéresse aux façons de dessiner ce genre de tracé, mais ce genre de tracé est le PDF. Un raccourci qui produit le même résultat ou un résultat très similaire est correct.
endolith
@endolith, Oh ouais, je comprends. Juste une question sur l'accent vraiment mis (c'est-à-dire quels types de raccourcis sont raisonnables).
datageist

Réponses:

7

Interpoler à plusieurs fois le taux d'origine (par exemple, 8x suréchantillonné). Cela vous permet de supposer un signal linéaire par morceaux. Ce signal aura très peu d'erreur par rapport à la résolution infinie, interpolation sin (x) / x continue de la forme d'onde.

Supposons que chaque paire de valeurs suréchantillonnées ait une ligne continue d'une valeur à la suivante. Utilisez toutes les valeurs entre. Cela vous donne une fine tranche horizontale de y1 à y2 à accumuler dans un fichier PDF à résolution arbitraire. Chaque tranche rectangulaire de probabilité doit être mise à l'échelle à une zone de 1 / néchantillons.

L'utilisation de la ligne entre les échantillons plutôt que l'échantillon lui-même empêche un PDF "spikey", même dans le cas où il existe une relation fondamentale entre la période d'échantillonnage et la forme d'onde.

Mark Borgerding
la source
J'ai écrit une fonction pour l'histogramme interpolé linéairement, mais c'est douteux. Connaissez-vous le code existant pour cela?
endolith
L'interpolation linéaire fait une énorme différence pour la plupart des formes d'onde, même sans suréchantillonnage. Le sinus de 1 kHz ressemble à présent à 997 Hz. Au lieu de simplement des lignes horizontales aux valeurs d'échantillon, ce sont maintenant des bandes de couleur horizontales entre elles. Avec un suréchantillonnage, les bandes sont également lissées. Avec le rééchantillonnage FFT et un certain chevauchement avec des morceaux adjacents, je devrais être en mesure de le faire atteindre les pics inter-échantillons réels. Je dois accélérer mon code d'histogramme interpolé, cependant ...
endolith
J'ai complètement réécrit mon script pour cela, et je pense que j'ai obtenu l'histogramme et l'antialiasing cette fois: gist.github.com/endolith/652d3ba1a68b629ed328
endolith
La dernière version est disponible sur github.com/endolith/scopeplot
endolith le
7

Ce que j'irais avec est essentiellement le "rééchantillonneur aléatoire" de Jason R, qui à son tour est une implémentation basée sur le signal pré-échantillonné de l'échantillonnage stochastique de Yoda.

J'ai utilisé une simple interpolation cubique à un point aléatoire entre chaque deux échantillons. Pour un son de synthé primitif (décroissance d'un signal carré saturé non limité en bande + même des harmoniques à un sinus), il ressemble à ceci:

Synthétiseur à rééchantillonnage aléatoire PDF

Comparons-le à une version plus échantillonnée,

entrez la description de l'image ici

et l'étrange avec le même échantillonnage mais sans interpolation.

entrez la description de l'image ici

Un artefact notable de cette méthode est le dépassement dans le domaine de type carré, mais c'est en fait ce à quoi ressemblerait le PDF du signal filtré sinc (comme je l'ai dit, mon signal n'est pas limité en bande) et représente bien mieux l'intensité sonore perçue que les pics, s'il s'agissait d'un signal audio.

Code (Haskell):

cubInterpolate vll vl v vr vrr vrrr x
    = v*lSpline x + vr*rSpline x
      + ((vr-vl) - (vrr-vll)/4)*ldSpline x
      + ((vrr-v) - (vrrr-vl)/4)*rdSpline x
     where lSpline x = rSpline (1-x)
           rSpline x = x*x * (3-2*x)
           ldSpline x = x * (1 + x*(x-2))
           rdSpline x = -ldSpline (1-x)

                   --  rand list   IN samples  OUT samples
stochasticAntiAlias :: [Double] -> [Double] -> [Double]
stochasticAntiAlias rs (lsll:lsl:lsc:lsr:lsrr:[]) = []
stochasticAntiAlias (r:rLst) (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)
    = ( cubInterpolate lsll lsl lsc lsr lsrr lsrrr r )
          : stochasticAntiAlias rLst (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)

rand list est une liste de variables aléatoires dans la plage [0,1].

à gauche
la source
1
Ça a l'air génial. +1 pour le code Haskell.
datageist
Oui, il doit dépasser les valeurs d'échantillonnage. En fait, j'avais prévu d'avoir une valeur de pic pour chaque colonne de pixels également, éventuellement dessinée différemment, en fonction du nombre maximal de pics entre échantillons et pas seulement du nombre maximal d'échantillons. Des formes d'onde comme flic.kr/p/7QAScX montrent pourquoi cela est nécessaire.
endolith
Par «version à échantillonnage supérieur», voulez-vous dire qu'elle est suréchantillonnée, mais toujours uniformément échantillonnée? Et ce sont les points bleus?
endolith
1
@endolith C'est tout simplement la forme d'onde d'origine calculée dans un taux d'échantillonnage plus élevé en premier lieu. Comme les points bleus représentent essentiellement un son échantillonné à 192 kHz, et les points jaunes les plus bas représentent un sous-échantillon naïvement fait à 24 kHz. Les points jaunes supérieurs en sont stochasticAntiAlias. Mais la version plus échantillonnée est en effet à taux uniforme dans les deux cas.
leftaroundabout
5

Alors que votre approche est théoriquement correcte (et doit être légèrement modifiée pour les fonctions non monotones), il est extrêmement difficile de calculer l'inverse d'une fonction générique. Comme vous le dites, vous devrez faire face aux points de branchement et aux coupures de branche, ce qui est faisable, mais vous ne le voudriez vraiment pas.

Comme vous l'avez déjà mentionné, l'échantillonnage régulier échantillonne le même ensemble de points et, en tant que tel, est très sensible aux mauvaises estimations dans les régions où il n'échantillonne pas (même si le critère de Nyquist est satisfait). Dans ce cas, l'échantillonnage sur une période plus longue n'aide pas non plus.

En général, lorsqu'il s'agit de fonctions de densité de probabilité et d'histogrammes, c'est une bien meilleure idée de penser en termes d' échantillonnage stochastique que d'échantillonnage régulier (voir la réponse liée pour une introduction). En échantillonnant stochastiquement, vous pouvez vous assurer que chaque point a une probabilité égale d'être «touché» et est une bien meilleure façon d'estimer le pdf.

f(x)=sin(20πx)+sin(100πx)fs=1000fN=1001000 échantillons (distribution uniforme) par seconde (je n'utilise pas Hz ici, car cela implique une signification différente) pendant 30 secondes donne l'intrigue à droite (même binning).

Vous pouvez facilement voir que bien qu'il soit bruyant, c'est une bien meilleure approximation du PDF réel que celui de droite qui montre des zéros à plusieurs intervalles et de grosses erreurs à plusieurs autres. En ayant un temps d'observation plus long, vous pouvez réduire la variance de celui de droite, pour finalement converger vers le PDF exact (ligne noire en pointillés) dans la limite des grandes observations.

entrez la description de l'image ici

Lorem Ipsum
la source
1
"il est extrêmement difficile de calculer l'inverse d'une fonction générique" Eh bien, ce n'est pas une fonction autant qu'une série d'échantillons, donc trouver l'inverse consiste simplement à échanger les coordonnées x et y des échantillons, puis à rééchantillonner pour les ajuster le nouveau système de coordonnées. Je ne peux pas changer l'échantillonnage de toute façon. Nous parlons de données préexistantes créées à l'aide d'un échantillonnage uniforme.
endolith
4

Estimation de la densité du noyau

Une façon d'estimer le PDF d'une forme d'onde est d'utiliser un estimateur de densité de noyau .

x(n)K(x)δ(xx(n))P^

P^(x)=n=0NK(xx(n))

Mise à jour: informations supplémentaires intéressantes.

Supposons que vous ayez votre signal pour , alors --- comme vous dites --- vous pouvez également connaître sa DFT :x(n)n=0,1,...,N1X(k)

X(k)=n=0N1x(n)eȷ2πnk/N

Pour que soit le coefficient de :X(k)eȷ2πnk/N

x(n)=1Nk=0N1X(k)eȷ2πnk/N

Donc, devinez ce que vous pourriez faire après la conversion de tous les PDF de chaque composant Fourier:

|X(k)|11x2

X(k)x(n)

Plus de réflexion cependant!

Peter K.
la source
J'y ai pensé, mais l'estimation de densité est utilisée pour estimer une fonction de densité de probabilité inconnue . En raison du théorème d'échantillonnage de Nyquist, la forme d'onde entière est connue, exactement, et la fonction de densité de probabilité exacte devrait également être connue. Je suis d'accord avec l'estimation si c'est un compromis vitesse / précision, mais il doit y avoir un moyen d'en extraire le PDF réel. Comme, une forme d'onde reconstruite peut être créée en mettant une fonction sinc à chaque échantillon et en les additionnant. Le PDF peut-il être créé en utilisant le PDF d'une fonction sinc comme noyau? Je ne pense pas que cela fonctionne comme ça.
endolith
Comme, je ne pense pas que cela résout le problème où les échantillons de signal sont un sous-multiple de la fréquence d'échantillonnage. Il ne prend pas en compte la forme d'onde reconstruite entre les échantillons, n'est-ce pas? Il brouille simplement chaque point du PDF pour essayer de combler les lacunes. J'ai eu un problème similaire en essayant de faire une estimation de la densité du noyau d'une trace GPS, car elle ne prend pas en compte les valeurs entre les échantillons.
endolith
4

Comme vous l'avez indiqué dans l'un de vos commentaires, il serait intéressant de pouvoir calculer l'histogramme du signal reconstruit en utilisant uniquement les échantillons et le PDF de la fonction sinc qui interpole les signaux à bande limitée. Malheureusement, je ne pense pas que cela soit possible car l'histogramme du sinc ne contient pas toutes les informations que possède le signal lui-même; toutes les informations sur les positions du domaine temporel où chaque valeur est rencontrée sont perdues. Cela rend impossible de modéliser la façon dont les versions échelonnées et temporisées du sinc se résumeraient, ce qui est ce que vous voudriez afin de calculer l'histogramme de la version "continue" ou suréchantillonnée du signal sans réellement faire le suréchantillonnage.

Je pense que vous vous retrouvez avec l'interpolation comme la meilleure option. Vous avez indiqué quelques problèmes qui vous ont empêché de vouloir faire cela, qui, je pense, peuvent être résolus:

  • Dépenses de calcul: il s'agit bien sûr toujours d'une préoccupation relative, en fonction de l'application spécifique pour laquelle vous souhaitez l'utiliser. Sur la base du lien que vous avez publié dans la galerie de rendus que vous avez collectés, je suppose que vous souhaitez le faire pour la visualisation des signaux audio. Que cela vous intéresse pour une application en temps réel ou hors ligne, je vous encourage à créer un prototype d'interpolateur efficace et à voir s'il est vraiment trop coûteux. Le rééchantillonnage polyphasé est un bon moyen de le faire de manière flexible (vous pouvez utiliser n'importe quel facteur rationnel).

  • π

Jason R
la source
Mais que faire si la forme d'onde est à 44,1 / π kHz? :) C'est pourtant un bon conseil. Existe-t-il une chose telle que le rééchantillonnage aléatoire? Ou vraiment, je suppose que ce qui fonctionnerait parfaitement serait de rééchantillonner de manière non uniforme, de sorte que les nouveaux échantillons s'intègrent parfaitement dans les bacs en dimension y, au lieu d'être régulièrement espacés dans la dimension x. Je ne sais pas s'il y a un moyen de le faire
endolith
2
Vous pouvez facilement implémenter un rééchantillonneur "aléatoire" en utilisant une structure Farrow. C'est un schéma qui permet un retard d'échantillonnage fractionnaire arbitraire en interpolant à l'aide de polynômes (souvent cubiques). Vous pouvez conserver un accumulateur de phase inter-échantillons, similaire à celui utilisé dans un NCO , qui est incrémenté par des fractions pseudo-aléatoires d'un intervalle d'échantillonnage pour chaque échantillon de sortie (rééchantillonné). La valeur de l'accumulateur est utilisée comme entrée de l'interpolateur Farrow, définissant la quantité de retard fractionnaire pour chaque sortie.
Jason R
Hmm, pour clarifier, Farrow est juste une version optimisée par processeur / mémoire de l'ancienne interpolation polynomiale régulière?
endolith
1
Oui. C'est juste une structure efficace pour implémenter un retard fractionnaire arbitraire basé sur des polynômes.
Jason R
L'interpolation cubique n'est cependant qu'une approximation. Je veux connaître les vrais pics intersample, et cela ne semble pas bien fonctionner sur les pics extrêmes: stackoverflow.com/questions/1851384/… En fait, il semble qu'une série infinie avec une discontinuité comme [..., -1, 1, -1, 1, 1, -1, 1, -1, ...] produira cependant un pic inter-échantillon infini, donc je ne sais pas dans quelle mesure cela importerait en pratique.
endolith
0

Vous devez lisser l'histogramme (cela donnera des résultats similaires à l'utilisation d'une méthode du noyau). Exactement comment le lissage doit être effectué a besoin d'expérimentation. Peut-être que cela pourrait aussi se faire par interpolation. En plus du lissage, je pense que vous obtiendrez également de meilleurs résultats si vous suréchantillonnez votre forme d'onde de manière à ce que la fréquence d'échantillonnage soit «significativement plus élevée» que la fréquence la plus élevée de votre entrée. Cela devrait aider dans le cas «délicat» où une onde sinusoïdale est liée à la fréquence d'échantillonnage de telle manière que seuls quelques casiers de l'histogramme sont remplis. Si elle est prise à l'extrême, une fréquence d'échantillonnage suffisamment élevée devrait vous donner de belles parcelles sans lissage. Le suréchantillonnage combiné à une sorte de lissage devrait donc donner de meilleurs graphiques.

Vous donnez un exemple de tonalité à 1 kHz, où l'intrigue n'est pas celle que vous attendez. Voici ma proposition (code Matlab / Octave)

pixels_vertical = 100;
% This needs to be tuned to your configuration and acceptance
upsampling_factor = 16*(pixels_vertical/100); 
fs_original = 48000;
fsine = 1000; % in Hz
fs_up = upsampling_factor*fs_original;
duration = 1; % in seconds
x = sin(2*pi*fsine*[0:duration*fs_up]/fs_up);
period_in_samples = fs_up/fsine;
hist_points = linspace(-1,1,pixels_vertical);
istart = 1;
iend   = period_in_samples;
pixel_values = hist(x(istart:iend), hist_points);
% smooth pixel values
[b,a] = butter(2,0.2);
pixel_values_smooth = filtfilt(b,a,pixel_values);
figure;hold on;
plot(hist_points, pixel_values);
plot(hist_points, pixel_values_smooth,'r');

Pour votre ton 1000Hz, vous obtenez ceci entrez la description de l'image ici

Ce que vous devez faire est de régler l'expression upsampling_factor selon vos préférences.

Toujours pas sûr à 100% de vos besoins. Mais en utilisant le principe de suréchantillonnage et de lissage ci-dessus, vous obtiendrez ceci pour la tonalité 1 kHz (faite avec Matlab). Notez que dans l'histogramme brut, il existe de nombreux bacs avec zéro hit.

entrez la description de l'image ici

niaren
la source
Oui, il a vraiment besoin d'un certain type d'interpolation dans le cadre de l'algorithme. Le lissage de l'histogramme seul ne le fera pas, car l'histogramme est constitué de points discrets, pas de la forme d'onde reconstruite. Le seul moyen de suréchantillonnage fonctionnerait si je le fais au point où il y a beaucoup plus d'échantillons que de pixels verticaux, mais c'est une méthode de force brute lourde qui prend beaucoup de temps.
endolith
ou calculer l' effet de l'interpolation sur la sortie sans réellement interpoler
endolith