TDE （Time Delay Estimation）时延估计

最近正在学习AEC，需要完成线性回声消除的一些简单的代码，第一个问题是需要求出输入信号与参考信号之间的时延。因为时延影响了后面AEC算法的效果，所以对整体的回声消除效果比较重要。

一、用互相关的方法完成时延估计

The generalized cross-correlation (GCC) is regarded as the most popular approach for estimating the time difference of arrival (TDOA) between the signals received at two sensors. The GCC of a pair of sensor signals is defined as the inverse Fourier transform of the weighted cross-power spectrum

上面的图片来自：Frequency-Sliding Generalized Cross-Correlation: A Sub-band Time Delay Estimation Approach ----Maximo Cobos, Senior Member, IEEE, Fabio Antonacci, Senior Member, IEEE, Luca Comanducci, Student Member, IEEE, and Augusto Sarti, Senior Member, IEEE

 

GCC TDE用FFT计算的流程：

 

 

互相关法：
                

广义互相关法：延时信息体现在互功率谱的相位上，而与幅度无关。

                

这种互功率谱加权方式，称为GCC-PHAT（PHAse Transformation），或
CSP（Cross Spectral Phase）。

             

百度百科对互相关的解析：  

        互相关函数是信号分析里的概念，表示的是两个时间序列之间的相关程度，即描述信号x(t),y(t)在任意两个不同时刻t1，t2的取值之间的相关程度。描述两个不同的信号之间的相关性时，这两个信号可以是随机信号，也可以是确知信号。

当计算出来的值等 于 1 的时候，相当于两个信号完美重合，因此可以利用这个信息来求两个信号的延时。

 

互相关 的 矩阵计算:

可以用来计算频域里的互相关

For two random vectors and  , each containing random elements whose expected value and variance exist, the cross-correlation matrix of X and Y  is defined by

and has dimensions .

Written component-wise:

The random vectors X and Y need not have the same dimension, and either might be a scalar value.

 

看完上面对互相关的定义和计算方法，可以看出，要求出两个信号的延迟。可以给定一段时长的信号（例如1S），然后做FFT来计算广义互相关，最后再做傅里叶逆变换找出最大的一个Rxy 的位置， 用这个位置来确定延时，精确到sample.

clear all;

y_fn = '..\audio_data\input.raw'
y_fid = fopen(y_fn,'r');
y = fread(y_fid, 'short', 'ieee-le');
y = y(1:4800) ;
fclose(y_fid);
%sound(y, 16000);

x_fn = '..\audio_data\ref.raw'
x_fid = fopen(x_fn,'r');
x = fread(x_fid, 'short', 'ieee-le');
x = x(1:4800) ;
fclose(x_fid);
%sound(x, 16000);

y_fd = fft(y);
x_fd = fft(x);


P_yx = (y_fd .* conj(x_fd));
P_yx_abs = abs(y_fd .* conj(x_fd));

%R_yx = ifft(P_yx);             %cross-corelation
R_yx = ifft(P_yx ./ P_yx_abs);  %generalized cross-corelation

R_yx = R_yx(1:4800);
[R_yx_max,m_idx] =  max(R_yx);
stem(1:4800, R_yx)
disp(m_idx);

当然matlab 里面已经有提供了gccphat的工具代码，可以直接食用