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

27

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:

wprowadź opis zdjęcia tutaj

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: wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

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):ddxarcsinx=11x2

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.

endolit
źródło
Czy masz przykład interesującej Cię funkcji wielowartościowej? Prawdopodobnie będziesz musiał to ocenić wzdłuż cięcia gałęzi, które najlepiej pasuje do twoich danych fizycznych.
Lorem Ipsum
Czy jesteś bardziej zainteresowany sposobami narysowania tego rodzaju fabuły, czy też fabuła jest tylko motywacją do pytania o obliczanie pliku PDF?
Datageist
@yoda: Cóż, powyższą funkcję dla fali sinusoidalnej można znaleźć, biorąc tylko pół cyklu, odwracając i przyjmując pochodną, ​​ponieważ każdy półcykl ma taki sam plik PDF jak następny. Ale aby uzyskać wartość dla całego arbitralnego sygnału audio, nie można przyjąć tego założenia. Myślę, że trzeba by go podzielić na „wycinki gałęzi”, wziąć kolejno każdy plik PDF i zsumować je wszystkie razem?
endolith
@datageist: Hmm. Interesują mnie sposoby narysowania tego rodzaju wykresu, ale taki wykres to PDF. Skrót, który daje taki sam lub bardzo podobny wynik, jest w porządku.
endolith
@endolith, Oh tak, rozumiem. Tylko pytanie o nacisk (naprawdę, jakie rodzaje skrótów są rozsądne).
Datageist

Odpowiedzi:

7

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.

Mark Borgerding
źródło
Napisałem funkcję dla histogramu interpolowanego liniowo, ale jest niejasna. Czy znasz istniejący kod do tego?
endolith,
Interpolacja liniowa stanowi ogromną różnicę dla większości przebiegów, nawet bez nadpróbkowania. Sinus 1 kHz wygląda teraz jak sinus 997 Hz. Zamiast tylko poziomych linii przy wartościach próbki, teraz są między nimi poziome pasy koloru. Dzięki oversamplingowi pasma również są wygładzane. Dzięki resamplingowi FFT i częściowemu nakładaniu się na sąsiednie fragmenty powinienem być w stanie sprawić, by uderzyło ono w prawdziwe szczyty między próbkami. Muszę jednak przyspieszyć interpolowany kod histogramu ...
endolith,
I całkowicie przepisał mój skrypt do tego, i myślę, że mam prawo histogramu i antialiasing ten czas: gist.github.com/endolith/652d3ba1a68b629ed328
endolit
Najnowsza wersja jest w github.com/endolith/scopeplot
endolit
7

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:

PDF z syntezą losową

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

wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

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].

po lewej stronie
źródło
1
Wyglada świetnie. +1 za kod Haskell.
Datageist
Tak, powinien przekraczać wartości próbek. Właściwie planowałem też mieć wartość szczytową dla każdej kolumny pikseli, być może narysowaną inaczej, w oparciu o maksimum pików między próbkami, a nie tylko maksimum próbek. Przebiegi, takie jak flic.kr/p/7QAScX, pokazują, dlaczego jest to konieczne.
endolith,
Pod pojęciem „wersja z wyższą próbką” masz na myśli to, że jest ona próbkowana w sposób upsamplowany, ale nadal próbkowana w sposób jednolity? I to są niebieskie kropki?
endolith,
1
@endolith To po prostu pierwotny kształt fali obliczony na podstawie wyższej częstotliwości próbkowania. Zasadniczo niebieskie punkty oznaczają dźwięk próbkowany z częstotliwością 192 kHz, a dolne żółte reprezentują naiwnie wykonane próbkowanie do 24 kHz. Są to górne żółte punkty stochasticAntiAlias. Ale wersja z wyższą próbką ma w obu przypadkach jednakową stawkę.
leftaroundabout
5

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πx)+sin(100πx)fs=1000fN=1001000

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.

wprowadź opis zdjęcia tutaj

Lorem Ipsum
źródło
1
„niezwykle trudno jest obliczyć odwrotność funkcji ogólnej” Cóż, nie jest to funkcja tyle co seria próbek, więc znalezienie odwrotności polega jedynie na zamianie współrzędnych xiy próbek, a następnie ponownym próbkowaniu w celu dopasowania nowy układ współrzędnych. I tak nie mogę zmienić próbkowania. Mówimy o wcześniej istniejących danych utworzonych przy użyciu jednolitego próbkowania.
endolith,
4

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(x)δ(xx(n))P^

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

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,...,N1X(k)

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

Tak więc jest współczynnikiem :X(k)eȷ2πnk/N

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

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

|X(k)|11x2

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!

Peter K.
źródło
Myślałem o tym, ale estymacja gęstości służy do oszacowania nieznanej funkcji gęstości prawdopodobieństwa. Z powodu twierdzenia Nyquista o próbkowaniu cały przebieg jest dokładnie znany, a dokładna funkcja gęstości prawdopodobieństwa również powinna być możliwa do poznania. Nic mi nie jest, jeśli chodzi o kompromis między szybkością a dokładnością, ale musi istnieć sposób na wyciągnięcie z niego faktycznego pliku PDF. Podobnie zrekonstruowany kształt fali można wykonać, umieszczając funkcję sinc na każdej próbce i sumując je razem. Czy plik PDF można utworzyć za pomocą pliku PDF funkcji sinc jako jądra? Nie sądzę, że tak to działa.
endolith
Nie sądzę, żeby to rozwiązało problem polegający na tym, że próbki sygnałów są podwielokrotnością częstotliwości próbkowania. Nie uwzględnia zrekonstruowanego kształtu fali między próbkami, prawda? Po prostu zamazuje każdy punkt w pliku PDF, aby spróbować wypełnić luki. Miałem podobny problem z próbą oszacowania gęstości jądra śladu GPS, ponieważ nie bierze on pod uwagę wartości między próbkami.
endolith
4

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ą.π

Jason R.
źródło
Ale co jeśli kształt fali wynosi 44,1 / π kHz? :) Ale to dobra rada. Czy istnieje coś takiego jak losowe ponowne próbkowanie? A tak naprawdę, sądzę, że idealnie działałaby próba nierównomiernego przeskalowania, tak aby nowe próbki idealnie pasowały do ​​pojemników w wymiarze y, zamiast równomiernie rozmieszczone w wymiarze x. Nie jestem pewien, czy jest na to sposób
endolith
2
Możesz łatwo wdrożyć „losowy” resampler za pomocą struktury Farrow. Jest to schemat, który pozwala na dowolne opóźnienie próbki ułamkowej poprzez interpolację za pomocą wielomianów (często sześciennych). Można utrzymać akumulator fazy między próbkami, podobny do stosowanego w NCO , który jest zwiększany o pseudolosowe ułamki okresu próbkowania dla każdej wyjściowej (ponownie próbkowanej) próbki. Wartość akumulatora służy jako dane wejściowe do interpolatora Farrow, określając wielkość opóźnienia ułamkowego dla każdego wyjścia.
Jason R
Hmm, dla wyjaśnienia, Farrow jest tylko zoptymalizowaną pod kątem procesorów / pamięci wersją zwykłej starej interpolacji wielomianowej?
endolith,
1
Tak. To po prostu wydajna struktura do implementacji arbitralnego opóźnienia ułamkowego opartego na wielomianach.
Jason R
Interpolacja sześcienna jest jednak tylko przybliżeniem. Chcę poznać prawdziwe szczyty między próbkami i wydaje się, że nie działa to dobrze na ekstremalnych pikach: stackoverflow.com/questions/1851384/... Właściwie wydaje się, że nieskończona seria z nieciągłością, jak [..., -1, 1, -1, 1, 1, -1, 1, -1, ...] wytworzy jednak nieskończony szczyt między próbkami, więc nie jestem pewien, jak bardzo to by miało znaczenie w praktyce.
endolith,
0

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 wprowadź opis zdjęcia tutaj

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.

wprowadź opis zdjęcia tutaj

niaren
źródło
Tak, naprawdę potrzebuje pewnego rodzaju interpolacji jako części algorytmu. Samo wygładzanie histogramu tego nie zrobi, ponieważ histogram składa się z dyskretnych punktów, a nie zrekonstruowanego kształtu fali. Jedynym sposobem, w jaki działałoby próbkowanie w górę, jest zrobienie tego do miejsca, w którym jest o wiele więcej próbek niż w pionowych pikselach, ale jest to metoda wykorzystująca brutalną siłę, która zajmuje dużo czasu.
endolith
lub obliczanie efektu interpolacji na wyjściu bez faktycznej interpolacji
endolith