Оценка временной задержки сигналов осциллографа с использованием взаимной корреляции

Я записал 2 сигнала от оскопа. Они выглядят так:

Я хочу измерить задержку между ними в Matlab. Каждый сигнал имеет 2000 отсчетов с частотой дискретизации 2001000,5.

Данные находятся в файле CSV. Это то, что я до сих пор.
Я удалил данные времени из файла CSV, чтобы в файле CSV были только уровни напряжения.

x1 = csvread('C://scope1.csv');
x2 = csvread('C://scope2.csv');  
cc = xcorr(x1,x2);
plot(cc);  

Это дает такой результат:

Из того, что я прочитал, мне нужно взять взаимную корреляцию этих сигналов, и это должно дать мне пик, связанный с временной задержкой. Однако, когда я беру взаимную корреляцию этих сигналов, я получаю пик в 2000 году, который, я знаю, не верен. Что я должен сделать с этими сигналами, прежде чем я буду взаимно коррелировать их? Просто ищу направление.

РЕДАКТИРОВАТЬ: после удаления смещения постоянного тока это результат, который я теперь получаю:

Есть ли способ убрать это, чтобы получить более определенную задержку?

РЕДАКТИРОВАТЬ 2: Вот файлы:
http://dl.dropbox.com/u/10147354/scope1col.csv
http://dl.dropbox.com/u/10147354/scope2col.csv

@NickS

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

В этом случае я считаю, что может сработать оценка времени значительной энергии сигналов. Конечно, «значительный» может или не может быть несколько субъективным, но я считаю, что, взглянув на ваши сигналы со статистической точки зрения, мы сможем количественно оценить «значительный» и перейдем оттуда.

Для этого я сделал следующее:

ШАГ 1: Рассчитайте огибающие сигнала:

Этот шаг прост, так как вычисляется абсолютное значение выходного преобразования Гильберта для каждого из ваших сигналов. Существуют и другие методы вычисления конвертов, но это довольно просто. Этот метод по существу вычисляет аналитическую форму вашего сигнала, другими словами, векторное представление. Когда вы берете абсолютное значение, вы разрушаете фазу и только после энергии.

Кроме того, поскольку мы проводим оценку времени задержки энергии ваших сигналов, такой подход оправдан.

ШАГ 2: Удаление шума с сохраняющих края нелинейных медиальных фильтров:

Это важный шаг. Цель здесь состоит в том, чтобы сгладить ваши энергетические оболочки, но без разрушения или сглаживания ваших краев и быстрого времени нарастания. На самом деле этому посвящено целое поле, но для наших целей мы можем просто использовать простой в реализации нелинейный медиальный фильтр . Медианная фильтрация. Это мощный метод, потому что, в отличие от средней фильтрации, медиальная фильтрация не обнулит ваши края, но в то же время «сгладит» ваш сигнал без значительного ухудшения важных краев, так как ни в коем случае не выполняется какая-либо арифметика с вашим сигналом (при условии, что длина окна нечетна). Для нашего случая здесь я выбрал медиальный фильтр с размером окна 25 сэмплов:

ШАГ 3: Удалить время: Построить функции оценки плотности ядра по Гауссу:

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

Что интуитивно понятно из рисунка выше? Несмотря на то, что шумовая энергия мала, она имеет преимущество в том, что она более «популярна». Напротив, в то время как огибающая сигнала, которая имеет энергию, является более энергичной, чем шум, она фрагментирована через пороги. Что если бы мы рассматривали «популярность» как меру энергии? Это то, что мы будем делать с (моей грубой) реализацией функции плотности ядра (KDE) с гауссовым ядром.

Для этого берется каждая выборка и строится гауссова функция, используя ее значение в качестве среднего значения и предварительно установленную полосу пропускания (дисперсию), выбранную априори. Установка дисперсии вашего гауссиана является важным параметром, но вы можете установить его на основе статистики шума, основанной на вашем приложении и типичных сигналах. (У меня есть только ваши 2 файла, чтобы уйти). Если затем мы построим оценку KDE, мы получим следующий график:

Вы можете думать о KDE как о непрерывной форме гистограммы, а о дисперсии - о ширине бина. Однако у него есть преимущество, заключающееся в том, что он гарантирует гладкий формат PDF, так что мы можем затем выполнить первое и второе производные исчисления. Теперь, когда у нас есть гауссовские KDE, мы можем видеть, где образцы шума достигают пика популярности. Помните, что ось X здесь представляет проекции наших данных на амплитудное пространство. Таким образом, мы можем видеть, в каких порогах шум является наиболее «энергичным», а те говорят нам, каких порогов следует избегать.

На втором графике берется первая производная гауссовых KDE, и мы выбираем абсциссу первого образца после первой производной после пика смеси гауссианов, чтобы достичь определенного значения, близкого к нулю. (Или сначала пересечение нуля). Мы можем использовать этот метод и быть «безопасными», потому что наш KDE был построен из гладких гауссианов с разумной шириной полосы, и была взята первая производная от этой гладкой и бесшумной функции. (Как правило, первые производные могут быть проблематичными в любом случае, кроме сигналов с высоким SNR, поскольку они увеличивают шум).

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

Таким образом, это дает выборок.$\delta{t} = 241$

Я надеюсь, что это помогло.

Есть несколько проблем, связанных с автокорреляцией

1. Огромное смещение постоянного тока (уже исправлено)
2. Временное окно: Matcorb xcorr () имеет раздражающее соглашение по существу «обнулять» сигнал на обоих концах, когда вы сдвигаете временную задержку. Т.е. окно данных является функцией временной задержки. Это создаст треугольную форму для любого стационарного сигнала (включая синусоиды). Лучше выбрать окно корреляции таким образом, чтобы размер окна плюс максимальный временной интервал вписывались в ваше общее окно данных, или нормализовать взаимную корреляцию по количеству не дополняемых выборок.
3. Эти два сигнала не выглядят особенно коррелированными для меня. Форма в некоторой степени похожа, но конкретное расстояние между пиками и провалами совершенно иное, поэтому я сомневаюсь, что даже правильная автокорреляция даст много понимания.

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

Как указали пикенеты, в этом случае пик в середине выходного сигнала указывает на 0 лагов. Смещение пика от средней точки - это ваше временное отставание.

РЕДАКТИРОВАТЬ: Меня беспокоит, что корреляция почти идеальный треугольник. Это указывает на то, что взаимная корреляция не приводит к нормализации мощности. Это дает несправедливый уклон меньшим лагам по сравнению с большими лагами. Я бы изменил ваш вызов xcorr на «cc = xcorr (x1, x2,« беспристрастный »);».

Имейте в виду, что это не идеальное решение, потому что результаты с большими задержками теперь более нестабильны, чем результаты с низкими задержками, потому что они основаны на меньшем количестве данных. Большой пик на конечностях может быть фальшивым по той же причине, по которой вы можете получить 100% головы и без хвостов всего за несколько бросков монет, хотя на многих бросках это крайне маловероятно.

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

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

Однако есть случаи, когда взаимная корреляция полезна. Скажите, что ваш второй сигнал был действительно оригинальной версией со сдвигом во времени, даже с добавлением некоторого дополнительного шума:

a = csvread('scope1col.csv');
a = a - mean(a);               % to remove DC offset
b = a(200:end) + sqrt(0.05)*randn(1801,1);
figure; subplot(211); plot(a); subplot(212); plot(b)

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

[c,lags] = xcorr(a,b);
igure; plot(lags,c); grid on; xlabel('Lag'); ylabel('Cross-correlation');

который показывает пик при правильной задержке 200 образцов. Корреляция может быть полезным инструментом для определения временной задержки при применении к наборам данных, которые содержат правильный тип сходства.

Основываясь на предложении Мухаммеда, я попытался написать сценарий Matlab. Однако я не могу сделать вывод, строит ли он распределение Гаусса на основе дисперсий, а затем принимает оценку KDE или выполняет оценку KDE с предположением Гаусса.

Также трудно понять, как он переводит время смещения KDE во временную область. Вот моя попытка этого. Любой пользователь, который заинтересован в использовании скрипта, может бесплатно обновить улучшенную версию, если это возможно.

%% Initialising data

Ws1 = data1;
Ws2 = data2;
mWs1 = nanmean(Ws1);
mWs2 = nanmean(Ws2);
sdWs1 = nanstd(Ws1);
sdWs2 = nanstd(Ws2);

%% Computing the signal envelopes
Ws1d = Ws1 - mWs1;
Ws2d = Ws2 - mWs2;
h1 = abs(hilbert(Ws1d));
h2 = abs(hilbert(Ws2d));
figure();
subplot(211)
plot([Ws1d, h1])
subplot(212)
plot([Ws2d, h2])

%% Denoise the signal with edge preserving nonlinear medial filtering
w = 25;
mf1 = medfilt1(h1, w);
mf2 = medfilt1(h2, w);
figure();
subplot(211)
plot(mf1)
subplot(212)
plot(mf2)

<%% Remove time: construct the gaussian kernel density estimation functions>
% Using the kde from Matlab central directly on the filtered data
data1 = mf1;
[bw1, den1, xmesh1, cdf1] = kde(data1, 2^14);
der1 = diff(den1);
data2 = mf2;
[bw2, den2, xmesh2, cdf2] = kde(data2, 2^14);
der2 = diff(den2);
figure();
plot([der1, der2]);
legend('Sig1', 'Sig2')

% the other method as explained in Muhammad's post
for i = 1:length(mf1)
gf1(:,i) = mf1(i) + sdWs1*randn(1000,1);
gf2(:,i) = mf2(i) + sdWs2*randn(1000,1);
end
[bwM1, denM1, xmeshM1, cdfM1] = kde(gf1(:,1), 2.^11);
dd1 = diff(denM1);
[bwM2, denM2, xmeshM2, cdfM2] = kde(gf2(:,1), 2.^11);
dd2 = diff(denM2);
figure();
plot([dd1, dd2]);
legend('Sig1', 'Sig2')