Мне нужно сравнить два сигнала с временным напряжением. Из-за особенностей источников этих сигналов одна из них может быть измененной по времени версией другой.
Как я могу найти, есть ли временной сдвиг? и если да, то сколько это.
Я делаю это в Python и хочу использовать библиотеки numpy/scipy.
scipy предоставляет функцию корреляции, которая будет хорошо работать для небольших входных данных, а также, если вы хотите некруговую корреляцию, означающую, что сигнал не будет распространяться. обратите внимание, что в mode='full' размер массива, возвращаемого signal.correlation, является суммой размеров сигналов минус один (т.е. len(a) + len(b) - 1), поэтому значение из argmax отключается на (размер сигнала -1 = 20) от того, что вы ожидаете.
from scipy import signal, fftpack
import numpy
a = numpy.array([0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0, 0, 0, 0])
b = numpy.array([0, 0, 0, 0, 0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0])
numpy.argmax(signal.correlate(a,b)) -> 16
numpy.argmax(signal.correlate(b,a)) -> 24
Два разных значения соответствуют тому, находится ли сдвиг в a или b.
Если вам нужна круговая корреляция и для большого размера сигнала, вы можете использовать теорему свертки/преобразования Фурье с оговоркой, что корреляция очень похожа, но не идентична свертке.
A = fftpack.fft(a)
B = fftpack.fft(b)
Ar = -A.conjugate()
Br = -B.conjugate()
numpy.argmax(numpy.abs(fftpack.ifft(Ar*B))) -> 4
numpy.argmax(numpy.abs(fftpack.ifft(A*Br))) -> 17
опять же, эти два значения соответствуют тому, интерпретируете ли вы сдвиг в a или сдвиг в b.
Отрицательное сопряжение связано со свертыванием одной из функций, но в корреляции это не происходит. Вы можете отменить переключение, изменив один из сигналов и затем взяв БПФ, или взяв БПФ сигнала, а затем взяв отрицательное сопряжение. то есть верно следующее: Ar = -A.conjugate() = fft(a[::-1])
Если один сдвинут во времени другим, вы увидите пик корреляции. Поскольку вычисление корреляции является дорогостоящим, лучше использовать БПФ. Итак, что-то вроде этого должно работать:
af = scipy.fft(a)
bf = scipy.fft(b)
c = scipy.ifft(af * scipy.conj(bf))
time_shift = argmax(abs(c))
Эта функция, вероятно, более эффективна для реальных сигналов. Он использует rfft и zero pads входы для мощности 2 достаточно большой, чтобы обеспечить линейную (то есть некруговую) корреляцию:
def rfft_xcorr(x, y):
M = len(x) + len(y) - 1
N = 2 ** int(np.ceil(np.log2(M)))
X = np.fft.rfft(x, N)
Y = np.fft.rfft(y, N)
cxy = np.fft.irfft(X * np.conj(Y))
cxy = np.hstack((cxy[:len(x)], cxy[N-len(y)+1:]))
return cxy
Возвращаемое значение - это длина M = len(x) + len(y) - 1 (взломанная вместе с hstack, чтобы удалить дополнительные нули от округления до степени 2). Неотрицательные лаги cxy[0], cxy[1], ..., cxy[len(x)-1], а отрицательные лага - cxy[-1], cxy[-2], ..., cxy[-len(y)+1].
Для согласования опорного сигнала, я бы вычислить rfft_xcorr(x, ref) и искать пика. Например:
def match(x, ref):
cxy = rfft_xcorr(x, ref)
index = np.argmax(cxy)
if index < len(x):
return index
else: # negative lag
return index - len(cxy)
In [1]: ref = np.array([1,2,3,4,5])
In [2]: x = np.hstack(([2,-3,9], 1.5 * ref, [0,3,8]))
In [3]: match(x, ref)
Out[3]: 3
In [4]: x = np.hstack((1.5 * ref, [0,3,8], [2,-3,-9]))
In [5]: match(x, ref)
Out[5]: 0
In [6]: x = np.hstack((1.5 * ref[1:], [0,3,8], [2,-3,-9,1]))
In [7]: match(x, ref)
Out[7]: -1
Это не надежный способ согласования сигналов, но это быстро и просто.
Это зависит от типа сигнала, который у вас есть (периодического?...), о том, имеют ли оба сигнала одинаковые амплитуды и какую точность вы ищете.
Функция корреляции, упомянутая highBandWidth, действительно может сработать для вас. Это достаточно просто, что вы должны попробовать.
Другим, более точным вариантом является тот, который я использую для высокоточного применения спектральной линии: вы моделируете свой "главный" сигнал с помощью сплайна и подбираете с ним сдвинутый по времени сигнал (хотя возможно масштабирование сигнала, если необходимо). Это дает очень точные сдвиги во времени. Одно из преимуществ этого подхода заключается в том, что вам не нужно изучать корреляционную функцию. Вы можете, например, легко создать сплайн с помощью interpolate.UnivariateSpline() (из SciPy). SciPy возвращает функцию, которая затем легко устанавливается с помощью optimize.leastsq().
Вот еще один вариант:
from scipy import signal, fftpack
def get_max_correlation(original, match):
z = signal.fftconvolve(original, match[::-1])
lags = np.arange(z.size) - (match.size - 1)
return ( lags[np.argmax(np.abs(z))] )