麦克风均匀线阵（ULA）计算导向矢量（steering vector）

麦克风均匀线阵（ULA）计算导向矢量（steering vector）

本文使用python实现麦克风均匀线阵的导向矢量的计算

公式参考Sharon Gannot 的综述:

[1] S. Gannot, E. Vincent, S. Markovich-Golan, and A. Ozerov, “A Consolidated Perspective on Multimicrophone Speech Enhancement and Source Separation,” IEEE/ACM Trans. Audio Speech Lang. Process., vol. 25, no. 4, pp. 692–730, Apr. 2017, doi: 10.1109/TASLP.2016.2647702.

麦克风均匀线阵如下图（图片来自于[1]）所示：

计算steering vector 的公式如下

\[\mathbf{d}\left( \theta ,\lambda \right) =\left[ 1,e^{-2j\pi \frac{\ell}{\lambda}\cos \left( \theta \right)},...,e^{-2j\pi \left( I-1 \right) \frac{\ell}{\lambda}\cos \left( \theta \right)} \right]^T \]

上式中，\(I\) 为麦克风数量，\(\ell\) 为麦克风间距，\(\lambda=\frac{c}{\nu_{f}}\)为波长，\(c\)为声速，\(\nu_{f}\) 为实际频率，即\(\nu_f=\frac{f}{F}f_s\)，\(f\in\{0,...,F-1\}\)，\(F\)为FFT长度，\(f_s\)为采样频率。

import numpy as np

#                                          |
#                                          o 3
#                                       /  o 2
#                                      / / o 1 mic #
#                                       / /|
#                                        / |
#                                       / θ|

def cal_theta(ref_pos, src_pos):
    if ref_pos[0, 2] == src_pos[0, 2]:
        theta = np.arccos((ref_pos[0, 1]-src_pos[0, 1])/np.linalg.norm(ref_pos-src_pos))
        return theta
    else:
        print('ref_pos\'s z-axis coordinate is not equal to src_pos z-axis coordinate!')


if __name__ == '__main__':

    # this code calculates the degree between reference position and source position                               
    ref_pos = np.array([3, 3, 1.2]).reshape(1, 3)
    src_pos = np.array([1.05, 1.95, 1.2]).reshape(1, 3)
    theta = cal_theta(ref_pos, src_pos)  
    
    theta = 6*np.pi/12
    print(theta*180/np.pi, end='')
    print('°')
    c = 343  # sound speed m/s
    mic_num = 4
    mic_space = 0.08 
    fs = 16000
    fft_size = 1024
    SV_mat = np.zeros((mic_num, fft_size//2 + 1), dtype=complex)  #
    for f in range(int(fft_size/2 + 1)):
        vf = f / fft_size * fs
        if vf == 0:
            for m in range(mic_num):
                SV_mat[m, f] = np.exp(-2j * np.pi * m * 0 * np.cos(theta))
        else:
            lamb = c / vf
            for m in range(mic_num):
                SV_mat[m, f] = np.exp(-2j*np.pi*m*mic_space/lamb*np.cos(theta))
    # SV_mat = SV_mat/np.linalg.norm(SV_mat, axis=0) # normalize