使用互相关估计示波器信号的时延

信息处理 matlab 互相关
2022-01-07 04:17:14

我已经记录了来自内窥镜的 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

4个回答

@尼克斯

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

在这种情况下,我相信可能有效的是对信号重要能量的时间估计。诚然,“显着”可以或不能有点主观,但我相信通过从统计角度查看您的信号,我们将能够量化“显着”并从那里开始。

为此,我做了以下工作:

第 1 步:计算信号包络:

这一步很简单,因为计算了每个信号的希尔伯特变换输出的绝对值。还有其他计算包络的方法,但这很简单。这种方法本质上是计算信号的分析形式,换句话说,就是相量表示。当您采用绝对值时,您正在破坏相位并且仅在能量之后。

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

在此处输入图像描述

第 2 步:使用保边非线性中值滤波器去噪:

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

在此处输入图像描述

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

如果你从侧面看上面的情节而不是正常的方式会发生什么?从数学上讲,这意味着,如果你将我们去噪信号的每个样本投影到 y 幅度轴上,你会得到什么?在这样做的过程中,可以说我们将设法消除时间,并且能够单独研究信号统计数据。

直观地从上图中弹出了什么?虽然噪音能量很低,但它的优势在于它更“流行”。相反,虽然具有能量的信号包络比噪声更有活力,但它在阈值上被分割。如果我们将“人气”视为能量的衡量标准会怎样?这就是我们将使用高斯内核实现内核密度函数(KDE)的(我的粗略)实现。

为了做到这一点,每个样本都被提取,并使用它的值作为平均值构造一个高斯函数,并先验地选择一个预设的带宽(方差)。设置高斯方差是一个重要参数,但您可以根据应用和典型信号的噪声统计数据进行设置。(我只有你的 2 个文件可以打开)。如果我们然后构造 KDE 估计,我们会得到以下图:

在此处输入图像描述

您可以将 KDE 视为直方图的连续形式,可以将方差视为您的 bin 宽度。然而,它具有保证平滑 PDF 的优点,然后我们可以对其执行一阶和二阶导数演算。现在我们有了高斯 KDE,我们可以看到噪声样本在哪里流行。请记住,这里的 x 轴代表我们的数据在幅度空间上的投影。因此,我们可以看到噪声在哪些阈值中最“活跃”,并告诉我们要避免哪些阈值。

在第二幅图中,取高斯KDE 的一阶导数,我们在高斯混合峰值后的一阶导数选取第一个样本的横坐标,以达到接近零的某个值。(或第一次过零)。我们可以使用这种方法并且是“安全的”,因为我们的 KDE 是由合理带宽的平滑高斯构造的,并且采用了这种平滑且无噪声的函数的一阶导数。(通常,除了高 SNR 信号之外,一阶导数可能会出现问题,因为它们会放大噪声)。

黑线显示了我们将明智地“分割”图像的阈值,这样我们就可以避免整个底噪。如果我们然后应用到我们的原始信号,我们会得到以下图表,黑线表示我们的信号能量的开始:

在此处输入图像描述

因此,这产生了一个个样本。δt=241

我希望这会有所帮助。

使用自相关执行此操作存在一些问题

  1. 巨大的直流偏移(已修复)
  2. 时间窗口:Matlab 的 xcorr() 具有令人讨厌的约定,即当您滑动时间延迟时,在两端基本上“零填充”信号。即数据窗口是时间滞后的函数。这将为任何固定信号(包括正弦波)创建一个三角形。更好的选择是选择一个相关窗口,以便窗口大小加上最大时间滞后适合您的总数据窗口,或者通过非填充样本的数量来标准化互相关。
  3. 这两个信号在我看来并不特别相关。形状有点相似,但峰和谷的具体间距却大不相同,所以我怀疑即使是适当的自相关也能在这里产生很多洞察力。

一个更简单的方法是使用阈值检测器来找到起点,并简单地使用这些点之间的差异作为延迟。

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

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

请注意,这不是一个完美的解决方案,因为大滞后结果现在比低滞后结果更不稳定,因为它们基于较少的数据。四肢的大峰值可能是虚假的,原因与您可以在几次抛硬币中获得 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 个样本的正确滞后处的峰值。当应用于包含正确相似性类型的数据集时,相关性可以成为确定时间延迟的有用工具。