Calcular el PDF de una forma de onda a partir de sus muestras

27

Hace un tiempo estaba intentando diferentes formas de dibujar formas de onda digitales , y una de las cosas que intenté fue, en lugar de la silueta estándar de la envolvente de amplitud, mostrarla más como un osciloscopio. Así es como se ve una onda sinusoidal y cuadrada en un osciloscopio:

ingrese la descripción de la imagen aquí

La forma ingenua de hacer esto es:

  1. Divida el archivo de audio en un fragmento por píxel horizontal en la imagen de salida
  2. Calcule el histograma de amplitudes de muestra para cada fragmento
  3. Trace el histograma por brillo como una columna de píxeles

Produce algo como esto: ingrese la descripción de la imagen aquí

Esto funciona bien si hay muchas muestras por fragmento y la frecuencia de la señal no está relacionada con la frecuencia de muestreo, pero no de otra manera. Si la frecuencia de la señal es un submúltiplo exacto de la frecuencia de muestreo, por ejemplo, las muestras siempre se producirán exactamente a las mismas amplitudes en cada ciclo y el histograma solo tendrá unos pocos puntos, aunque la señal reconstruida real exista entre estos puntos. Este pulso sinusoidal debe ser tan suave como el anterior a la izquierda, pero no es porque sea exactamente 1 kHz y las muestras siempre ocurren alrededor de los mismos puntos:

ingrese la descripción de la imagen aquí

Intenté el muestreo para aumentar el número de puntos, pero no resuelve el problema, solo ayuda a suavizar las cosas en algunos casos.

Entonces, lo que realmente me gustaría es una forma de calcular el PDF verdadero (probabilidad vs amplitud) de la señal reconstruida continua a partir de sus muestras digitales (amplitud vs tiempo). No sé qué algoritmo usar para esto. En general, el PDF de una función es la derivada de su función inversa .

PDF de sin (x):rereXarcsinX=11-X2

Pero no sé cómo calcular esto para ondas donde la inversa es una función de valores múltiples , o cómo hacerlo rápido. Divídalo en ramas y calcule el inverso de cada uno, tome las derivadas y sume todas. Pero eso es bastante complicado y probablemente haya una forma más simple.

Este "PDF de datos interpolados" también es aplicable a un intento que hice para hacer una estimación de la densidad del núcleo de una pista GPS. Debería haber sido en forma de anillo, pero debido a que solo miraba las muestras y no consideraba los puntos interpolados entre las muestras, el KDE parecía más una joroba que un anillo. Si las muestras son todo lo que sabemos, entonces esto es lo mejor que podemos hacer. Pero las muestras no son todo lo que sabemos. También sabemos que hay un camino entre las muestras. Para GPS, no hay una reconstrucción perfecta de Nyquist como la hay para audio de banda ilimitada, pero la idea básica aún se aplica, con algunas conjeturas en la función de interpolación.

endolito
fuente
¿Tiene un ejemplo de una función multivalor que le interese? Probablemente tendrá que evaluarlo a lo largo de un corte de rama que tenga más sentido para sus datos físicos.
Lorem Ipsum
¿Está más interesado en formas de dibujar ese tipo de trama, o la trama es solo una motivación para la pregunta sobre el cálculo del PDF?
Datageist
@yoda: Bueno, la función anterior para la onda sinusoidal se encuentra tomando solo medio ciclo, invirtiendo y tomando la derivada, porque cada medio ciclo tiene el mismo PDF que el siguiente. Pero para obtener el valor de una señal de audio arbitraria completa, no puede hacer esa suposición. ¿Creo que necesitaría dividirlo en "cortes de rama", tomar el PDF de cada uno y sumarlos todos juntos?
Endolith
@datageist: Hmm. Estoy interesado en formas de dibujar ese tipo de trama, pero ese tipo de trama es el PDF. Un acceso directo que produce el mismo resultado o muy similar está bien.
Endolith
@endolith, Oh sí, lo entiendo. Solo una pregunta sobre el énfasis realmente (es decir, qué tipos de atajos son razonables).
Datageist

Respuestas:

7

Interpolar varias veces la tasa original (por ejemplo, 8x sobremuestreado). Esto le permite asumir una señal lineal por partes. Esta señal tendrá muy poco error en comparación con la resolución infinita, interpolación sin (x) / x continua de la forma de onda.

Suponga que cada par de valores sobremuestreados tiene una línea continua de un valor al siguiente. Use todos los valores entre. Esto le proporciona un corte horizontal delgado de y1 a y2 para que se acumule en un PDF de resolución arbitraria. Cada segmento rectangular de probabilidad debe escalarse a un área de muestras 1 / n.

El uso de la línea entre muestras en lugar de la muestra misma evita un PDF "puntiagudo", incluso en el caso de que haya una relación fundamental entre el período de muestreo y la forma de onda.

Mark Borgerding
fuente
He escrito una función para el histograma interpolado linealmente, pero es dudoso. ¿Conoces el código existente para esto?
endolito
La interpolación lineal hace una gran diferencia para la mayoría de las formas de onda, incluso sin el sobremuestreo. El seno de 1 kHz se parece principalmente al seno de 997 Hz ahora. En lugar de solo líneas horizontales en los valores de muestra, ahora son bandas horizontales de color entre ellas. Con el sobremuestreo, las bandas también se suavizan. Con el remuestreo FFT y cierta superposición con fragmentos adyacentes, debería ser capaz de hacer que llegue a los picos reales entre muestras. Aunque necesito hacer mi código de histograma interpolado más rápido ...
endolito
Reescribí completamente mi script para esto, y creo que obtuve el histograma y el antialiasing en este momento: gist.github.com/endolith/652d3ba1a68b629ed328
endolith
La última versión está en github.com/endolith/scopeplot
endolith
7

Con lo que me gustaría ir es esencialmente el "remuestreador aleatorio" de Jason R, que a su vez es una implementación basada en señales previamente muestreadas del muestreo estocástico de Yoda.

He usado interpolación cúbica simple a un punto aleatorio entre cada dos muestras. Para un sonido de sintetizador primitivo (descomposición de una señal cuadrada saturada sin límite de banda + incluso armónicos a un seno) se ve así:

PDF de sintetizador con muestreo aleatorio

Comparémoslo con una versión de muestra más alta,

ingrese la descripción de la imagen aquí

y el extraño con la misma frecuencia de muestreo pero sin interpolación.

ingrese la descripción de la imagen aquí

El artefacto notable de este método es el sobreimpulso en el dominio cuadrado, pero en realidad es así como se vería el PDF de la señal filtrada sinc (como dije, mi señal no tiene límite de banda) y representa mucho mejor el volumen percibido. que los picos, si esto fuera una señal de audio.

Código (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 es una lista de variables aleatorias en el rango [0,1].

a la izquierda
fuente
1
Parece increíble. +1 para el código Haskell.
Datageist
Sí, debería sobrepasar los valores de la muestra. De hecho, también planeé tener un valor máximo para cada columna de píxeles, posiblemente dibujado de manera diferente, basado en los picos máximos entre muestras y no solo en las muestras máximas. Las formas de onda como flic.kr/p/7QAScX muestran por qué esto es necesario.
Endolith
Por "versión de muestra superior", ¿quiere decir que está muestreado en exceso, pero que todavía está muestreado de manera uniforme? ¿Y esos son los puntos azules?
Endolith
1
@endolith Es simplemente la forma de onda original calculada en una frecuencia de muestreo más alta en primer lugar. Esencialmente, como los puntos azules representan un sonido muestreado a 192 kHz, y los amarillos más bajos representan una disminución de la muestra ingenuamente hecha a 24 kHz. Los puntos amarillos superiores son stochasticAntiAliasde esto. Pero la versión de muestra más alta es de hecho una tasa uniforme en ambos casos.
Leftaroundabout
5

Si bien su enfoque es teóricamente correcto (y debe modificarse ligeramente para funciones no monótonas), es extremadamente difícil calcular el inverso de una función genérica. Como dices, tendrás que lidiar con los puntos de ramificación y los cortes de ramificación, lo cual es factible, pero en serio no querrás hacerlo.

Como ya mencionó, el muestreo regular muestrea el mismo conjunto de puntos y, como tal, es altamente susceptible a estimaciones pobres en regiones donde no se muestrea (incluso si se cumple el criterio de Nyquist). En este caso, el muestreo por un período más largo tampoco ayuda.

En general, cuando se trata de funciones de densidad de probabilidad e histogramas, es una idea mucho mejor pensar en términos de muestreo estocástico que el muestreo regular (vea la respuesta vinculada para una introducción). Al tomar muestras estocásticamente, puede asegurarse de que cada punto tenga la misma probabilidad de ser "alcanzado" y que sea una forma mucho mejor de estimar el pdf.

F(X)=pecado(20πX)+pecado(100πX)Fs=1000Fnorte=1001000 las muestras (distribución uniforme) por segundo (no estoy usando Hz aquí, porque eso implica un significado diferente) durante 30 segundos da la gráfica a la derecha (mismo binning).

Puede ver fácilmente que, aunque es ruidoso, es una aproximación mucho mejor al PDF real que el de la derecha que muestra ceros en varios intervalos y grandes errores en varios otros. Al tener un tiempo de observación más largo, puede reducir la varianza en el de la derecha, eventualmente convergiendo al PDF exacto (línea negra discontinua) en el límite de observaciones grandes.

ingrese la descripción de la imagen aquí

Lorem Ipsum
fuente
1
"es extremadamente difícil calcular el inverso de una función genérica" ​​Bueno, esta no es una función sino una serie de muestras, por lo que encontrar el inverso es simplemente intercambiar las coordenadas x e y de las muestras y luego volver a muestrear para ajustar El nuevo sistema de coordenadas. No puedo cambiar el muestreo de todos modos. Estamos hablando de datos preexistentes creados con muestreo uniforme.
Endolith
4

Estimación de densidad del núcleo

Una forma de estimar el PDF de una forma de onda es usar un estimador de densidad del núcleo .

X(norte)K(X)δ(X-X(norte))PAGS^

PAGS^(X)=norte=0 0norteK(X-X(norte))

Actualización: información adicional interesante.

X(norte)norte=0 0,1,...,norte-1X(k)

X(k)=norte=0 0norte-1X(norte)mi-ȷ2πnortek/ /norte

X(k)miȷ2πnortek/ /norte

X(norte)=1nortek=0 0norte-1X(k)miȷ2πnortek/ /norte

Por lo tanto, adivine qué es lo que puede ser para involucrar todos los archivos PDF de cada componente de Fourier:

El |X(k)El |11-X2

X(k)X(norte)

¡Sin embargo, se requiere más pensamiento!

Peter K.
fuente
Pensé en eso, pero la estimación de densidad se usa para estimar una función de densidad de probabilidad desconocida . Debido al teorema de muestreo de Nyquist, toda la forma de onda se conoce exactamente, y la función de densidad de probabilidad exacta también debe ser conocida. Estoy de acuerdo con la estimación si es una compensación entre velocidad y precisión, pero debe haber una forma de obtener el PDF real. Como, una forma de onda reconstruida se puede hacer al poner una función sinc en cada muestra y sumarlas. ¿Se puede crear el PDF utilizando el PDF de una función sinc como núcleo? No creo que funcione así.
Endolith
Como, no creo que esto resuelva el problema donde las muestras de señal son un submúltiplo de la frecuencia de muestreo. No tiene en cuenta la forma de onda reconstruida entre muestras, ¿verdad? Simplemente difumina cada punto en el PDF para tratar de llenar los vacíos. Tuve un problema similar al intentar hacer una estimación de la densidad del núcleo de un rastreo GPS, porque no tiene en cuenta los valores entre muestras.
Endolith
4

Como indicó en uno de sus comentarios, sería atractivo poder calcular el histograma de la señal reconstruida utilizando solo las muestras y el PDF de la función sinc que interpola las señales de límite de banda. Desafortunadamente, no creo que esto sea posible porque el histograma del sinc no tiene toda la información que tiene la señal; toda la información sobre las posiciones en el dominio del tiempo donde se encuentra cada valor se pierde. Esto hace que sea imposible modelar cómo se sumarían las versiones escaladas y retardadas de tiempo de la sinc, que es lo que desearía para calcular el histograma de la versión "continua" o muestreada hacia arriba de la señal sin hacer realmente el muestreo ascendente.

Creo que te queda la interpolación como la mejor opción. Indicaste un par de problemas que te impidieron querer hacer esto, lo que creo que puede abordarse:

  • Gastos computacionales: Por supuesto, esto es siempre una preocupación relativa, dependiendo de la aplicación específica para la que desea usar esto. Según el enlace que publicó en la galería de representaciones que ha recopilado, supongo que desea hacer esto para visualizar las señales de audio. Ya sea que esté interesado en esto para una aplicación en tiempo real o fuera de línea, le animo a que realice un prototipo de un interpolador eficiente y vea si es realmente demasiado costoso. El remuestreo polifásico es una buena manera de hacer esto que es flexible (puede usar cualquier factor racional).

  • π

Jason R
fuente
Pero, ¿qué pasa si la forma de onda está a 44.1 / π kHz? :) Sin embargo, este es un buen consejo. ¿Existe algo como el muestreo aleatorio? O realmente, supongo que lo que funcionaría perfectamente sería volver a muestrear de manera no uniforme, de modo que las nuevas muestras encajen perfectamente en contenedores en la dimensión y, en lugar de estar espaciados uniformemente en la dimensión x. No estoy seguro de si hay una manera de hacerlo
endolito
2
Puede implementar fácilmente un muestreador "aleatorio" utilizando una estructura Farrow. Es un esquema que permite un retraso arbitrario de la muestra fraccional al interpolar utilizando polinomios (a menudo cúbicos). Podría mantener un acumulador de fase entre muestras, similar al utilizado en un NCO , que se incrementa mediante fracciones pseudoaleatorias de un intervalo de muestreo para cada muestra de salida (remuestreada). El valor del acumulador se utiliza como entrada para el interpolador Farrow, que define la cantidad de retraso fraccional para cada salida.
Jason R
Hmm, para aclarar, ¿Farrow es solo una versión optimizada de procesador / memoria de la vieja interpolación polinómica normal?
Endolith
1
Sí. Es solo una estructura eficiente para implementar retrasos fraccionarios arbitrarios basados ​​en polinomios.
Jason R
Sin embargo, la interpolación cúbica es solo una aproximación. Quiero saber verdaderos picos entre muestras, y no parece funcionar bien en picos extremos: stackoverflow.com/questions/1851384/… En realidad, parece que una serie infinita con una discontinuidad como [..., -1, Sin embargo, 1, -1, 1, 1, -1, 1, -1, ..., producirán un pico entre muestras infinito, así que no estoy seguro de cuánto importará esto en la práctica.
endolito el
0

Necesita suavizar el histograma (esto arrojará resultados similares a los de usar un método kernel). Exactamente cómo se debe realizar el alisado necesita experimentación. Tal vez también podría hacerse por interpolación. Además del suavizado, creo que también obtendrá mejores resultados si sube la muestra de su forma de onda de tal manera que la frecuencia de muestreo sea 'significativamente mayor' que la frecuencia más alta en su entrada. Esto debería ayudar en el caso 'complicado' en el que una onda sinusoidal está relacionada con la frecuencia de muestreo de tal manera que solo se rellenan unos pocos contenedores en el histograma. Si se lleva al extremo, una frecuencia de muestreo lo suficientemente alta debería proporcionar buenas parcelas sin suavizar. Por lo tanto, el muestreo ascendente combinado con algún tipo de suavizado debería producir mejores parcelas.

Usted da un ejemplo de un tono de 1 kHz, donde la trama no es la esperada. Aquí está mi propuesta (código 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');

Para tu tono de 1000Hz obtienes esto ingrese la descripción de la imagen aquí

Lo que debe hacer es ajustar la expresión upsampling_factor a su preferencia.

Todavía no estoy 100% seguro exactamente cuáles son sus requisitos. Pero utilizando el principio anterior de muestreo y suavizado obtendrá esto para el tono de 1 kHz (hecho con Matlab). Tenga en cuenta que en el histograma sin formato hay muchos contenedores con cero aciertos.

ingrese la descripción de la imagen aquí

niaren
fuente
Sí, realmente necesita algún tipo de interpolación como parte del algoritmo. Alisar el histograma solo no lo hará, porque el histograma es de puntos discretos, no la forma de onda reconstruida. La única forma en que el muestreo ascendente funcionaría es si lo hago hasta el punto en que haya muchas más muestras que píxeles verticales, pero ese es un método de fuerza bruta fuerte que lleva mucho tiempo.
Endolith
o calcular el efecto de interpolar en la salida sin realmente interpolar
endolith