Найти задержку сигнала или фазы из взаимной корреляции

#python #numpy #scipy #signal-processing #cross-correlation

#python #numpy #scipy #обработка сигналов #взаимная корреляция

Вопрос:

Я использую python, но это общий вопрос (больше связанный с алгоритмами и т. Д.), И поэтому я пропускаю некоторые шаги, чтобы понять суть вопроса:

Я генерирую синусоидальный сигнал следующим образом:

 import math as m
signal = [m.sin(2*m.pi*1*(t/n-d)) for t in range(n)]
  

Итак, синусоидальный сигнал, нормализованный таким образом, что частота равна 1, а время идет от 0 до 1 секунды (так что в основном это простой один цикл синусоидальной волны). Существует также термин задержки d, который задерживает сигнал (вызывает сдвиг фазы). N — это только количество выборок

Я также создаю другой сигнал с другой задержкой. Допустим, я использую задержку 0 для первого сигнала и задержку x для второго сигнала (я сокращаю предыдущий для ясности):

 signal1 = signal(delay=0)
signal2 = signal(delay=x)
  

и затем я делаю корреляцию:

 from scipy import signal as sgn
corr11 = sgn.correlate(signal1, signal1, mode = 'full')
corr12 = sgn.correlate(signal1, signal2, mode = 'full')
  

Я также знаю, что задержка сигнала коррелирует с максимумом точки корреляции, поэтому я вынимаю две точки:

 import numpy as np

a1 = np.argmax(corr11)
a2 = np.argmax(corr12)
  

Итак, я обнаружил, что корреляция сигнала с самим собой имеет максимальный пик в середине корреляционного массива (или графика / функции). Но другой пик странный:

  • При задержке 0 и 1: a2 совпадает с a1
  • При задержке 0,5: расстояние a2 от a1 равно 0,5 от a1 (инвертированный сигнал)
  • При задержке 0,28328: a2 составляет 0,75 от a1
  • При задержке 0,1: a2 составляет 0,90888 от a1

Итак, вопрос в том, как задержка d соотносится с местоположением пика после корреляции сигналов?

Комментарии:

1. Что такое nop ? Кроме того, sgn здесь имеется в виду scipy.signal ?

2. Nop была функцией, которая ничего не делала (не выполняла никаких операций), а sgn — scipy.signal да. Я исправил эти ошибки

Ответ №1:

Похоже, что задержка примерно равна (a1 - a2) / n . Однако я думаю, что ответ несколько искажен тем фактом, что а) вы используете только синусоидальную волну с одним периодом и б) вы используете конечное число точек данных (очевидно). Чтобы получить более точный ответ для случая синусоидальной волны с одним периодом, вы, вероятно, захотите получить математическое определение корреляции и выполнить необходимую интеграцию с правильными ограничениями (но я не уверен, что SO — правильное место, чтобы обратиться за помощью с интеграцией).

Вот автономный скрипт, который отображает сигналы и корреляции, что, мы надеемся, даст больше интуиции. ПРИМЕЧАНИЕ: приведенное выше приближение, по-видимому, становится более точным, когда вы повторяете количество периодов синусоидальной волны. Например, при 100 периодах и 100000 точках данных приведенное выше приближение (измененное здесь как n_repeats * (a1 - a2) / n ), кажется, становится намного более точным.

Скрипт

 import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# Set parameters

# x = 0.5
x = 0.28328
# x = 0.25
# x = 0.1
# n = 100000
# n_repeats = 100
n = 1000
n_repeats = 1

# Get correlations
t = np.linspace(0, n_repeats, n)

sin_delay = lambda delay: np.sin(2.0 * np.pi * (t - delay))

signal1 = sin_delay(delay=0)
signal2 = sin_delay(delay=x)

corr11 = signal.correlate(signal1, signal1, mode = 'full')
corr12 = signal.correlate(signal1, signal2, mode = 'full')

a1 = np.argmax(corr11)
a2 = np.argmax(corr12)

# Print output
print(a1, a2, x, n_repeats * (a1 - a2) / n)

# Make plots
plt.figure()
plt.plot(signal1, "r")
plt.plot(signal2, "b")
plt.title("Signals, delay = {:.3f}".format(x))
plt.legend(["Original signal", "Delayed signal"], loc="upper right")
plt.grid(True)
plt.savefig("Signals")
plt.figure()
plt.plot(corr11, "r")
plt.plot(corr12, "b")
plt.title("Correlations, delay = {:.3f}".format(x))
plt.legend(["Auto-correlation", "Cross-correlation"], loc="upper right")
plt.grid(True)
plt.savefig("Correlations")
  

Вывод на консоль с n = 1000, n_repeats = 1

 999 749 0.28328 0.25
  

Вывод на консоль с n = 100000, n_repeats = 100

 99999 99716 0.28328 0.283
  

Вывод изображений с n = 1000, n_repeats = 1

введите описание изображения здесь
введите описание изображения здесь