信号之间的时延估计

说明

这篇文章的目的是为了探讨对存在明显时延关系的不同路通道内信号的时延估计。之前在一个声源定位的课题中遇到了相关技术，其中TOA/TDOA-DOA算法的关键在于对不同通道之间的信号时延的估计，因此这里有必要总结一下。

1. 互相关的常见计算方法

在实际应用中，通常使用信号之间的互相关CC（Cross-Correlation）来计算时延，互相关函数的定义如下：
$$
R\left( \tau \right) = E\left[ x_1\left( m \right) \cdot x_2\left( {m + \tau } \right) \right] \tag{1}
$$
其中 ${ {x_1}\left( m \right)}$ 和 ${ {x_2}\left( m \right)}$ 是两个不同的连续信号。而对于经过ADC采样后得到的 ${N}$ 点离散信号 ${ {x_1}\left[ m \right]}$，若假设其对应的时间轴为 ${\left[ {0,N - 1} \right]}$，则上式可以转化为：

$$
\begin{align}
R\left[ n \right] & = E\left[ { {x_1}\left[ m \right] \cdot {x_2}\left[ {m + n} \right]} \right] \\
& = \sum\limits_{m = - n}^{N - 1} { {x_1}\left[ m \right] \cdot {x_2}\left[ {m + n} \right]}
\end{align} \tag{2}
$$

特别注意，信号之间的相关函数可以利用线性卷积来计算，二者的表达式如下所示，这里以连续时间信号为例：
$$
f\left( \tau \right) = x\left( \tau \right) * y\left( \tau \right) = \int\limits_{ - \infty }^{ + \infty } {x\left( t \right) \cdot y\left( {\tau - t} \right)dt} \tag{3}
$$

$$
R\left( \tau \right) = x\left( \tau \right) \star y\left( \tau \right) = \int\limits_{ - \infty }^{ + \infty } {x\left( t \right)^* \cdot y\left( {t + \tau } \right)dt} \tag{4}
$$

其中，计算互相关需要进行取共轭操作，而互相关可以由线性卷积表达出来：

$$
R\left( \tau \right) = \int\limits_{ - \infty }^{ + \infty } {x{ {\left( t \right)}^*} \cdot y\left( {t + \tau } \right)dt} = \int\limits_{ - \infty }^{ + \infty } {x{ {\left( { - t} \right)}^*} \cdot y\left( { - t + \tau } \right)dt} \tag{5}
$$

$$
R\left( \tau \right) = x\left( \tau \right) \star y\left( \tau \right) = R\left( \tau \right) = x\left( -\tau \right)^* * y\left( \tau \right) \tag{6}
$$
即两个信号的互相关等于第一个信号经过翻折共轭后与第二个信号的线性卷积，因此当两个信号的长度分别为 ${M}$ 点和 ${N}$ 点时，最终计算所得互相关函数的长度为 ${(M+N-1)}$ 点。特别的，当两个信号均为 ${N}$ 点长时，其互相关结果为 ${(2*N-1)}$ 点长。

利用线性卷积计算互相关时，可以借助线性卷积和圆周卷积的关系来实现利用频域FFT的方式计算互相关，这种方法相对于传统的时域计算法可以极大地减小计算量，因此在开发底层的互相关计算时经常用到。应该注意到，这种频域计算法实际上也是维纳-辛钦定理的一种体现：由于平稳信号的自相关函数与其功率谱密度互为傅里叶变换对，因此相应的，在频域利用信号的互功率谱可以用来计算时域互相关函数，这也是 (7) 式在频域形式下的物理意义。
$$
P\left( w \right) = F_X^*\left( w \right) \cdot F_Y\left( w \right)\tag{7}
$$

$$
P\left( w \right) = \int\limits_{ - \infty }^{ + \infty } {R\left( \tau \right) \cdot {e^{ - jw\tau } }d\tau }\tag{8}
$$

$$
R\left( \tau \right) = \frac{1}{ {2\pi } }\int\limits_{ - \infty }^{ + \infty } {P\left( w \right) \cdot {e^{jw\tau } }dw}\tag{9}
$$

以下利用MATLAB进行仿真验证，其中选择的原始数据如下所示：

%% 互相关的不同计算方法
clear;
close all;
clc;

% 假定原始信号
x1 = [0,0,1,2,3,7,9,8,0,0];
x2 = [1,2,3,7,9,8,0,0,0,0]; % 可以看到超前2个数据点
N = length(x2);
time = -N+1:N-1;

可以发现，第一路信号相较第一路信号的时延为 -2/point，以下分别利用时域法、频域法以及MATLAB自带的 ${xocrr()}$ 函数计算互相关，代码如下：

% 传统时域法
xcorrTime = zeros(2*N-1,1);
m = 0;
for i = -(N-1):N-1
    m = m+1;
    for t = 1:N
        if 0<(i+t)&&(i+t)<=N
            xcorrTime(m) = xcorrTime(m) + x2(t)*x1(t+i);
        end 
    end
end
xcorrTime = xcorrTime'/N;

% 频域法
Nfft = length(x1)+length(x2)-1;
xcorrfft = fftshift(ifft(fft(x1,Nfft).*conj(fft(x2,Nfft))));

% Matlab调用xcorr计算
xcorrMat = xcorr(x1,x2,'biased');
[k,ind] = max(xcorrMat);

结果如下所示：

Fig 1. 不同方法计算互相关

可以发现三者除了在幅度上存在比例关系外，输出结果的时延关系保持相同。特别注意，MATLAB中的xcorr()函数本身就是利用频域方法编写的，现从xcorr.m中将核心代码摘录如下，以供底层移植参考：

function c = crosscorr(x,y,maxlag)
% Compute cross-correlation for vector inputs. Output is clipped based on
% maxlag but not padded if maxlag >= max(size(x,1),size(y,1)).
X = fft(x,m2,1);
Y = fft(y,m2,1);
if isreal(x) && isreal(y)
    c1 = real(ifft(X.*conj(Y),[],1));
else
    c1 = ifft(X.*conj(Y),[],1);
end
% Keep only the lags we want and move negative lags before positive
% lags.
c = [c1(m2 - mxl + (1:mxl)); c1(1:mxl+1)];

2. GCC-PHAT

理论上使用传统的CC方法已经可以得到时延值，但是考虑到实际的信号会有噪声，且低信噪比会导致互相关函数的峰值不够明显，在判断极值的时候造成误差。为了得到具有更陡峭极值的互相关函数，一般在频域使用一个加权函数来白化输入信号，这就是经典的广义互相关方法（generalized cross-correlation，GCC）。

$$
\widetilde R\left( \tau \right) = \frac{1}{ {2\pi } }\int\limits_{ - \infty }^{ + \infty } {A\left( w \right) \cdot P\left( w \right) \cdot {e^{jw\tau } }dw}\tag{10}
$$

通过查阅资料可得，常见的加权函数如下所示：

Fig 2. 常见的加权函数

这里仅对PHAT加权下的GCC进行分析，此时的加权函数为：

$$
A(w) = { {\rm{1} } \over {\left| { {G_{XY} }\left( w \right)} \right|} } = {1 \over {X\left( w \right) \cdot Y{ {\left( w \right)}^*} } }\tag{11}
$$

理想的加权输出结果为：
$$
R\left( \tau \right) = a \cdot \delta \left( {\tau - {\tau _{XY} } } \right)\tag{12}
$$

利用MATLAB仿真如下所示：

% 频域计算GCC-PHAT
Gss = fft(x1,Nfft).*conj(fft(x2,Nfft));
Gss = Gss./abs(Gss);
% Gss = exp(1i*angle(Gss)); % 这种方式也可以
xcorrGcc = fftshift(ifft(Gss));

与其它方法之间的结果对比如下：

Fig 3. GCC-PHAT计算互相关对比与其它方法对比

可以发现GCC-PHAT结果中仅在期望的时延 2/points 处产生了一个非零值，其它时刻都保持为零，因此更能够凸显出时延的相对关系。以下是对GCC-PHAT算法的理解：

PHAT算法的出发点是，利用CC求解时延时，关心的地方在于两个通道中信号的相位关系，因此利用GCC求解出的互相关结果的幅度没有实际的意义，只需确定时延对应的位置即可。所以PHAT将互功率谱进行了幅度归一化，即在每个频点上幅度分量等价，而相位仍保持原来的信息，这样进行反傅里叶变换后，可以在时延处获得较大的峰值。

从信号平坦性和频谱之间的联系分析：由于时域和频域具有对偶性，当进行PHAT加权后，互功率谱变得更加平坦，同时这也意味着频谱分量相对更加丰富，因此时域上只有出现较为陡峭的冲激才能够与这种“平坦”的信号频谱相对应。

其实，信号进行了PHAT加权后，实质上相当于进行了白化滤波，倘若信号转换为了理想的白噪声形式，则时域中必然会在零延时处出现冲激，这是由傅里叶变换决定的。然而仅进行了幅度上的归一化尚且不能认为是理想的白化，被保留的相位信息仍占据主导因素，因此会在准确的时延处产生更大的冲激。因此两路具有确定延时关系的信号在进行GCC-PHAT后会获得两个冲激峰值。如果两路信号不具备延时关系，则它们的相关性很小，因此进行PHAT加权而产生的频谱更加趋向于理想白噪声，因此它们的时域结果也趋向于在零时延处产生较大峰值。

这样会引入一个模糊问题，即考虑两路信号之间的延时很小或两路信号不存在相关性两种情况，GCC-PHAT在理论上会获得几乎一致的结果，即都仅仅会在零时延附近产生一个峰值，从而难以区分。这是一个值得思考的问题，需要结合实际的信号之间的互相关情况进行分析，这里暂且搁置争议，不再深入讨论。