Calculando o PDF de uma forma de onda a partir de suas amostras

27

Há um tempo atrás, eu estava tentando maneiras diferentes de desenhar formas de onda digitais , e uma das coisas que tentei foi, em vez da silhueta padrão do envelope de amplitude, exibi-lo mais como um osciloscópio. É assim que uma onda senoidal e quadrada se parece em um escopo:

insira a descrição da imagem aqui

A maneira ingênua de fazer isso é:

  1. Divida o arquivo de áudio em um pedaço por pixel horizontal na imagem de saída
  2. Calcular o histograma das amplitudes da amostra para cada bloco
  3. Plote o histograma por brilho como uma coluna de pixels

Produz algo como isto: insira a descrição da imagem aqui

Isso funciona bem se houver muitas amostras por bloco e a frequência do sinal não estiver relacionada à frequência de amostragem, mas não o contrário. Se a frequência do sinal for um submúltiplo exato da frequência de amostragem, por exemplo, as amostras sempre ocorrerão exatamente nas mesmas amplitudes em cada ciclo e o histograma será apenas alguns pontos, mesmo que o sinal reconstruído real exista entre esses pontos. Esse pulso senoidal deve ser tão suave quanto o esquerdo acima, mas não é porque é exatamente 1 kHz e as amostras sempre ocorrem nos mesmos pontos:

insira a descrição da imagem aqui

Tentei fazer upsampling para aumentar o número de pontos, mas isso não resolve o problema, apenas ajuda a facilitar as coisas em alguns casos.

Então, o que eu realmente gostaria é uma maneira de calcular o verdadeiro PDF (probabilidade versus amplitude) do sinal reconstruído contínuo de suas amostras digitais (amplitude versus tempo). Não sei qual algoritmo usar para isso. Em geral, o PDF de uma função é a derivada de sua função inversa .

PDF do sin (x):ddxarcsinx=11x2

Mas não sei como calcular isso para ondas em que o inverso é uma função com vários valores , ou como fazê-lo rapidamente. Dividi-lo em galhos e calcular o inverso de cada um, pegar as derivadas e somar todas juntas? Mas isso é bastante complicado e provavelmente existe uma maneira mais simples.

Este "PDF de dados interpolados" também é aplicável a uma tentativa que fiz de estimar a densidade do núcleo de uma trilha GPS. Deveria ter a forma de um anel, mas como ele estava apenas olhando as amostras e não considerando os pontos interpolados entre as amostras, o KDE parecia mais uma corcunda do que um anel. Se as amostras são tudo o que sabemos, é o melhor que podemos fazer. Mas as amostras não são tudo o que sabemos. Também sabemos que existe um caminho entre as amostras. Para o GPS, não existe uma reconstrução Nyquist perfeita como a do áudio com banda ilimitada, mas a idéia básica ainda se aplica, com algumas suposições na função de interpolação.

endólito
fonte
Você tem um exemplo de uma função de vários valores em que está interessado? Você provavelmente terá que avaliá-lo ao longo de um corte de ramificação que faça mais sentido para seus dados físicos.
Lorem Ipsum
Você está mais interessado em maneiras de desenhar esse tipo de plotagem, ou a plotagem é apenas motivação para a pergunta sobre o cálculo do PDF?
datageist
@yoda: Bem, a função acima para onda senoidal é encontrada fazendo apenas meio ciclo, invertendo e obtendo a derivada, porque cada meio ciclo tem o mesmo PDF que o próximo. Mas, para obter o valor de um sinal de áudio arbitrário inteiro, você não pode fazer essa suposição. Eu acho que você precisaria dividi-lo em "cortes de galhos", pegar o PDF de cada um deles e somar todos eles?
endolith
@datageist: Hmm. Estou interessado em maneiras de desenhar esse tipo de enredo, mas esse tipo de enredo é o PDF. Um atalho que produz o mesmo ou muito semelhante resultado está ok.
endolith
@ Endolith, Oh sim, eu entendo. Apenas uma pergunta sobre ênfase realmente (ou seja, que tipos de atalhos são razoáveis).
datageist

Respostas:

7

Interpole para várias vezes a taxa original (por exemplo, 8x sobre-amostragem). Isso permite que você assuma um sinal linear por partes. Este sinal terá muito pouco erro em comparação com a resolução infinita, interpolação contínua de sin (x) / x da forma de onda.

Suponha que cada par de valores superamostrados tenha uma linha contínua de um valor para o próximo. Use todos os valores entre. Isso fornece uma fatia horizontal fina de y1 a y2 a ser acumulada em um PDF de resolução arbitrária. Cada fatia retangular de probabilidade deve ser dimensionada para uma área de 1 / nsamples.

O uso da linha entre amostras em vez da própria amostra impede um PDF "pontudo", mesmo no caso de existir uma relação fundamental entre o período de amostragem e a forma de onda.

Mark Borgerding
fonte
Eu escrevi uma função para o histograma linearmente interpolado, mas é desonesto. Você conhece o código existente para isso?
endolith
A interpolação linear faz uma enorme diferença para a maioria das formas de onda, mesmo sem a super amostragem. O seno de 1 kHz se parece principalmente com o seno de 997 Hz agora. Em vez de apenas linhas horizontais nos valores da amostra, agora são faixas horizontais de cores entre elas. Com a super amostragem, as bandas também são suavizadas. Com a reamostragem da FFT e alguma sobreposição com os pedaços adjacentes, devo conseguir atingir os picos reais entre as amostras. Eu preciso tornar meu código de histograma interpolado mais rápido, ...
endolith 10/10
Eu reescrevi completamente o meu script para isso, e eu acho que eu tenho o histograma e antialiasing certo desta vez: gist.github.com/endolith/652d3ba1a68b629ed328
endolith
A versão mais recente está em github.com/endolith/scopeplot
endolith
7

O que eu diria é essencialmente o "reamostrador aleatório" de Jason R, que por sua vez é uma implementação baseada em sinal pré-amostrado da amostragem estocástica de yoda.

Eu usei interpolação cúbica simples para um ponto aleatório entre cada duas amostras. Para um som de sintetizador primitivo (decaindo de um sinal quadrado saturado e sem banda ilimitada + até harmônicos para um seno), é assim:

PDF sintetizado com reamostragem aleatória

Vamos compará-lo com uma versão de amostra mais alta,

insira a descrição da imagem aqui

e o estranho com a mesma amostra, mas sem interpolação.

insira a descrição da imagem aqui

O artefato notável desse método é o overshoot no domínio quadrado, mas na verdade é assim que o PDF do sinal filtrado por sinc (como eu disse, meu sinal não é ilimitado por banda) também se pareceria e representa a sonoridade percebida muito melhor que os picos, se este fosse um sinal de áudio.

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 é uma lista de variáveis ​​aleatórias no intervalo [0,1].

leftaroundabout
fonte
11
Parece incrível. +1 para o código Haskell.
datageist
Sim, ele deve ultrapassar os valores da amostra. Na verdade, planejei ter um valor de pico para cada coluna de pixel, possivelmente desenhado de maneira diferente, com base nos picos máximos entre amostras e não apenas nas amostras máximas. Formas de onda como flic.kr/p/7QAScX mostram por que isso é necessário.
endolith 12/09
Com "versão com amostra mais alta", você quer dizer que ela é ampliada, mas ainda com amostragem uniforme? E esses são os pontos azuis?
endolith
11
@endolith É simplesmente a forma de onda original calculada em uma taxa de amostragem mais alta em primeiro lugar. Essencialmente, os pontos azuis representam um som amostrado em 192 kHz e os amarelos mais baixos representam uma amostra reduzida ingenuamente para 24 kHz. Os pontos amarelos superiores são stochasticAntiAliasdisso. Mas a versão com maior amostra é, de fato, uma taxa uniforme nos dois casos.
usar o seguinte código
5

Embora sua abordagem seja teoricamente correta (e precise ser levemente modificada para funções não monotônicas), é extremamente difícil calcular o inverso de uma função genérica. Como você diz, terá que lidar com pontos de ramificação e cortes de ramificação, o que é factível, mas você seriamente não gostaria.

Como você já mencionou, a amostragem regular faz a amostragem do mesmo conjunto de pontos e, como tal, é altamente suscetível a estimativas ruins em regiões onde não é amostrada (mesmo que o critério de Nyquist seja satisfeito). Nesse caso, a amostragem por um período mais longo também não ajuda.

Em geral, ao lidar com funções de densidade de probabilidade e histogramas, é uma idéia muito melhor pensar em termos de amostragem estocástica do que a amostragem regular (consulte a resposta vinculada para uma introdução). Ao fazer uma amostragem estocástica, você pode garantir que cada ponto tenha a mesma probabilidade de ser "atingido" e seja uma maneira muito melhor de estimar o pdf.

f(x)=pecado(20πx)+pecado(100πx)fs=1000fN=1001000 amostras (distribuição uniforme) por segundo (não estou usando Hz aqui, porque isso implica um significado diferente) por 30 segundos, fornece o gráfico à direita (mesma classificação).

Você pode ver facilmente que, embora seja barulhento, é uma aproximação muito melhor do PDF real do que o da direita, que mostra zeros em vários intervalos e erros grandes em vários outros. Ao ter um tempo de observação mais longo, é possível reduzir a variação da direita, eventualmente convergindo para o PDF exato (linha preta tracejada) no limite de grandes observações.

insira a descrição da imagem aqui

Lorem Ipsum
fonte
11
"é extremamente difícil calcular o inverso de uma função genérica" ​​Bem, isso não é uma função, mas sim uma série de amostras, portanto, encontrar o inverso é apenas trocar as coordenadas xey das amostras e depois reamostrar para ajustar o novo sistema de coordenadas. Não posso alterar a amostragem de qualquer maneira. Estamos falando de dados pré-existentes criados usando amostragem uniforme.
endolith
4

Estimativa de densidade do kernel

Uma maneira de estimar o PDF de uma forma de onda é usar um estimador de densidade do kernel .

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

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

Atualização: informações adicionais interessantes.

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

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

X(k)eȷ2πnk/N

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

Portanto, adivinhe o que você precisa para reunir todos os PDFs de cada componente Fourier:

|X(k)|1 11 1-x2

X(k)x(n)

Mais pensamento é necessário!

Peter K.
fonte
Pensei nisso, mas a estimativa de densidade é usada para estimar uma função de densidade de probabilidade desconhecida . Por causa do teorema da amostragem de Nyquist, toda a forma de onda é conhecida exatamente e a função exata da densidade de probabilidade também deve ser conhecida. Estou bem com a estimativa se é uma troca de velocidade versus precisão, mas deve haver uma maneira de obter o PDF real dele. Assim, uma forma de onda reconstruída pode ser criada colocando uma função sinc em cada amostra e somando-as. O PDF pode ser criado usando o PDF de uma função sinc como o kernel? Eu não acho que funciona assim.
endolith
Tipo, eu não acho que isso resolva o problema em que as amostras de sinal são submúltiplas da frequência de amostragem. Não leva em consideração a forma de onda reconstruída entre as amostras, leva? Ele apenas desfoca cada ponto do PDF para tentar preencher as lacunas. Eu tive um problema semelhante ao tentar fazer uma estimativa da densidade do kernel de um rastreamento GPS, porque ele não leva em conta os valores entre as amostras.
endolith
4

Como você indicou em um de seus comentários, seria atraente poder calcular o histograma do sinal reconstruído usando apenas as amostras e o PDF da função sinc que interpola sinais ilimitados de banda. Infelizmente, acho que isso não é possível porque o histograma do sinc não possui todas as informações que o sinal em si possui; todas as informações nas posições no domínio do tempo em que cada valor é encontrado são perdidas. Isso torna impossível modelar como as versões em escala e com atraso de tempo do sinc se somariam, o que você desejaria para calcular o histograma da versão "contínua" ou com amostragem ampliada do sinal sem realmente fazer o amostragem ascendente.

Acho que você fica com a interpolação como a melhor opção. Você indicou alguns problemas que o impediram de fazer isso, que eu acho que pode ser resolvido:

  • Despesas computacionais: é claro que sempre é uma preocupação relativa, dependendo do aplicativo específico para o qual você deseja usá-lo. Com base no link que você postou na galeria de renderizações coletadas, suponho que você queira fazer isso para a visualização de sinais de áudio. Se você está interessado nisso para um aplicativo em tempo real ou offline, recomendamos que você protótipo de um interpolador eficiente e veja se ele é realmente muito caro. A reamostragem polifásica é uma boa maneira de fazer isso de maneira flexível (você pode usar qualquer fator racional).

  • π

Jason R
fonte
Mas e se a forma de onda estiver em 44,1 / π kHz? :) Este é um bom conselho, no entanto. Existe algo como reamostragem aleatória? Ou, na verdade, acho que o que funcionaria perfeitamente seria reamostrar de maneira não uniforme, de modo que as novas amostras caibam perfeitamente nos compartimentos na dimensão y, em vez de serem espaçadas uniformemente na dimensão x. Não tenho certeza se existe uma maneira de fazer isso #
1179
2
Você pode implementar facilmente um reamostrador "aleatório" usando uma estrutura Farrow. É um esquema que permite atraso arbitrário de amostra fracionária interpolando usando polinômios (geralmente cúbicos). Você pode manter um acumulador de fase entre amostras, semelhante ao usado em um NCO , que é incrementado por frações pseudo-aleatórias de um intervalo de amostragem para cada amostra de saída (reamostrada). O valor do acumulador é usado como uma entrada para o interpolador Farrow, definindo a quantidade de atraso fracionário para cada saída.
Jason R
Hmm, para esclarecer, o Farrow é apenas uma versão otimizada para processador / memória da interpolação polinomial antiga e regular?
endolith 12/09
11
Sim. É apenas uma estrutura eficiente para implementar o atraso fracionário arbitrário baseado em polinômios.
Jason R
A interpolação cúbica é apenas uma aproximação, no entanto. Quero conhecer verdadeiros picos entre amostras e isso não parece funcionar bem em picos extremos: stackoverflow.com/questions/1851384/… Na verdade, parece que uma série infinita com uma descontinuidade como [..., -1, 1, -1, 1, 1, -1, 1, -1, ...] produzirá um pico infinito entre amostras, portanto, não tenho certeza do quanto isso importaria na prática.
endolith
0

Você precisa suavizar o histograma (isso produzirá resultados semelhantes aos do método kernel). Exatamente como a suavização deve ser realizada, é necessário experimentar. Talvez isso também possa ser feito por interpolação. Além da suavização, acredito que você também obterá melhores resultados se fizermos uma ampliação de sua forma de onda de modo que a frequência de amostragem seja 'significativamente maior' do que a frequência mais alta da sua entrada. Isso deve ajudar no caso "complicado", em que uma onda senoidal está relacionada à frequência de amostragem, de modo que apenas alguns compartimentos no histograma sejam preenchidos. Se levada ao extremo, uma taxa de amostragem suficientemente alta deve fornecer gráficos agradáveis ​​sem suavização. Portanto, a ampliação de amostras combinada com algum tipo de suavização deve gerar melhores plotagens.

Você dá um exemplo de um tom de 1kHz, onde a plotagem não é a esperada. Aqui está a minha proposta (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 o seu tom de 1000Hz, você obtém insira a descrição da imagem aqui

O que você precisa fazer é ajustar a expressão upsampling_factor de acordo com sua preferência.

Ainda não tem 100% de certeza exatamente quais são seus requisitos. Mas, usando o princípio acima de upsampling e suavização, você obtém isso para o tom de 1kHz (feito com o Matlab). Observe que no histograma bruto existem muitos compartimentos com zero de acertos.

insira a descrição da imagem aqui

niaren
fonte
Sim, ele realmente precisa de algum tipo de interpolação como parte do algoritmo. Suavizar o histograma sozinho não serve, porque o histograma possui pontos discretos, não a forma de onda reconstruída. A única maneira de o upsampling funcionar é se eu fizer isso no ponto em que existem muito mais amostras do que pixels verticais, mas esse é um método de força bruta pesada que leva muito tempo.
endolith
ou calculando o efeito de interpolação sobre a saída sem realmente interpolando
endolith