Найти сдвиг во времени между двумя одинаковыми формами волны
Я должен сравнить два сигнала времени-напряжения. Из-за особенностей источников этих сигналов одна из них может быть сдвинутой во времени версией другой.
Как я могу узнать, есть ли сдвиг во времени? и если да, то сколько это стоит.
Я делаю это на Python и хочу использовать библиотеки numpy/scipy.
5 ответов
scipy предоставляет функцию корреляции, которая отлично подойдет для небольшого входа, а также, если вы хотите некруговую корреляцию, означающую, что сигнал не будет распространяться. обратите внимание, что в mode='full', размер массива, возвращаемого signal.correlation, является суммой размеров входного сигнала - 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 и нулевые площадки для входов до степени 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
Это не надежный способ сопоставления сигналов, но он быстрый и простой.
Вот еще один вариант:
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))] )
Это зависит от типа вашего сигнала (периодического?...), от того, имеют ли оба сигнала одинаковую амплитуду, и от того, какую точность вы ищете.
Функция корреляции, упомянутая highBandWidth, действительно может работать для вас. Это достаточно просто, чтобы вы могли попробовать.
Другой, более точный вариант, который я использую для высокоточной подгонки спектральной линии: вы моделируете свой "мастер" сигнал с помощью сплайна и подгоняете смещенный по времени сигнал (при этом, возможно, масштабируя сигнал, если это необходимо). Это дает очень точные временные сдвиги. Одним из преимуществ этого подхода является то, что вам не нужно изучать корреляционную функцию. Например, вы можете легко создать сплайн с помощью interpolate.UnivariateSpline() (от SciPy). SciPy возвращает функцию, которая затем легко устанавливается optimize.leastsq().
Цитата
(Очень поздний ответ), чтобы найти временной сдвиг между двумя сигналами: используйте свойство временного сдвига FT, чтобы сдвиги могли быть короче, чем интервал между выборками, затем вычислите квадратичную разницу между сдвигом во времени сигнала и эталоном форма волны. Это может быть полезно, когда у вас есть n сдвинутых сигналов с множеством сдвигов, например, n приемников, равномерно распределенных для одной и той же входящей волны. Вы также можете скорректировать дисперсию, заменив статический временной сдвиг функцией частоты.
Код выглядит так:
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft, fftshift, fftfreq
from scipy import signal
# generating a test signal
dt = 0.01
t0 = 0.025
n = 512
freq = fftfreq(n, dt)
time = np.linspace(-n * dt / 2, n * dt / 2, n)
y = signal.gausspulse(time, fc=10, bw=0.3) + np.random.normal(0, 1, n) / 100
Y = fft(y)
# time-shift of 0.235; could be a dispersion curve, so y2 would be dispersive
Y2 = Y * np.exp(-1j * 2 * np.pi * freq * 0.235)
y2 = ifft(Y2).real
# scan possible time-shifts
error = []
timeshifts = np.arange(-100, 100) * dt / 2 # could be dispersion curves instead
for ts in timeshifts:
Y2_shifted = Y2 * np.exp(1j * 2 * np.pi * freq * ts)
y2_shifted = ifft(Y2_shifted).real
error.append(np.sum((y2_shifted - y) ** 2))
# show the results
ts_final = timeshifts[np.argmin(error)]
print(ts_final)
Y2_shifted = Y2 * np.exp(1j * 2 * np.pi * freq * ts_final)
y2_shifted = ifft(Y2_shifted).real
plt.subplot(221)
plt.plot(time, y, label="y")
plt.plot(time, y2, label="y2")
plt.xlabel("time")
plt.legend()
plt.subplot(223)
plt.plot(time, y, label="y")
plt.plot(time, y2_shifted, label="y_shifted")
plt.xlabel("time")
plt.legend()
plt.subplot(122)
plt.plot(timeshifts, error, label="error")
plt.xlabel("timeshifts")
plt.legend()
plt.show()