时延估计之互相关方法(Cross-Correlation Method, CC)

时延估计之互相关方法(Cross-Correlation Method, CC)

信号模型

在一个声学环境中，假如声场中只有一个麦克风，则传声器的输出信号可以写为：
$$x_n\left(k\right)={\alpha }_{n}s\left(k-D_n\right)+b_n\left(k\right),$$
其中${\alpha }_n$满足$0<{\alpha }_n<1$为一个衰减因子（由于声传播过程中的能量衰减），$D_n$为从声源到麦克风之间的传播时间，$s(k)$为语音信号，通常被认为是宽带信号，$b_n(k)$为高斯随机噪声，通常认为声源信号与噪声信号是不相关的。在这个模型下，第$i$个声源与$j$个声源之间的到达时间差(TDOA)可以写作：
$${\tau }_{ij}=D_j-D_i,$$
其中$i,j=1,2,\cdots ,N,$且$i\ne j$，而时延估计(TDE)的目的就是得到${\tau }_{ij}$的估计$\widehat{{\tau }_{ij}}$。

互相关方法(Cross-Correlation Method)

假如现在有两个传感器，每个传感器的输出分别用$x_1(k)$和$x_2(k)$表示，则两个信号之间的互相关函数可以写作：
$$r_{x_1{x}_2}\left(p\right)=E\left[x_1\left(k\right)x_2\left(k+p\right)\right].$$
则由上面的第一个式子，$x_1(k)$和$x_2(k)$表示为：
$$x_1\left(k\right)={\alpha }_{1}s\left(k-D_1\right)+b_1\left(k\right),$$
$$x_2\left(k+p\right)={\alpha }_{2}s\left(k-D_2+p\right)+b_2\left(k+p\right),$$
$$\begin{gathered} r_{x_1{x}_2}\left(p\right)=E\left[\left[{\alpha }_{1}s\left(k-D_1\right)+b_1\left(k\right)\right]\left[{\alpha }_{2}s\left(k-D_2+p\right)+b_2\left(k+p\right)\right]\right] \\ =E\left[{\alpha }_1{\alpha }_{2}s\left(k-D_1\right)s\left(k-D_2+p\right)+{\alpha }_{1}s\left(k-D_1\right)b_2\left(k+p\right)+{\alpha }_2{b}_1\left(k\right)s\left(k-D_2+p\right)+b_1\left(k\right)b_2\left(k+p\right)\right] \\ ={\alpha }_1{\alpha }_2{r}_{ss}\left(p+D_1-D_2\right)+{\alpha }_1{r}_{s{b}_2}\left(p+D_1\right)+{\alpha }_2{r}_{s{b}_1}\left(p-D_2\right)+r_{b_1{b}_2}\left(p\right). \end{gathered}$$
通过上面的假设：$b_n(k)$是高斯序列，且与源信号和其他传感器上的噪声不相关，则上式可以简化为：
$$r_{x_1{x}_2}\left(p\right)={\alpha }_1{\alpha }_2{r}_{ss}\left(p+D_1-D_2\right),$$
当$p=D_2-D_1$ 时，互相关函数$r_{x_1{x}_2}$达到最大，所以通过互相关函数，我们可以得到$x_1(k)$与$x_2(k)$之间的TDOA为：
$${\hat{\tau }}_{12}=\arg \underset{p}{\max }{r}_{x_1{x}_2}\left(p\right),$$
其中$p\in [-{\tau }_{\max },{\tau }_{\max }]$.

在实际的应用中我们可以通过借助FFT等高效计算相关函数：
$${\hat{r}}_{x_1{x}_2}\left(p\right)=\frac{1}{K}\sum _{k^{\prime }=0}^{K-1}{X}_1\left({\omega }_{k^{\prime }}\right)X_2^{\ast }\left({\omega }_{k^{\prime }}\right)e^{j{\omega }_{k^{\prime }}p},$$

Matlab代码

互相关方法函数的代码：

function [tau] = CC_function(X_mat,block_size)
%CC_FUNCTION : Cross-Correlation method for time delay estimation
%   input:
%           Xmat: multi-channel signal matrix,
%           block_size: snapshot length
%   output:
%           tau: TDE parameters 
M = size(X_mat,2);
data_length = size(X_mat,1);
% block_size = 2048;
x1_vect = X_mat(1:block_size,1);
x2_vect = X_mat(1:block_size,2);
%% Method.1 Using Matlab XCORR funciotn
% CCF = xcorr(x1_vect,x2_vect);
% half_CCF = CCF(ceil(length(CCF)/2)+1:end); %% remove r_xy(<=0)


%% Method.2 Using FFT
x1_fft = fft(x1_vect);
x2_fft = fft(x2_vect);
CCF = ifft((x1_fft.*conj(x2_fft)));
half_CCF = CCF(2:end); %% remove r_xy(0) 

%% time delay parameter
[value,tau] = max(half_CCF);
end

这里用了两种方法来计算互相关，一种是Matlab自带的xcorr函数，另一种是使用FFT来计算互相关函数。

测试代码：

clc
clear

% 读取音频文件
[h,fs] = audioread('yourpath\RWCP_type1_rir_circle_ane_imp320.wav');

% 对h的第一列求绝对值最大值和对应的索引
[value1,D_1] = max(abs(h(:,1)));

% 对h的第四列求绝对值最大值和对应的索引
[value2,D_2] = max(abs(h(:,4)));

% 计算实际的时间差
real_TDE = abs(D_2 - D_1);

% 读取语音文件
[x,fs] = audioread('yourpath\your speech file.wav',[1,fs]);

% 将x与h的第一列进行卷积
x1_vect = conv(x,h(:,1));

% 将x与h的第四列进行卷积
x2_vect = conv(x,h(:,4));

% 构建X矩阵，包含x1_vect和x2_vect
X_mat = [x1_vect, x2_vect];

% 定义块大小为1024
block_size = 1024;

% 调用CC_function函数计算估计的时间差
[estimated_TDE] = CC_function(X_mat, block_size);

测试结果