Как реализовать взаимную корреляцию, чтобы доказать, что два аудиофайла похожи?

58

Я должен сделать взаимную корреляцию двух аудио файлов, чтобы доказать, что они похожи. Я взял БПФ двух звуковых файлов, и их значения спектра мощности находятся в отдельных массивах.

Как мне продолжить кросс-корреляцию и доказать, что они похожи? Есть ли лучший способ сделать это? Любые основные идеи будут полезны для меня, чтобы изучить и применить его.

Лорем Ипсум
источник
Приведена взаимная корреляция двух случайных сигнальных векторов. Как реализовать обратное, чтобы получить два вектора в MATLAB. Джон Мухехе

Ответы:

56

Кросс-корреляция и свертка тесно связаны. Короче говоря, чтобы сделать свертку с БПФ, вы

  1. обнуляйте входные сигналы (добавьте нули в конец, чтобы хотя бы половина волны была «пустой»)
  2. принять БПФ обоих сигналов
  3. умножить результаты вместе (поэлементное умножение)
  4. сделать обратное БПФ

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

Вам необходимо выполнить заполнение нулями, потому что метод 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)))

что проще с вашим аппаратным / программным обеспечением. Для автокорреляции (взаимной корреляции сигнала с самим собой) лучше сделать комплексное сопряжение, потому что тогда вам нужно только рассчитать БПФ один раз.

Если сигналы действительны, вы можете использовать реальные БПФ (RFFT / IRFFT) и сэкономить половину времени вычислений, рассчитав только половину спектра.

Также вы можете сэкономить время вычислений, дополнив его большим размером, для которого оптимизировано FFT (например, 5-гладким числом для FFTPACK, ~ 13-гладким числом для FFTW или степенью 2 для простой аппаратной реализации).

Вот пример в Python корреляции БПФ по сравнению с корреляцией грубой силы: https://stackoverflow.com/a/1768140/125507

Это даст вам функцию взаимной корреляции, которая является мерой сходства против смещения. Чтобы получить смещение, при котором волны «выровнены» друг с другом, в корреляционной функции будет пик:

пик в корреляционной функции

Значение x пика - это смещение, которое может быть отрицательным или положительным.

Я видел только это, чтобы найти смещение между двумя волнами. Вы можете получить более точную оценку смещения (лучше, чем разрешение ваших выборок), используя параболическую / квадратичную интерполяцию на пике.

Чтобы получить значение подобия между -1 и 1 (отрицательное значение, указывающее, что один из сигналов уменьшается по мере увеличения другого), вам необходимо масштабировать амплитуду в соответствии с длиной входов, длиной FFT, вашей конкретной реализацией FFT. масштабирование и т. д. Автокорреляция волны с самим собой даст вам значение максимально возможного соответствия.

Обратите внимание, что это будет работать только на волнах, имеющих одинаковую форму. Если они были сэмплированы на другом оборудовании или имело некоторый добавленный шум, но в остальном они по-прежнему имеют одинаковую форму, это сравнение будет работать, но если форма волны была изменена фильтрацией или фазовыми сдвигами, они могут звучать одинаково, но выиграли не коррелирует также.

эндолиты
источник
3
Заполнение нулями должно быть не менее N = size (a) + size (b) -1, предпочтительно округлено до степени 2. Чтобы получить значение от -1 до 1, разделите на norm (a) * norm (b) ), который дает косинус угла между двумя векторами в N-пространстве для данного лага (т.е. круговое смещение по модулю N). На крайних лагах не так много перекрывающихся выборок (только один на дальнем крайнем расстоянии), поэтому деление на norm (a) * norm (b) сместит эти корреляции в сторону 0 (т. Е. Покажет их относительную ортогональность в N-пространстве) ,
Eryk Sun
1
Я думаю, что в описании может быть ошибка. Разве не следует умножать БПФ вместе термин за термином, чтобы получить БПФ свертки сигналов, а не БПФ взаимной корреляции ? Насколько я понимаю, чтобы получить БПФ кросс-корреляции, необходимо использовать комплексное сопряжение одного из векторов БПФ в терминах умножения перед каждым взятием iFFT.
Дилип Сарватэ
@DilipSarwate: Да, вы правы. Вы также можете повернуть один сигнал во времени, который я добавил к ответу.
Эндолит
1
«Почему в аппаратном обеспечении трудно изменить время?» Во многих случаях данные хранятся в систолических массивах в ожидании того, что вычисления являются локальными , т.е. , хранящийся в ячейке, взаимодействует только со своими ближайшими соседями . Посылка в ячейку # и отправка в ячейку # , и делать это для всех возрастающих затрат, Электропроводка задержки (и , следовательно , уменьшает максимальную достижимую частоту), а также, потому что все провода должны пересекаться друг с другом, создает проблемы маршрутизации. Следует избегать , если это возможно, и в этом случае, это можно избежать.i x [ ± i ] x [ i ] ( N - i ) x [ N - i ] i ix[i]ix[±i]x[i](Ni)x[Ni]ii
Дилип Сарвейт
1
@ Leo поэлементное умножение. Массив n-by-1 x n-by-1 array = n-by-1 массив Я назвал это «выборка за выборкой» в ответе.
эндолит
17

Корреляция - это способ выразить сходство двух временных рядов (аудиосэмплов в вашем случае) в одном числе. Это адаптация ковариации, которая реализуется следующим образом:

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

Корреляция - это нормализованная версия ковариации, которая представляет собой ковариацию, деленную на произведение стандартных отклонений обоих временных рядов. Корреляция даст 0, когда нет корреляции (полностью не похожая) и 1 для полной корреляции (полностью похожая).

Вы можете представить, что два образца звука могут быть похожими, но не синхронизированы. Вот где появляется взаимная корреляция . Вы вычисляете корреляцию между временными рядами, в которых один из них сдвинут на одну выборку:

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

Затем найдите максимальное значение в corrсерии, и все готово. (или прекратите, если вы нашли достаточную корреляцию) Конечно, это немного больше. Вы должны реализовать стандартное отклонение, и вы должны сделать некоторое управление памятью и реализовать сдвиг времени. Если все ваши аудиосэмплы равны по длине, вы можете обойтись без нормализации ковариации и продолжить вычислять кросс-ковариацию.

Интересное отношение к вашему предыдущему вопросу : анализ Фурье - это всего лишь адаптация кросс-ковариации. Вместо того, чтобы сдвигать один временной ряд и вычислять ковариации с другим сигналом, вы вычисляете ковариации между одним сигналом и количеством (син) синусоидальных волн с разными частотами. Все это основано на том же принципе.

Сообщество
источник
1
Вы упомянули, что 0 - нет корреляции, а 1 - общая корреляция. Я просто хочу отметить, что -1 полностью отрицательно коррелируется. Например, -1 означает, что выборка 1 противоположна выборке 2. Если вы думаете об этом на графике X, Y, это линия с положительным наклоном по сравнению с линией с отрицательным наклоном. И когда вы приближаетесь к 0, линия становится «толще».
Келленджб
@kellenjb, да, но я бы, наверное, это сказал, величина корреляции, что вас, вероятно, интересует. 1 или -1 означают, что сигналы напрямую влияют друг на друга.
Кортук
14

При обработке сигнала взаимная корреляция (xcorr в MATLAB) является операцией свертки с обращением одной из двух последовательностей. Поскольку обращение времени соответствует комплексному сопряжению в частотной области, вы можете использовать ДПФ для вычисления взаимной корреляции следующим образом:

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

где N = размер (х) + размер (у) - 1 (предпочтительно округленный до степени 2) - длина ДПФ.

Умножение ДПФ эквивалентно круговой свертке во времени. Нулевое заполнение обоих векторов до длины N удерживает циклически сдвинутые компоненты y от перекрытия с x, что делает результат идентичным линейной свертке x и обращенному по времени y.

Отставание 1 - это правое круговое смещение y, а отставание от -1 - левое круговое смещение. Кросс-корреляция - это просто последовательность точечных произведений для всех лагов. Основываясь на стандартном порядке fft, они будут в массиве, к которому можно получить доступ следующим образом. Индексы от 0 до размера (x) -1 являются положительными лагами. Индексы N-размера (y) от +1 до N-1 являются отрицательными лагами в обратном порядке. (В Python к отрицательным лагам можно обращаться с помощью отрицательных индексов, таких как R_xy [-1].)

Вы можете думать о дополненных нулями x и y как о N-мерных векторах. Точечное произведение x и y для данного отставания равно |x|*|y|*cos(theta). Нормы x и y постоянны для круговых сдвигов, поэтому их разделение оставляет только изменяющийся косинус угла тета. Если x и y (для данного лага) ортогональны в N-пространстве, корреляция равна 0 (тета = 90 градусов). Если они коллинеарны, значение равно 1 (положительно коррелируется) или -1 (отрицательно коррелируется, то есть тета = 180 градусов). Это приводит к взаимной корреляции, нормированной к единице:

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

Это можно сделать беспристрастным, пересчитав нормы только для перекрывающихся частей, но тогда вы можете также выполнить все вычисления во временной области. Также вы увидите разные варианты нормализации. Вместо того, чтобы быть нормализованным к единице, иногда взаимная корреляция нормализуется с помощью M (смещения), где M = max (размер (x), размер (y)) или M- | m | (непредвзятая оценка m-го лага).

Для максимальной статистической значимости среднее значение (смещение постоянного тока) должно быть удалено перед вычислением корреляции. Это называется кросс-ковариацией (xcov в 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
источник
Означает ли это, что окончательный размер массива должен быть 2*size (a) + size(b) - 1или 2*size (b) + size (a) - 1? Но в любом случае два дополненных массива имеют разные размеры. Каковы последствия заполнения с большим количеством нулей?
@RobertK Массив кросс-корреляции должен иметь длину не менее суммы длин a и b (минус один), как говорит Эриксун в своем ответе. Для простоты длина часто принимается равной двойной длине более длинного вектора (иногда округляется до следующей большей степени , чтобы использовать эффективное БПФ). Выбор помогает, когда клиент с опозданием решает, что он также хочет автокорреляцию более длинного вектора. Одним из следствий заполнения нулями слишком большого количества нулей является дополнительное вычисление, но это может быть улучшено более эффективными реализациями FFT. 2
Дилип Сарватэ
@RobertKJ: Вы скольжение bвдоль a, с одним выходом в смену, с минимальным перекрытием одного образца. Это дает size(a)положительные лаги и size(b) - 1отрицательные лаги. Используя обратное преобразование произведения N-точечных ДПФ, 0сквозные индексы size(a)-1являются положительными лагами, а N-size(b)+1сквозные индексы N-1- отрицательными лагами в обратном порядке.
Eryk Sun
3

если вы используете Matlab, попробуйте функцию взаимной корреляции:

c= xcorr(x,y)

Вот документация Matlab:

xcorrоценивает кросс-корреляционную последовательность случайного процесса. Автокорреляция рассматривается как особый случай.

...

c = xcorr(x,y)возвращает последовательность взаимной корреляции в векторе длины 2 * N-1, где xи y- Nвекторы длины ( N > 1). Если xи yне имеют одинаковую длину, более короткий вектор дополняется нулями до длины более длинного вектора.

корреляция http://www.mathworks.com/help/toolbox/signal/ref/eqn1263487323.gif

smashtastic
источник
Ссылка, кажется, не работает.
Даниэль
2

Быстрый и простой способ сравнить аудио файлы. Возьмите аудиофайл, сделайте копию, в галку, вставьте их рядом, в 2 стереоканала, инвертируйте фазу на одной из стереотреков, выровняйте оба файла в начале в режиме масштабирования, убедитесь, что оба файла имеют одинаковую амплитуду в начале, затем воспроизводите, если есть полная тишина, тогда оба файла идентичны, если есть разница, вы услышите это довольно четко !.

user31971
источник
1

Как большинство здесь написали, вы должны использовать корреляцию.

Просто учтите 2 фактора:

  1. Если громкость масштабируется по-другому, вы должны нормализовать корреляцию.
  2. Если есть масштабирование времени, вы можете использовать Dynamic Time Warping.
Дэвид
источник
1

Для непериодических сигналов (размер (y) -1) необходимо вычесть из индекса R_xy, чтобы получить фактическое запаздывание.

N = размер (х) + размер (у) - 1;

лаги = [0, N] - (размер (у) - 1);

Патрик
источник
0

Самый простой способ найти разницу, IMO, - вычесть два аудиосигнала во временной области. Если они равны, результат в каждый момент времени будет равен нулю. Если они не равны, разница между ними останется после вычитания, и вы можете слушать ее напрямую. Быстрая оценка того, насколько они похожи, будет среднеквадратичным значением этой разницы. Это часто делается в микшировании и мастеринге аудио, например, чтобы услышать разницу между MP3 и WAV-файлами. (Инверсия фазы одного сигнала и добавление их - это то же самое, что вычитание. Этот метод используется, когда это делается в программном обеспечении DAW.) Они должны быть идеально выровнены по времени, чтобы это работало. Если это не так, вы можете разработать алгоритм их выравнивания, такой как обнаружение первых десяти пиков, вычисление среднего смещения пиков и смещение одного сигнала.

Преобразование в частотную область и сравнение спектров мощности сигналов, которые вы предлагаете, игнорирует некоторую информацию во временной области. Например, звук, воспроизводимый в обратном направлении, будет иметь тот же спектр при воспроизведении вперед. Таким образом, два очень разных аудиосигнала могут иметь одинаковый спектр.

Мартин Вандепас
источник