Berechnen der PDF-Datei einer Wellenform aus ihren Beispielen

27

Vor einiger Zeit habe ich verschiedene Möglichkeiten ausprobiert , um digitale Wellenformen zu zeichnen , und eines der Dinge, die ich versucht habe, war, sie anstelle der Standard-Silhouette der Amplitudenhüllkurve eher wie ein Oszilloskop anzuzeigen. So sehen Sinus und Rechteck auf einem Oszilloskop aus:

Bildbeschreibung hier eingeben

Der naive Weg dazu ist:

  1. Teilen Sie die Audiodatei im Ausgabebild in einen Teil pro horizontalem Pixel auf
  2. Berechnen Sie das Histogramm der Probenamplituden für jeden Block
  3. Zeichnen Sie das Histogramm nach Helligkeit als Pixelspalte

Es erzeugt so etwas: Bildbeschreibung hier eingeben

Dies funktioniert gut, wenn es viele Samples pro Chunk gibt und die Frequenz des Signals nicht mit der Sampling-Frequenz zusammenhängt, aber nicht anders. Wenn die Signalfrequenz beispielsweise ein genaues Submultiple der Abtastfrequenz ist, treten die Abtastungen in jedem Zyklus immer mit genau den gleichen Amplituden auf, und das Histogramm besteht nur aus wenigen Punkten, obwohl das tatsächlich rekonstruierte Signal zwischen diesen Punkten vorhanden ist. Dieser Sinuspuls sollte so glatt sein wie der obige, aber nicht, weil er genau 1 kHz beträgt und die Abtastungen immer um die gleichen Punkte erfolgen:

Bildbeschreibung hier eingeben

Ich habe Upsampling versucht, um die Anzahl der Punkte zu erhöhen, aber es behebt das Problem nicht, sondern hilft nur, die Situation in einigen Fällen zu glätten.

Was ich also wirklich möchte, ist eine Möglichkeit, die wahre PDF (Wahrscheinlichkeit gegen Amplitude) des kontinuierlich rekonstruierten Signals aus seinen digitalen Abtastwerten (Amplitude gegen Zeit) zu berechnen . Ich weiß nicht, welchen Algorithmus ich dafür verwenden soll. Im Allgemeinen ist das PDF einer Funktion die Ableitung ihrer Umkehrfunktion .

PDF von sin (x): ddxarcsinx=11x2

Aber ich weiß nicht, wie ich das für Wellen berechnen soll, bei denen die Inverse eine mehrwertige Funktion ist , oder wie ich es schnell machen soll. Brechen Sie es in Zweige auf und berechnen Sie die Umkehrung von jedem, nehmen Sie die Ableitungen und addieren Sie sie alle zusammen? Aber das ist ziemlich kompliziert und es gibt wahrscheinlich einen einfacheren Weg.

Dieses "PDF mit interpolierten Daten" ist auch auf einen Versuch anwendbar, eine Schätzung der Kerndichte einer GPS-Spur durchzuführen. Es hätte ringförmig sein sollen, aber da es nur die Samples betrachtete und die interpolierten Punkte zwischen den Samples nicht berücksichtigte, sah das KDE eher wie ein Buckel als wie ein Ring aus. Wenn wir nur die Proben kennen, ist dies das Beste, was wir tun können. Aber die Proben sind nicht alles, was wir wissen. Wir wissen auch, dass es einen Pfad zwischen den Samples gibt. Für GPS gibt es keine perfekte Nyquist-Rekonstruktion wie für bandbegrenztes Audio, aber die Grundidee gilt weiterhin, mit einigen Vermutungen in der Interpolationsfunktion.

Endolith
quelle
Haben Sie ein Beispiel für eine mehrwertige Funktion, die Sie interessiert? Sie müssen es wahrscheinlich entlang eines Astschnitts auswerten, der für Ihre physischen Daten am sinnvollsten ist.
Lorem Ipsum
Interessieren Sie sich mehr für Möglichkeiten zum Zeichnen eines solchen Diagramms, oder ist das Diagramm nur eine Motivation für die Frage nach der Berechnung des PDF-Dokuments?
Datageist
@yoda: Nun, die obige Funktion für die Sinuswelle wird gefunden, indem nur ein halber Zyklus genommen, invertiert und die Ableitung genommen wird, da jeder Halbzyklus das gleiche PDF wie der nächste hat. Um den Wert für ein beliebiges Audiosignal zu erhalten, können Sie diese Annahme jedoch nicht treffen. Ich denke, Sie müssten es in "Verzweigungsschnitte" unterteilen, die PDF-Dateien nacheinander abrufen und alle zusammenfassen?
Endolith
@datageist: Hmm. Ich bin an Möglichkeiten interessiert, diese Art von Plot zu zeichnen, aber diese Art von Plot ist das PDF. Eine Verknüpfung, die das gleiche oder ein sehr ähnliches Ergebnis liefert, ist in Ordnung.
Endolith
@ Endolith, Oh ja, ich verstehe. Nur eine Frage zur Betonung (dh welche Arten von Verknüpfungen sind vernünftig).
Datageist

Antworten:

7

Interpolieren Sie auf ein Mehrfaches der ursprünglichen Rate (z. B. 8x überabgetastet). Dadurch können Sie ein stückweise lineares Signal annehmen. Dieses Signal weist verglichen mit der kontinuierlichen sin (x) / x-Interpolation der Wellenform mit unendlicher Auflösung einen sehr geringen Fehler auf.

Angenommen, jedes Paar überabgetasteter Werte weist eine durchgehende Linie von einem Wert zum nächsten auf. Verwenden Sie alle Werte zwischen. Auf diese Weise erhalten Sie eine dünne horizontale Schicht von y1 bis y2, die in einer PDF-Datei mit beliebiger Auflösung akkumuliert werden kann. Jede rechteckige Wahrscheinlichkeitsscheibe muss auf einen Bereich von 1 / nsamples skaliert werden.

Die Verwendung der Linie zwischen Samples anstelle des Samples selbst verhindert ein "spitzes" PDF, selbst wenn eine grundlegende Beziehung zwischen der Sampling-Periode und der Wellenform besteht.

Mark Borgerding
quelle
Ich habe eine Funktion für das linear interpolierte Histogramm geschrieben, aber es ist zweifelhaft. Kennen Sie den vorhandenen Code dafür?
Endolith
Die lineare Interpolation macht für die meisten Wellenformen einen großen Unterschied, auch ohne Überabtastung. 1 kHz-Sinus sieht jetzt meistens wie 997 Hz-Sinus aus. Anstatt nur horizontale Linien bei den Abtastwerten, werden jetzt horizontale Farbstreifen zwischen ihnen angezeigt. Mit Oversampling werden auch die Bänder geglättet. Mit FFT-Resampling und einigen Überlappungen mit benachbarten Chunks sollte es mir möglich sein, die echten Intersample-Peaks zu treffen. Ich muss meinen interpolierten Histogrammcode allerdings schneller machen ...
Endolith
Ich habe mein Skript dafür komplett neu geschrieben und ich glaube, ich habe das Histogramm und das Antialiasing diesmal richtig eingestellt
endolith
Die neueste Version finden Sie unter github.com/endolith/scopeplot
endolith
7

Was ich mitmachen würde, ist im Wesentlichen Jason Rs "Random Resampler", der wiederum eine vorab abgetastete signalbasierte Implementierung von Yodas stochastischem Sampling ist.

Ich habe eine einfache kubische Interpolation auf einen zufälligen Punkt zwischen zwei Samples angewendet. Für einen primitiven Synth-Sound (der von einem gesättigten, nicht bandbegrenzten, quadratischen Signal + geraden Harmonischen zu einem Sinus abfällt) sieht das so aus:

Zufällig neu abgetastetes Synth-PDF

Vergleichen wir es mit einer höher abgetasteten Version,

Bildbeschreibung hier eingeben

und der seltsame mit der gleichen Abtastrate, aber ohne Interpolation.

Bildbeschreibung hier eingeben

Bemerkenswertes Artefakt dieser Methode ist das Überschwingen im quadratischen Bereich, aber so würde auch das PDF des sinc-gefilterten Signals (wie gesagt, mein Signal ist nicht bandbegrenzt) aussehen und die wahrgenommene Lautstärke viel besser wiedergeben als die Spitzen, wenn dies ein Audiosignal wäre.

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 ist eine Liste von Zufallsvariablen im Bereich [0,1].

links herum
quelle
1
Sieht toll aus. +1 für den Haskell-Code.
Datageist
Ja, es sollte die Beispielwerte überschreiten. Eigentlich wollte ich auch für jede Pixelspalte einen Spitzenwert haben, der möglicherweise auf der Grundlage der maximalen Zwischenabtastspitzen und nicht nur der maximalen Abtastwerte unterschiedlich gezeichnet ist. Wellenformen wie flic.kr/p/7QAScX zeigen, warum dies notwendig ist.
Endolith
Mit "höher abgetastete Version" meinen Sie, dass sie hochabgetastet, aber immer noch gleichmäßig abgetastet ist? Und das sind die blauen Punkte?
Endolith
1
@endolith Es ist einfach die ursprüngliche Wellenform, die in erster Linie mit einer höheren Abtastrate berechnet wurde. Im Wesentlichen wie die blauen Punkte einen mit 192 kHz abgetasteten Klang darstellen und die untersten gelben einen auf 24 kHz naiv durchgeführten Downsample. stochasticAntiAliasDavon sind die oberen gelben Punkte . Aber die Version mit der höheren Stichprobe ist in beiden Fällen in der Tat eine einheitliche Rate.
links ungefähr
5

Während Ihr Ansatz theoretisch korrekt ist (und für nicht-monotone Funktionen leicht modifiziert werden muss), ist es äußerst schwierig, die Inverse einer generischen Funktion zu berechnen. Wie Sie sagen, müssen Sie sich mit Verzweigungspunkten und Verzweigungsschnitten befassen, was machbar ist, aber Sie würden es ernsthaft nicht wollen.

Wie Sie bereits erwähnt haben, werden bei der regelmäßigen Stichprobenerhebung die gleichen Punkte erfasst. Daher kann es in Regionen, in denen keine Stichprobenerfassung durchgeführt wird, zu schlechten Schätzungen kommen (auch wenn das Nyquist-Kriterium erfüllt ist). In diesem Fall hilft auch die Probenahme über einen längeren Zeitraum nicht.

Wenn Sie sich mit Wahrscheinlichkeitsdichtefunktionen und Histogrammen befassen, ist es im Allgemeinen eine viel bessere Idee, in Bezug auf stochastische Stichproben zu denken als in Bezug auf reguläre Stichproben (siehe die verknüpfte Antwort für eine Einführung). Durch stochastisches Abtasten können Sie sicherstellen, dass jeder Punkt die gleiche Wahrscheinlichkeit hat, "getroffen" zu werden, und eine viel bessere Methode zum Schätzen des PDFs darstellt.

f(x)=Sünde(20πx)+Sünde(100πx)fs=1000fN=1001000 Samples (gleichmäßige Verteilung) pro Sekunde (ich verwende hier nicht Hz, da dies eine andere Bedeutung impliziert) für 30 Sekunden ergeben die Grafik auf der rechten Seite (gleiches Binning).

Sie können leicht erkennen, dass es, obwohl es verrauscht ist, eine viel bessere Annäherung an die tatsächliche PDF-Datei darstellt als die auf der rechten Seite, die in mehreren Intervallen Nullen und in mehreren anderen große Fehler anzeigt. Wenn Sie eine längere Beobachtungszeit haben, können Sie die Varianz in der rechten verringern und schließlich an der Grenze großer Beobachtungen zur exakten PDF-Datei (gestrichelte schwarze Linie) konvergieren.

Bildbeschreibung hier eingeben

Lorem Ipsum
quelle
1
"Es ist extrem schwierig, die Inverse einer generischen Funktion zu berechnen." Nun, dies ist keine Funktion, sondern eine Reihe von Samples. Wenn Sie also die Inverse finden, tauschen Sie einfach die x- und y-Koordinaten der Samples und passen Sie sie erneut an das neue Koordinatensystem. Ich kann das Sampling sowieso nicht ändern. Es handelt sich um bereits vorhandene Daten, die mithilfe einer einheitlichen Stichprobe erstellt wurden.
Endolith
4

Kernel-Dichteschätzung

Eine Möglichkeit, das PDF einer Wellenform zu schätzen, besteht in der Verwendung eines Kernel-Dichteschätzers .

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

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

Update: Interessante Zusatzinformationen.

x(n)n=0,1,...,N-1X(k)

X(k)=n=0N-1x(n)e-ȷ2πnk/N

X(k)eȷ2πnk/N

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

So eine Vermutung an , was Sie aftermight alle die PDF - Dateien der einzelnen Fourier - Komponente zusammen zu falten sein:

|X(k)|11-x2

X(k)x(n) .

Es sind jedoch weitere Überlegungen erforderlich!

Peter K.
quelle
Ich habe darüber nachgedacht, aber die Dichteschätzung wird zum Schätzen einer unbekannten Wahrscheinlichkeitsdichtefunktion verwendet. Aufgrund des Nyquist-Abtasttheorems ist die gesamte Wellenform genau bekannt, und die genaue Wahrscheinlichkeitsdichtefunktion sollte ebenfalls bekannt sein. Ich bin mit der Schätzung einverstanden, ob es sich um einen Kompromiss zwischen Geschwindigkeit und Genauigkeit handelt, aber es muss eine Möglichkeit geben, das eigentliche PDF herauszuholen. Ebenso kann eine rekonstruierte Wellenform hergestellt werden, indem bei jedem Abtastwert eine Sinusfunktion gesetzt und zusammenaddiert wird. Kann das PDF erstellt werden, indem das PDF einer sinc-Funktion als Kernel verwendet wird? Ich denke nicht, dass es so funktioniert.
Endolith
Ich glaube nicht, dass dies das Problem löst, bei dem die Signalabtastwerte ein Teil der Abtastfrequenz sind. Es berücksichtigt nicht die rekonstruierte Wellenform zwischen Samples, oder? Es verwischt nur jeden Punkt in der PDF, um Lücken zu füllen. Ich hatte ein ähnliches Problem beim Versuch, eine Kernel-Dichteschätzung für eine GPS-Spur durchzuführen, da die Werte zwischen den Abtastwerten nicht berücksichtigt werden.
Endolith
4

Wie Sie in einem Ihrer Kommentare angedeutet haben, wäre es attraktiv, das Histogramm des rekonstruierten Signals nur mit den Samples und dem PDF der sinc-Funktion zu berechnen, die bandbegrenzte Signale interpoliert. Leider glaube ich nicht, dass dies möglich ist, da das Histogramm des Sinc nicht alle Informationen enthält, die das Signal selbst enthält. Alle Informationen zu den Zeitbereichspositionen, an denen jeder Wert angetroffen wird, gehen verloren. Dies macht es unmöglich zu modellieren, wie sich die skalierte und die zeitverzögerte Version des sinc summieren würden, was Sie wünschen würden, um das Histogramm der "kontinuierlichen" oder der aufwärts abgetasteten Version des Signals zu berechnen, ohne dies tatsächlich zu tun Upsampling.

Ich denke, Sie haben die Interpolation als beste Option. Sie haben einige Probleme angedeutet, die Sie daran gehindert haben, dies zu tun. Ich denke, sie können behoben werden:

  • Rechenaufwand: Dies ist natürlich immer ein relatives Problem, abhängig von der spezifischen Anwendung, für die Sie dies verwenden möchten. Basierend auf dem Link, den Sie zur Galerie der von Ihnen gesammelten Renderings gepostet haben, gehe ich davon aus, dass Sie dies für die Visualisierung von Audiosignalen tun möchten. Unabhängig davon, ob Sie sich für eine Echtzeit- oder eine Offline-Anwendung interessieren, möchte ich Sie ermutigen, einen Prototyp eines effizienten Interpolators zu erstellen und zu prüfen, ob dieser wirklich zu kostspielig ist. Polyphasen-Resampling ist eine gute Methode, die flexibel ist (Sie können jeden rationalen Faktor verwenden).

  • π

Jason R
quelle
Was aber, wenn die Wellenform bei 44,1 / π kHz liegt? :) Dies ist jedoch ein guter Rat. Gibt es so etwas wie zufälliges Resampling? Oder wirklich, ich denke, was perfekt funktionieren würde, wäre eine ungleichmäßige Neuabtastung, so dass die neuen Abtastwerte perfekt in Bins in der y-Dimension passen, anstatt in der x-Dimension gleichmäßig beabstandet zu sein. Ich bin mir nicht sicher, ob es eine Möglichkeit gibt, dies zu tun
Endolith
2
Sie können problemlos einen "zufälligen" Resampler mit einer Farrow-Struktur implementieren. Es ist ein Schema, das durch Interpolation mit Polynomen (oft kubisch) eine willkürliche Verzögerung für gebrochene Abtastwerte ermöglicht. Sie könnten einen Inter-Sample-Phasenakkumulator, ähnlich dem in einem NCO verwendeten , beibehalten, der für jedes ausgegebene (neu abgetastete) Sample um pseudozufällige Bruchteile eines Abtastintervalls erhöht wird. Der Wert des Akkumulators wird als Eingabe für den Farrow-Interpolator verwendet und definiert den Betrag der Teilverzögerung für jeden Ausgang.
Jason R
Hmm, um es zu verdeutlichen, ist Farrow nur eine prozessor- / speicheroptimierte Version der regulären alten Polynominterpolation?
Endolith
1
Ja. Es ist nur eine effiziente Struktur zum Implementieren einer auf Polynomen basierenden willkürlichen Bruchverzögerung.
Jason R
Die kubische Interpolation ist jedoch nur eine Annäherung. Ich möchte echte Intersample-Peaks kennen, und es scheint bei extremen Peaks nicht gut zu funktionieren: stackoverflow.com/questions/1851384/… Eigentlich scheint es, dass eine unendliche Reihe mit einer Diskontinuität wie [..., -1, 1, -1, 1, 1, -1, 1, -1, ...] erzeugt jedoch einen unendlichen Intersample-Peak. Ich bin mir also nicht sicher, wie wichtig dies in der Praxis ist.
Endolith
0

Sie müssen das Histogramm glätten (dies führt zu ähnlichen Ergebnissen wie mit einer Kernelmethode). Wie genau die Glättung durchgeführt werden muss, muss experimentiert werden. Möglicherweise könnte dies auch durch Interpolation erfolgen. Neben der Glättung werden Sie meiner Meinung nach auch bessere Ergebnisse erzielen, wenn Sie Ihre Wellenform so upsamplen, dass die Sampling-Frequenz „erheblich höher“ ist als die höchste Frequenz in Ihrem Eingang. Dies sollte in dem "schwierigen" Fall helfen, in dem eine Sinuswelle so mit der Abtastfrequenz zusammenhängt, dass nur wenige Bins im Histogramm gefüllt werden. Im Extremfall sollte eine ausreichend hohe Abtastrate zu schönen Plots ohne Glättung führen. Daher sollte Upsampling in Kombination mit einer Art Glättung zu besseren Plots führen.

Sie geben ein Beispiel für einen 1-kHz-Ton an, bei dem die Darstellung nicht Ihren Erwartungen entspricht. Hier ist mein Vorschlag (Matlab / Octave-Code)

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');

Für Ihren 1000Hz-Ton erhalten Sie dies Bildbeschreibung hier eingeben

Sie müssen lediglich den Ausdruck upsampling_factor nach Ihren Wünschen anpassen.

Immer noch nicht 100% sicher, was genau Ihre Anforderungen sind. Bei Verwendung des oben beschriebenen Prinzips des Upsamplings und der Glättung erhalten Sie dies für den 1-kHz-Ton (mit Matlab erstellt). Beachten Sie, dass das Rohhistogramm viele Bins mit Null Treffern enthält.

Bildbeschreibung hier eingeben

Niaren
quelle
Ja, es braucht wirklich eine Art Interpolation als Teil des Algorithmus. Das Glätten des Histogramms allein reicht nicht aus, da das Histogramm aus diskreten Punkten besteht und nicht aus der rekonstruierten Wellenform. Das Upsampling funktioniert nur, wenn es so weit geht, dass wesentlich mehr Samples als vertikale Pixel vorhanden sind. Dies ist jedoch eine sehr brachiale Methode, die viel Zeit in Anspruch nimmt.
Endolith
oder Berechnen der Auswirkung der Interpolation auf die Ausgabe, ohne tatsächlich zu interpolieren
Endolith