我已经从显微镜记录了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.csvhttp://dl.dropbox.com/u/10147354 /scope2col.csv

评论

确切地说,您是如何进行互相关的?在回答您的直接问题时,互相关之前不需要对信号做任何事情,尽管在某些情况下,首先进行滤波有助于消除可能使结果失真的噪声。

请发布您已使用的代码,更重要的是发布互相关信号的曲线图。一些工具/库将(lag = 0)得分放在图表的中间;我不记得Matlab是否这样做。

@pichenettes:更新后的帖子

@JimClay:更新后的帖子

@NickS。如果信号完全对齐,则在cc图的中间会出现一个峰值。因此,在2000年达到峰值就意味着没有延迟。现在让我们说您有10个采样的延迟,这意味着signal2与signal1相比有10个采样。这只会使您的抄送高峰从2000年移至2010年(或1990年)。因此,您的时间延迟对应于您的实际峰值位置,即MINUS2000。

#1 楼

@NickS

由于无法确定图中的第二个信号实际上是第一个信号的唯一延迟版本,因此,除了经典的互相关之外,还必须尝试其他方法。这是因为如果您的信号是彼此的延迟版本,则互相关(CC)仅仅是最大似然估计器。在这种情况下,它们显然也不是,它们的非平稳性也不用说了。

在这种情况下,我认为可能有用的是对信号的大量能量进行时间估计。诚然,“重要”在某种程度上可以是主观的,但我相信,从统计学的角度看待您的信号,我们将能够量化“重要”并从那里去。

为此,我执行了以下操作:

步骤1:计算信号包络:

此步骤很简单,因为绝对值计算每个信号的希尔伯特变换的输出。还有其他计算信封的方法,但这很简单。此方法本质上计算信号的解析形式,即相量表示。当您取绝对值时,您将破坏阶段并且只有在消耗能量之后。

此外,由于我们要对您的信号能量进行时延估算,因此这种方法是必要的。



步骤2:使用保留边缘的非线性中间滤波器对噪声进行除噪:

这是重要的一步。此处的目的是使能量包络平滑,但不破坏或平滑边缘和快速上升时间。实际上有一个专门的领域,但是出于我们这里的目的,我们可以简单地使用易于实现的非线性Medial滤波器。 (中值过滤)。这是一种强大的技术,因为与均值滤波不同,中间滤波不会使边缘变零,但同时会“平滑”信号,而不会明显降低重要的边缘,因为在任何时候都不会对信号执行任何算术运算(假设窗口长度为奇数)。对于我们这里的情况,我选择了一个窗口大小为25个样本的中间滤波器:



第3步:删除时间:构造高斯核密度估计函数:

如果不按常规方式看上面的图,会发生什么?从数学上讲,这意味着,如果将每个去噪信号样本投影到y振幅轴上,将会得到什么?通过这样做,我们将设法节省时间,并且能够仅研究信号统计。

直觉上看是什么?尽管噪声能量较低,但它具有“受欢迎”的优点。相反,尽管具有能量的信号包络比噪声更有能量,但它在阈值之间分散。如果我们将“人气”视为一种能量度量,该怎么办?这就是我们(粗略地)使用高斯内核实现内核密度函数(KDE)的方式。

为此,将采集每个样本,并使用其值作为平均值构造高斯函数,并选择先验带宽作为预设带宽(方差)。设置高斯方差是一个重要参数,但是您可以基于基于应用程序和典型信号的噪声统计信息来设置它。 (我只有两个文件要关闭)。如果然后构造KDE估计,则会得到以下图:



您可以将KDE看作直方图的连续形式,可以这么说。方差作为您的bin宽度。但是,这样做的好处是可以保证平滑的PDF,然后我们可以对它们执行一阶和二阶微积分。现在我们有了高斯KDE,我们可以看到噪声采样在哪里达到峰值。请记住,这里的x轴表示我们的数据在幅度空间上的投影。因此,我们可以看到噪声中最“充满活力”的阈值,而那些告诉我们应避免使用哪些阈值。

在第二个图中,获取了高斯KDE的一阶导数,我们选择了高斯混合峰之后的一阶导数之后的第一个样本的横坐标,以获得一定值接近零。 (或第一个零交叉)。我们可以使用这种方法并且是“安全的”,因为我们的KDE是由具有合理带宽的平滑高斯构成的,并且采用了该平滑且无噪声的函数的一阶导数。 (通常,一阶导数会在高SNR信号以外的任何情况下产生问题,因为它们会放大噪声。)

黑线表明了在什么阈值下将图像“分段”是明智的,这样我们避免了整个本底噪声。然后,如果我们将其应用于原始信号,则会得到以下图,其中黑线表示信号能量的开始:



因此得出一个$ \ delta {t} = 241 $个样本。

希望对您有所帮助。

评论


$ \ begingroup $
哇。非常感谢。这些都是我将要开始研究的新技术。有什么办法可以看看您使用的matlab代码吗?
$ \ endgroup $
–尼克·西纳斯(Nick Sinas)
2012年4月14日在3:12



$ \ begingroup $
因此,我在Matlab中完成了步骤1和2,结果与您的相符,但是步骤3出现了问题。您使用了什么功能?
$ \ endgroup $
–尼克·西纳斯(Nick Sinas)
2012年4月14日15:40在

$ \ begingroup $
@NickS。询问,然后您会收到一封电子邮件给我,我可以将您使用的代码发送给您。
$ \ endgroup $
–太空
2012年4月14日在17:48

$ \ begingroup $
@Mohammed能否请您发布代码以估计时间延迟。我已就此事向您发送了电子邮件,请帮助
$ \ endgroup $
–user4856
13年6月20日在9:15

#2 楼

使用自相关操作时会遇到一些问题。


巨大的DC偏移(已修复)
时间窗口:Matlab的xcorr()具有令人讨厌的约定,即基本上“零填充”了滑动时滞时,两端都发出信号。即数据窗口是时间滞后的函数。这将为任何平稳信号(包括正弦波)创建一个三角形。更好的选择是选择一个相关窗口,以使窗口大小加上最大时滞适合您的总数据窗口,或者通过未填充样本的数量对互相关进行归一化。
这两个信号看起来并不特别相关对我来说。形状有些相似,但峰和谷的特定间距却大不相同,因此我怀疑即使适当的自相关也会在此处产生很多见识。

一种更简单的方法是使用阈值检测器查找起点,只需使用这些点之间的差作为延迟即可。

#3 楼

如小插图所示,在这种情况下,输出中间的峰值表示0滞后。峰值与中点的偏移量就是您的时间滞后。

编辑:与我相关的是,这种相关性几乎是一个完美的三角形。这向我表明,互相关没有进行功率归一化。相对于较大的滞后,这给较小的滞后带来了不公平的偏差。我会将您的xcorr调用修改为“ cc = xcorr(x1,x2,'unbiased');“。

,请注意,这不是一个完美的解决方案,因为滞后结果现在更多了比低滞后结果不稳定,因为它们基于更少的数据。四肢的大峰可能是虚假的,原因相同,只需几次抛硬币就可以实现100%正反,而多次抛则极不可能发生。

评论


$ \ begingroup $
信号没有延迟吗?
$ \ endgroup $
–尼克·西纳斯(Nick Sinas)
2012年4月6日在19:12

$ \ begingroup $
我不确定-高峰在哪里?我可以看到它在中间,但是不清楚它实际上在中间。此外,还有一个电源标准化问题,我将在编辑答案时解决。
$ \ endgroup $
– Jim Clay
2012年4月6日19:15

$ \ begingroup $
'unbiased'参数无疑使它看起来更好。我期望的更多。我将继续调查。谢谢。
$ \ endgroup $
–尼克·西纳斯(Nick Sinas)
2012年4月6日19:50

$ \ begingroup $
@JimClay也许Nick S正在关联其信号的包络而不是实际信号(Nick是真的吗?)。我想象这将产生(大致)这个三角形。
$ \ endgroup $
–太空
2012年4月6日19:56

$ \ begingroup $
@NickS。穆罕默德(Mohammad)的评论让我看起来并意识到,您有一个巨大的DC偏移量,使您的结果混乱。从两个信号中减去均值,然后对它们运行xcorr。我会先尝试不带“无偏”选项的情况。
$ \ endgroup $
– Jim Clay
2012年4月6日在20:01

#4 楼

正如其他人指出的那样,看来您是根据对问题的最后编辑而意识到的,似乎互相关不会为您提供一个很好的估计显示数据集的时间延迟的方法。相关性通过在一个时间间隔内使一个时间序列在另一个时间序列上滑动并在每个时间间隔计算两个序列之间的内积来度量两个时间序列之间的形状相似性。当两个序列在​​质量上相似或彼此“关联”时,结果将具有很大的幅度。这类似于当两个向量指向相同方向时两个向量的内积最大。

您显示的数据存在的问题是(至少对于我们可以看到的摘录)的形状似乎没有太多相似之处。可以无延迟地应用一个信号使其看起来与另一个信号一样,这正是通过计算它们的互相关来进行的操作。

在某些情况下,相关性是有用的。假设您的第二个信号确实是原始信号的时移版本,即使添加了一些额外的噪音:

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个样本的正确滞后出现的峰值。当应用于包含正确类型的相似性的数据集时,相关性可能是确定时间延迟的有用工具。

评论


$ \ begingroup $
关于我还能做什么的任何想法?可能是互相关之外的另一种技术,或者是某种类型的滤波器?谢谢。
$ \ endgroup $
–尼克·西纳斯(Nick Sinas)
2012年4月8日在6:16

$ \ begingroup $
@NickS。我也查看了它们,它们并不是彼此的延迟副本。话虽这么说,您是否想估算能量的延迟?我认为在这种情况下,信号的VS延迟会更有意义。如果您告诉我们有关您正在运行的基础渠道/实验的更多信息,我们可以为您提供有关可能路径的更多信息。
$ \ endgroup $
–太空
2012年4月9日在16:02

$ \ begingroup $
@Mohammad谢谢。下方的通道是钢。关于如何估计能量延迟的任何想法?
$ \ endgroup $
–尼克·西纳斯(Nick Sinas)
2012年4月9日在23:42

$ \ begingroup $
@Mohammad您是否认为信号的失真可能是某种混响类型,可以通过滤波消除?
$ \ endgroup $
–尼克·西纳斯(Nick Sinas)
2012年4月9日在23:44

$ \ begingroup $
@NickS。可能会有一些混响清理技巧(我不知道这些技巧是如何实现的),但是我已经拼凑了一些简单的方法,如果您想看一看,这将是一种能量估算器。
$ \ endgroup $
–太空
2012年4月12日在16:25

#5 楼

根据穆罕默德的建议,我尝试编写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')