Jak wdrożyć korelację krzyżową, aby udowodnić, że dwa pliki audio są podobne?

58

Muszę wykonać korelację krzyżową dwóch plików audio, aby udowodnić, że są one podobne. Wziąłem FFT dwóch plików audio i mam ich wartości widma mocy w osobnych tablicach.

Jak powinienem kontynuować ich korelację krzyżową i udowodnić, że są do siebie podobne? Czy jest na to lepszy sposób? Wszelkie podstawowe pomysły pomogą mi się nauczyć i zastosować.

— Lorem Ipsum
źródło

Biorąc pod uwagę korelację krzyżową dwóch losowych wektorów sygnałowych. Jak zaimplementować odwrotność, aby uzyskać dwa wektory w MATLAB. John Muhehe

56

Korelacja krzyżowa i splot są ze sobą ściśle powiązane. Krótko mówiąc, aby przeprowadzić splot z FFT, ty

zeruj piki sygnałów wejściowych (dodaj zera na końcu, aby co najmniej połowa fali była „pusta”)
weź FFT obu sygnałów
pomnóż wyniki razem (mnożenie elementarne)
wykonaj odwrotną FFT

conv(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros))

Musisz uzupełnić zero, ponieważ metoda FFT jest w rzeczywistości korelacją krzyżową kołową , co oznacza, że sygnał zawija się na końcach. Więc dodajesz wystarczającą liczbę zer, aby pozbyć się nakładania, aby zasymulować sygnał od zera do nieskończoności.

Aby uzyskać korelację krzyżową zamiast splotu, musisz albo odwrócić czas jednego z sygnałów przed wykonaniem FFT, albo wziąć złożoną koniugat jednego z sygnałów po FFT:

corr(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros[reversed]))
corr(a, b) = ifft(fft(a_and_zeros) * conj(fft(b_and_zeros)))

cokolwiek jest łatwiejsze ze sprzętem / oprogramowaniem. W przypadku autokorelacji (korelacji krzyżowej sygnału z samym sobą) lepiej jest wykonać koniugat złożony, ponieważ wtedy wystarczy obliczyć FFT tylko raz.

Jeśli sygnały są prawdziwe, możesz użyć prawdziwych FFT (RFFT / IRFFT) i zaoszczędzić połowę czasu obliczeniowego, obliczając tylko połowę widma.

Możesz także zaoszczędzić czas obliczeń, wypełniając do większego rozmiaru, dla którego FFT jest zoptymalizowany (takiego jak 5-gładka liczba dla FFTPACK, ~ 13-gładka liczba dla FFTW lub potęga 2 dla prostej implementacji sprzętowej).

Oto przykład w Pythonie korelacji FFT w porównaniu z korelacją siłową: https://stackoverflow.com/a/1768140/125507

To da ci funkcję korelacji krzyżowej, która jest miarą podobieństwa do przesunięcia. Aby uzyskać przesunięcie, w którym fale są „wyrównane” ze sobą, w funkcji korelacji wystąpi szczyt:

Wartość x piku jest przesunięciem, które może być ujemne lub dodatnie.

Widziałem to tylko w celu znalezienia przesunięcia między dwiema falami. Można uzyskać dokładniejsze oszacowanie przesunięcia (lepsze niż rozdzielczość próbek), stosując interpolację paraboliczną / kwadratową na szczycie.

Aby uzyskać wartość podobieństwa między -1 a 1 (wartość ujemna wskazująca, że jeden z sygnałów maleje wraz ze wzrostem drugiego), należy skalować amplitudę zgodnie z długością wejść, długością FFT, konkretną implementacją FFT skalowanie itp. Autokorelacja fali z samym sobą daje wartość maksymalnego możliwego dopasowania.

Pamiętaj, że będzie to działać tylko na falach o tym samym kształcie. Jeśli próbkowano je na innym sprzęcie lub dodano szumy, ale poza tym nadal mają ten sam kształt, to porównanie będzie działać, ale jeśli kształt fali został zmieniony przez filtrowanie lub przesunięcia fazowe, mogą brzmieć tak samo, ale wygrał też nie są skorelowane.

— endolit
źródło

3

Wypełnienie zerowe powinno wynosić co najmniej N = rozmiar (a) + rozmiar (b) -1, najlepiej zaokrąglony w górę do potęgi 2. Aby uzyskać wartość między -1 a 1, podziel przez normę (a) * normę (b ), który podaje cosinus kąta między dwoma wektorami w przestrzeni N dla danego opóźnienia (tj. moduł przesunięcia kołowego N). W skrajnych opóźnieniach nie ma wielu nakładających się próbek (tylko jedna w skrajnie skrajnej), więc podzielenie przez normę (a) * normę (b) spowoduje przesunięcie tych korelacji w kierunku 0 (tj. Pokazanie ich względnej ortogonalności w przestrzeni N) .

— Eryk Sun

1

Myślę, że może być błąd w opisie. Czy mnożenie FFT razem termin po terminie nie powinno dawać FFT splotu sygnałów, a nie FFT korelacji krzyżowej ? Jak rozumiem, aby uzyskać FFT korelacji krzyżowej, konieczne jest zastosowanie złożonej koniugatu jednego z wektorów FFT w multiplikacjach termin po semestrze przed podjęciem iFFT.

— Dilip Sarwate

@DilipSarwate: Tak, masz rację. Możesz także odwrócić jeden sygnał w kierunku czasu, który dodałem do odpowiedzi.

— endolith

1

„Dlaczego odwrócenie czasu jest trudne w sprzęcie?” W wielu przypadkach dane są przechowywane w układach skurczowych w oczekiwaniu, że obliczenia są lokalne , tj. , przechowywane w tej komórce, oddziałuje tylko z najbliższymi sąsiadami . Wysyłanie do komórki # i wysyłanie do komórki # , i robienie tego dla wszystkich zwiększenie kosztów okablowania, opóźnień okablowania (a tym samym zmniejsza maksymalne osiągalne taktowanie zegara), a także dlatego, że wszystkie przewody muszą się przecinać, co powoduje problemy z routingiem. Należy unikać, jeśli to możliwe, w tym przypadku, to jest do uniknięcia.

x [i]

$x[i]$

i

$i$

x [\pm i]

$x[\pm i]$

x [i]

$x[i]$

(N - i)

$(N-i)$

x [N - i]

$x[N-i]$

i

$i$

i

$i$

— Dilip Sarwate

1

@Leo Mnożenie elementów. tablica n-po-1 x tablica n-po-1 = tablica n-po-1 W odpowiedzi nazwałem to „próbka po próbce”.

— endolith,

17

Korelacja jest sposobem na wyrażenie podobieństwa dwóch szeregów czasowych (w twoim przypadku próbek audio) w jednym numerze. Jest to adaptacja kowariancji, która jest realizowana w następujący sposób:

period = 1/sampleFrequency;
covariance=0;

for (iSample = 0; iSample<nSamples; iSample++)
    covariance += (timeSeries_1(iSample)*timeSeries_2(iSample))/period;
    //Dividing by `period` might not even be necessary

Korelacja jest znormalizowaną wersją kowariancji, która jest kowariancją podzieloną przez iloczyn standardowych odchyleń obu szeregów czasowych. Korelacja da 0, gdy nie będzie korelacji (całkowicie nie podobne) i 1 dla korelacji całkowitej (całkowicie podobne).

Można sobie wyobrazić, że dwie próbki dźwięku mogą być podobne, ale nie są zsynchronizowane. Właśnie tutaj pojawia się korelacja krzyżowa . Obliczasz korelację między szeregami czasowymi, w których jedna z nich została przesunięta o jedną próbkę:

for (iShift=0; iShift<nSamples; iShift++)
    xcorr(iShift) = corr(timeSeries_1, timeSeries_2_shifted_one_sample);

Następnie wyszukaj maksymalną wartość w corrserii i gotowe. (lub przestań, jeśli znajdziesz wystarczającą korelację) Oczywiście jest coś więcej. Musisz zaimplementować standardowe odchylenie i wykonać zarządzanie pamięcią i wprowadzić zmiany czasu. Jeśli wszystkie twoje próbki dźwięku są równej długości, możesz obejść się bez normalizacji kowariancji i przystąp do obliczania kowariancji krzyżowej.

Fajny związek z twoim wcześniejszym pytaniem : analiza Fouriera jest tylko adaptacją kowariancji krzyżowej. Zamiast przesuwać jedną serię czasową i obliczać kowariancje z drugim sygnałem, oblicza się kowariancje między jednym sygnałem a liczbą fal (ko) sinusoidalnych o różnych częstotliwościach. Wszystko opiera się na tej samej zasadzie.

— Społeczność
źródło

1

Wspomniałeś, że 0 to brak korelacji, a 1 to korelacje całkowite. Chcę tylko zauważyć, że -1 jest całkowicie ujemnie skorelowane. Podobnie jak w -1, oznacza to, że próbka 1 jest przeciwieństwem próbki 2. Jeśli pomyślisz o tym na wykresie X, Y, jest to linia o dodatnim nachyleniu w stosunku do linii o ujemnym nachyleniu. A gdy zbliżasz się do zera, linia staje się „grubsza”.

— Kellenjb

@kellenjb, tak, ale pewnie bym to powiedział, wielkość korelacji, czym prawdopodobnie jesteś zainteresowany. 1 lub -1 oznaczają, że sygnały bezpośrednio na siebie wpływają.

— Kortuk

14

W przetwarzaniu sygnału korelacja krzyżowa (xcorr w MATLAB) jest operacją splotu z odwróconą jedną z dwóch sekwencji. Ponieważ odwrócenie czasu odpowiada złożonej koniugacji w dziedzinie częstotliwości, można użyć DFT do obliczenia korelacji krzyżowej w następujący sposób:

R_xy = ifft(fft(x,N) * conj(fft(y,N)))

gdzie N = rozmiar (x) + rozmiar (y) - 1 (najlepiej w zaokrągleniu do potęgi 2) jest długością DFT.

Mnożenie DFT jest równoważne zwojowi kołowemu w czasie. Zero wypełniania obu wektorów do długości N zapobiega nakładaniu się kołowo przesuniętych komponentów y na x, co czyni wynik identycznym z liniowym splotem xi odwróconym czasem y.

Opóźnienie 1 jest przesunięciem kołowym w prawo o y, natomiast opóźnienie -1 jest przesunięciem kołowym w lewo. Korelacja krzyżowa jest po prostu sekwencją iloczynu punktowego dla wszystkich opóźnień. W oparciu o standardowe porządkowanie fft będą one w tablicy, do której można uzyskać dostęp w następujący sposób. Wskaźniki od 0 do rozmiaru (x) -1 są dodatnimi opóźnieniami. Wskaźniki Rozmiar N (y) +1 do N-1 są opóźnieniami ujemnymi w odwrotnej kolejności. (W Pythonie do ujemnych opóźnień można wygodnie uzyskać dostęp za pomocą ujemnych wskaźników, takich jak R_xy [-1].)

Możesz myśleć o zerach x i y jako wektorach N-wymiarowych. Iloczyn iloczynu xiy dla danego opóźnienia wynosi |x|*|y|*cos(theta). Normy xiy są stałe dla przesunięć kołowych, więc ich podzielenie pozostawia tylko zmienny cosinus kąta theta. Jeśli xiy (dla danego opóźnienia) są ortogonalne w przestrzeni N, korelacja wynosi 0 (tj. Theta = 90 stopni). Jeśli są współliniowe, wartość wynosi 1 (korelacja dodatnia) lub -1 (korelacja ujemna, tj. Theta = 180 stopni). Prowadzi to do korelacji krzyżowej znormalizowanej do jedności:

R_xy = ifft(fft(x,N) * conj(fft(y,N))) / (norm(x) * norm(y))

Można to uczynić bezstronnym poprzez ponowne obliczenie norm tylko dla nakładających się części, ale wtedy równie dobrze możesz wykonać całe obliczenia w dziedzinie czasu. Zobaczysz także różne wersje normalizacji. Zamiast znormalizować do jedności, czasami korelacja krzyżowa jest znormalizowana przez M (tendencyjne), gdzie M = max (rozmiar (x), rozmiar (y)) lub M- | m | (bezstronna ocena m-tego opóźnienia).

Aby uzyskać maksymalne znaczenie statystyczne, średnią (błąd systematyczny DC) należy usunąć przed obliczeniem korelacji. Nazywa się to kowariancją krzyżową (xcov w MATLAB):

x2 = x - mean(x)
y2 = y - mean(y)
phi_xy = ifft(fft(x2,N) * conj(fft(y2,N))) / (norm(x2) * norm(y2))

— Eryk Sun
źródło

Czy to oznacza, że ostateczny rozmiar tablicy powinien wynosić 2*size (a) + size(b) - 1lub 2*size (b) + size (a) - 1? Ale w obu przypadkach dwa wyściełane tablice mają różne rozmiary. Jaka jest konsekwencja wypełniania zbyt dużą liczbą zer?

@RobertK Tablica korelacji krzyżowej musi mieć długość co najmniej sumę długości aib (minus jeden), jak mówi eryksun w swojej odpowiedzi. Dla uproszczenia często przyjmuje się, że długość jest dwa razy większa niż dłuższy wektor (czasami zaokrąglana w górę do następnej większej potęgi w celu użycia wydajnego FFT). Wybór pomaga, gdy klient z opóźnieniem decyduje, że chce także autokorelacji dłuższego wektora. Jedną konsekwencją wypełnienia zbyt dużą liczbą zer jest dodatkowe obliczenie, ale można to poprawić dzięki bardziej wydajnym implementacjom FFT.

2

$2$

— Dilip Sarwate

@RobertKJ: Jesteś przesuwne bwraz az jednym wyjściem na zmianę, minimum zachodzenia jednej próbce. Daje to size(a)opóźnienia dodatnie i size(b) - 1ujemne. Wykorzystując odwrotną transformatę iloczynu N-punktowych DFT, indeksy 0przechodzące size(a)-1są opóźnieniami dodatnimi, a wskaźniki N-size(b)+1przechodzące N-1są opóźnieniami ujemnymi w odwrotnej kolejności.

— Eryk Sun

3

jeśli używasz Matlaba, wypróbuj funkcję korelacji krzyżowej:

c= xcorr(x,y)

Oto dokumentacja Matlab:

xcorrszacuje sekwencję korelacji krzyżowej losowego procesu. Autokorelacja jest traktowana jako szczególny przypadek.

...

c = xcorr(x,y)zwraca sekwencję korelacji krzyżowej w wektorze długości 2 * N-1, gdzie xi ysą Nwektorami długości ( N > 1). Jeśli xi ynie są tej samej długości, krótszy wektor jest dopełniany zerami do długości dłuższego wektora.

korelacja http://www.mathworks.com/help/toolbox/signal/ref/eqn1263487323.gif

— smashtastic
źródło

Link wydaje się być zepsuty.

— Danijel

2

Szybki i prosty sposób porównywania plików audio. Weź plik audio, zrób kopię, od razu wklej je obok siebie, w 2 kanałach stereo, odwróć fazę na jednej ze ścieżek stereo, wyrównaj oba pliki na początku w trybie powiększenia, upewnij się, że oba pliki mają na początku tę samą amplitudę, a następnie odtwarzaj, jeśli panuje całkowita cisza, to oba pliki są identyczne, jeśli jest jakaś różnica, usłyszysz to całkiem wyraźnie !.

— użytkownik31971
źródło

1

Jak większość tutaj napisała, powinieneś użyć korelacji.

Wystarczy wziąć pod uwagę 2 czynniki:

Jeśli wolumin jest skalowany inaczej, należy znormalizować korelację.
W przypadku skalowania czasu możesz użyć dynamicznego dopasowania czasu.

— David
źródło

1

Dla sygnałów nieokresowych (wielkość (y) -1) należy odjąć od indeksu R_xy, aby uzyskać rzeczywiste opóźnienie.

N = rozmiar (x) + rozmiar (y) - 1;

opóźnienia = [0, N] - (rozmiar (y) - 1);

— Patrick
źródło

0

Najłatwiejszym sposobem znalezienia różnicy, IMO, jest odjęcie dwóch sygnałów audio w dziedzinie czasu. Jeśli są równe, wynik w każdym punkcie czasu będzie wynosił zero. Jeśli nie są równe, różnica między nimi pozostanie po odjęciu i możesz słuchać bezpośrednio. Szybką miarą ich podobieństwa byłaby wartość RMS tej różnicy. Często odbywa się to podczas miksowania i masteringu dźwięku, aby usłyszeć różnicę na przykład pliku MP3 w stosunku do WAV. (Odwrócenie fazy jednego sygnału i dodanie go jest tym samym, co odejmowanie. Jest to metoda stosowana, gdy odbywa się to w oprogramowaniu DAW). Muszą być idealnie wyrównane czasowo, aby to zadziałało. Jeśli nie są, możesz opracować algorytm do ich wyrównywania, taki jak wykrywanie dziesięciu najlepszych pików, obliczanie średniego przesunięcia pików i przesuwanie jednego sygnału.

Przekształcenie w dziedzinę częstotliwości i porównanie widm mocy sygnałów, tak jak proponujesz, ignoruje niektóre informacje w dziedzinie czasu. Na przykład dźwięk odtwarzany do tyłu miałby to samo widmo, gdy odtwarzany był do przodu. Zatem dwa bardzo różne sygnały audio mogą mieć dokładnie to samo widmo.

— Martin Vandepas
źródło