Obliczanie pliku PDF kształtu fali na podstawie jego próbek

Jakiś czas temu próbowałem różnych sposobów rysowania cyfrowych przebiegów , a jedną z rzeczy, które próbowałem, zamiast standardowej sylwetki obwiedni amplitudy, było wyświetlenie jej bardziej jak oscyloskop. Tak wygląda fala sinusoidalna i prostokątna na lunecie:

Naiwnym sposobem na to jest:

1. Podziel plik audio na jeden fragment na piksel poziomy w obrazie wyjściowym
2. Oblicz histogram amplitud próbek dla każdej porcji
3. Wykreśl histogram według jasności jako kolumny pikseli

Wytwarza coś takiego:

Działa to dobrze, jeśli na próbkę jest dużo próbek, a częstotliwość sygnału nie jest powiązana z częstotliwością próbkowania, ale nie inaczej. Jeśli na przykład częstotliwość sygnału jest dokładnie podwielokrotnością częstotliwości próbkowania, próbki zawsze będą występować z dokładnie tymi samymi amplitudami w każdym cyklu, a histogram będzie tylko kilkoma punktami, nawet jeśli rzeczywisty zrekonstruowany sygnał istnieje między tymi punktami. Ten impuls sinusoidalny powinien być tak gładki, jak powyżej po lewej, ale nie jest tak, ponieważ ma dokładnie 1 kHz, a próbki zawsze pojawiają się wokół tych samych punktów:

Próbowałem upsamplingu, aby zwiększyć liczbę punktów, ale to nie rozwiązuje problemu, po prostu pomaga w niektórych przypadkach wygładzić.

Tak więc naprawdę chciałbym sposób na obliczenie prawdziwego PDF (prawdopodobieństwo vs amplituda) ciągłego zrekonstruowanego sygnału z jego próbek cyfrowych (amplituda w funkcji czasu). Nie wiem, jakiego algorytmu użyć do tego. Zasadniczo plik PDF funkcji jest pochodną jej funkcji odwrotnej .

Plik PDF sin (x):$\frac{d}{dx} \arcsin x = \frac{1}{\sqrt{1-x^2}}$

Ale nie wiem, jak to obliczyć dla fal, gdzie odwrotność jest funkcją wielowartościową , ani jak to zrobić szybko. Podziel go na gałęzie i oblicz odwrotność każdego z nich, weź pochodne i zsumuj je wszystkie razem? Ale to dość skomplikowane i prawdopodobnie istnieje prostszy sposób.

Ten „PDF interpolowanych danych” ma również zastosowanie do podjętej przeze mnie próby oszacowania gęstości jądra ścieżki GPS. Powinien był mieć kształt pierścienia, ale ponieważ patrzył tylko na próbki i nie uwzględniał interpolowanych punktów między próbkami, KDE wyglądało bardziej jak garb niż pierścień. Jeśli próbki są wszystkim, co wiemy, jest to najlepsze, co możemy zrobić. Ale próbki to nie wszystko, co wiemy. Wiemy również, że między próbkami istnieje ścieżka. W przypadku GPS nie ma idealnej rekonstrukcji Nyquista, tak jak w przypadku pasma z ograniczeniem pasma, ale podstawowa idea nadal obowiązuje, z pewną domysłem w funkcji interpolacji.

Interpoluj do wielokrotności pierwotnej stawki (np. 8-krotny oversampling). Umożliwia to przyjęcie częściowego sygnału liniowego. Ten sygnał będzie miał bardzo mały błąd w porównaniu z nieskończoną rozdzielczością, ciągłą interpolacją sin (x) / x kształtu fali.

Załóżmy, że każda para nadpróbkowanych wartości ma linię ciągłą od jednej wartości do następnej. Użyj wszystkich wartości pomiędzy. Daje to jeden cienki poziomy przekrój od y1 do y2, który można zgromadzić w pliku PDF o dowolnej rozdzielczości. Każdy prostokątny wycinek prawdopodobieństwa musi zostać przeskalowany do obszaru 1 / npróbek.

Użycie linii między próbkami zamiast samej próbki zapobiega powstawaniu „kolczastego” pliku PDF, nawet w przypadku, gdy istnieje fundamentalny związek między okresem próbkowania a kształtem fali.

To, co wybrałbym, to w zasadzie „przypadkowy resampler” Jasona R.

Użyłem prostej interpolacji sześciennej do jednego losowego punktu między każdą z dwóch próbek. Dla prymitywnego dźwięku syntezatora (rozpadającego się z nasyconego nie pasmowego sygnału kwadratowego + nawet harmonicznych do sinusoidy) wygląda to tak:

Porównajmy to z wersją o wyższej próbce,

i ten dziwny z tym samym próbkowaniem, ale bez interpolacji.

Znaczącym artefaktem tej metody jest przekroczenie w domenie kwadratowej, ale tak właśnie wyglądałby PDF sygnału filtrowanego przez filtr cynkowy (jak powiedziałem, mój sygnał nie jest ograniczony pasmem) i znacznie lepiej reprezentuje postrzeganą głośność niż szczyty, jeśli byłby to sygnał audio.

Kod (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 to lista zmiennych losowych z zakresu [0,1].

Chociaż twoje podejście jest teoretycznie poprawne (i wymaga nieznacznej modyfikacji w przypadku funkcji niemonotonicznych), niezwykle trudno jest obliczyć odwrotność funkcji ogólnej. Jak mówisz, będziesz musiał radzić sobie z punktami gałęzi i cięciami gałęzi, co jest wykonalne, ale poważnie nie chciałbyś.

Jak już wspomniałeś, regularne próbkowanie próbkuje ten sam zestaw punktów i jako takie jest bardzo podatne na słabe szacunki w regionach, w których nie próbkuje (nawet jeśli kryterium Nyquista jest spełnione). W tym przypadku próbkowanie przez dłuższy czas również nie pomaga.

Zasadniczo, mając do czynienia z funkcjami gęstości prawdopodobieństwa i histogramami, o wiele lepszym pomysłem jest myślenie w kategoriach prób stochastycznych niż regularne pobieranie próbek (patrz wprowadzenie w odpowiedzi na link) Próbując stochastycznie, możesz upewnić się, że każdy punkt ma równe prawdopodobieństwo „trafienia” i jest znacznie lepszym sposobem oszacowania pliku pdf.

Oto przykład: Rozważ funkcję . Teraz, gdy próbuję to z częstotliwością próbkowania Hz (częstotliwość Nyquista, Hz), wynikową gęstością prawdopodobieństwa jest wykres po lewej stronie (401 równych przedziałów między -2 a 2). Nie ma znaczenia, czy próbkuję przez 10 sekund, czy 100. Nadal pozostaje taki sam. Z drugiej strony, próbkowanie stochastycznie z prędkością próbek (rozkład równomierny) na sekundę (nie używam tutaj Hz, ponieważ to sugeruje inne znaczenie) przez 30 sekund daje wykres po prawej stronie (ten sam podział na grupy).$f(x)=\sin(20\pi x) + \sin(100\pi x)$$f_s=1000$$f_N=100$$1000$

Można łatwo zauważyć, że chociaż jest głośny, jest o wiele lepszym przybliżeniem do rzeczywistego pliku PDF niż ten po prawej, który pokazuje zera w kilku odstępach i duże błędy w kilku innych. Mając dłuższy czas obserwacji, możesz zmniejszyć wariancję w tym po prawej stronie, ostatecznie zbliżając się do dokładnego pliku PDF (przerywana czarna linia) na granicy dużych obserwacji.

Szacowanie gęstości jądra

Jednym ze sposobów oszacowania pliku PDF kształtu fali jest użycie estymatora gęstości jądra .

To pobiera wszystkie twoje próbki i funkcję jądra (np. Gaussa), i przekształca je w (delta Diraca) aby oszacować plik PDF, :$x(n)$$K(\mathrm{x})$$\delta(\mathrm{x} - x(n))$$\hat{P}$

Aktualizacja: Interesujące dodatkowe informacje.

Załóżmy, że masz sygnał dla , a następnie --- jak mówisz --- możesz także znać jego DFT :$x(n)$$n= 0,1,...,N-1$$X(k)$

Tak więc jest współczynnikiem :$X(k)$$e^{\jmath 2\pi n k /N}$

Więc zgadnij, co masz zamiar splatać razem wszystkie pliki PDF każdego komponentu Fouriera:

Nie uwzględnia to jednak sposobu, w jaki faza przyczynia się (lub nie) do dodania w .x ( n $X(k)$$x(n)$

Jednak potrzeba więcej myśli!

Jak wskazano w jednym z twoich komentarzy, byłoby atrakcyjne obliczenie histogramu zrekonstruowanego sygnału przy użyciu tylko próbek i pliku PDF funkcji sinc interpolującej sygnały o ograniczonym paśmie. Niestety nie sądzę, aby było to możliwe, ponieważ histogram sinc nie zawiera wszystkich informacji, które zawiera sam sygnał; wszystkie informacje o pozycjach w dziedzinie czasu, w których napotyka się każdą wartość, są tracone. To uniemożliwia modelowanie, w jaki sposób sumowane skalowane i opóźnione w czasie wersje sinc mogłyby się sumować, co jest potrzebne, aby obliczyć histogram „ciągłej” lub próbkowanej w górę wersji sygnału bez faktycznego wykonywania próbkowanie w górę.

Myślę, że najlepszą opcją jest interpolacja. Wskazałeś kilka problemów, które uniemożliwiły ci zrobienie tego, co moim zdaniem można rozwiązać:

• Koszt obliczeniowy: Jest to oczywiście zawsze względny problem, w zależności od konkretnej aplikacji, do której chcesz go użyć. W oparciu o link, który opublikowałeś w galerii renderingów, które zebrałeś, zakładam, że chcesz to zrobić w celu wizualizacji sygnałów audio. Niezależnie od tego, czy jesteś zainteresowany aplikacją w czasie rzeczywistym, czy offline, zachęcam do prototypowania wydajnego interpolatora i zobaczenia, czy jest to naprawdę zbyt kosztowne. Ponowne próbkowanie wielofazowe jest dobrym sposobem na to, aby być elastycznym (można użyć dowolnego racjonalnego czynnika).

• Odchylenie od składników okresowych racjonalnie związanych z częstotliwością próbkowania: Chociaż nie możesz tego całkowicie wyeliminować, powinieneś być w stanie nieco to złagodzić interpolując przez „dziwny” czynnik: zamiast próbkowania w górę o 4, spróbuj 71/18 ( tylko przykład). Będzie to nieco bardziej skomplikowana struktura, ale nadal można ją skutecznie wdrożyć. Zapewni to bardziej równomierny rozkład próbek w okresach komponentów o częstotliwościach związanych z częstotliwością próbkowania. Alternatywnie użyj schematu ponownego próbkowania, który pozwala wybrać dowolny współczynnik ponownego próbkowania, a następnie ponownie próbkować według (w przybliżeniu) liczby niewymiernej, takiej jak . Można to zrobić skutecznie za pomocą interpolatora Farrow , który wykorzystuje interpolację wielomianową.$\pi$

Musisz wygładzić histogram (da to podobne wyniki jak przy użyciu metody jądra). Dokładne wykonanie wygładzania wymaga eksperymentów. Może można to również zrobić przez interpolację. Sądzę, że oprócz wygładzania uzyskasz również lepsze wyniki, jeśli upsamplujesz przebieg w taki sposób, że częstotliwość próbkowania jest „znacznie wyższa” niż najwyższa częstotliwość na wejściu. Powinno to pomóc w „trudnym” przypadku, gdy fala sinusoidalna jest powiązana z częstotliwością próbkowania w taki sposób, że zapełnia się tylko kilka przedziałów na histogramie. Jeśli doprowadzi się do skrajności, wystarczająco wysoka częstotliwość próbkowania powinna dać ładne wykresy bez wygładzania. Tak więc upsampling w połączeniu z pewnego rodzaju wygładzaniem powinien dać lepsze wykresy.

Podajesz przykład tonu 1kHz, w którym wykres nie jest zgodny z oczekiwaniami. Oto moja propozycja (kod 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');

Dostajesz to dla swojego tonu 1000 Hz

To, co musisz zrobić, to dostosować wyrażenie upsampling_factor do swoich preferencji.

Nadal nie jestem w 100% pewien, jakie dokładnie są twoje wymagania. Ale stosując powyższą zasadę upsamplowania i wygładzania, dostaniesz to dla tonu 1kHz (wykonanego za pomocą Matlaba). Zauważ, że na surowym histogramie jest wiele pojemników z zerowymi trafieniami.