波束形成（BF）从算法仿真到工程源码实现-第四节-最小方差无失真响应波束形成（MVDR）

一、概述

        本节我们讨论最 小 方 差 无 失 真 响 应 (Minimum Variance Distortionless Response, MVDR)波束形成算法，包括原理分析及代码实现。 更多资料和代码可以进入https://t.zsxq.com/qgmoN ，同时欢迎大家提出宝贵的建议，以共同探讨学习。

二 、原理分析

        1969年，J.Capon提出了最小方差无失真响应(Minimum Variance Distortionless Response, MVDR)波束形成算法。该算法是应用得最为广泛的自适应波束形成方法之一。

        原理:在期望信号无失真的的束条件下，选择合适的滤波器系数，使得阵列输出的平均功率最小化。MVDR的权重优化问题可以表示为

                ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​ 

其中，为目标信号导向矢量，表示声源方向和麦克风之间的传递函数，可以通过纯净语音信号达到每个麦克风的不同延迟时间计算得到。为空间信号相关协方差矩阵，根据快拍次数估计得到。当在时间上彼此不相关的个噪声信号从不同方向到达麦克风麦克风时，空间相关协方差矩阵被定义为:
        ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​                

运用拉格朗日乘子法计算得到最优权重为:

        ​​​​​​​        ​​​​​​​        ​​​​​​​                        

因为约束条件是纯净语音信号无失真，即纯净语音信号是保持不变的，为了使得输出的方差最小化，仅仅只需要让噪声信号最小化。所以上式信号相关矩阵,可以用噪声相关矩阵替换。

        当为非相干噪声场的相关矩阵，则MVDR退化为延迟求和波束形成器，当
为散射噪声场的相关矩阵，则MVDR退化为超指向性波束形成器。

        MVDR是普遍采用的波束形成典型算法，在复杂环境下，由于协方差矩阵计算的不精确性算法会导致性能急剧下降。后来的研究者提出了许多基于对角加载的解决方法。这些方法解决了对角加载值不易确定且无法通过样本更改自动调整的问题，使协方差矩阵误差问题得到一定的改善。但是，这些算法相对比较复杂，效率较低。

三、代码仿真

import numpy as np
import soundfile as sf
import scipy
import matplotlib.pyplot as plt


fft_size = 256
freq_bin = 129

def calculate_circular_array_steering_vector(angle, r=0.0463, N=6, fs=16000, fft_size=256, c=343):
    steering_vector = np.zeros((N, fft_size//2 + 1), dtype=complex)
    for f in range(int(fft_size/2+1)):
        for n in range(N):
            frequency = fs * f / fft_size
            if frequency == 0:
                phase_delay = 0
                steering_vector[n, f] = np.exp(1j * phase_delay)
            else:
                lambda_val = c / frequency
                theta_mic = -2 * np.pi * n / N + 2 * np.pi
                theta_signal = np.pi * angle / 180
                phase_delay = 2 * np.pi * np.cos(theta_signal - theta_mic) * r / lambda_val
                steering_vector[n, f] = np.exp(1j*phase_delay)

    return steering_vector


def calculate_circular_array_steering_vector_anticlockwise(angle, r=0.0463, N=6, fs=16000, fft_size=256, c=343):
    steering_vector = np.zeros((N, fft_size // 2 + 1), dtype=complex)
    for f in range(int(fft_size / 2 + 1)):
        for n in range(N):
            frequency = fs * f / fft_size
            if frequency == 0:
                phase_delay = 0
                steering_vector[n, f] = np.exp(1j * phase_delay)
            else:
                lambda_val = c / frequency
                theta_mic = 2 * np.pi * n / N
                theta_signal = np.pi * angle / 180
                phase_delay = 2 * np.pi * np.cos(theta_signal - theta_mic) * r / lambda_val
                steering_vector[n, f] = np.exp(1j * phase_delay)

    return steering_vector

def mvdr(Rxx, a, data):
    beamformer = np.zeros((6, freq_bin), dtype=complex)
    for i in range(freq_bin):
        Rxx_i = Rxx
        # Rxx_i = np.reshape(Rxx[i, :, :], (6, 6))
        a_i = a[:, i]
        Rxx_i_inv = np.linalg.pinv(Rxx_i)
        fenzi = np.matmul(Rxx_i_inv, a_i)
        fenmu = np.matmul(np.conjugate(a_i).transpose(), Rxx_i_inv)
        fenmu = np.matmul(fenmu, a_i) + 1e-6
        beamformer[:, i] = fenzi / fenmu
    data1 = np.multiply(np.conjugate(beamformer), data)
    data2 = np.sum(data1, axis=0) / 6
    return data2
def main():
    # 读取WAV文件
    data, samplerate = sf.read('output/simulate_role1_0_t60_0.2_role2_180_t60_0.2.wav')

    # 定义帧长和帧移
    frame_length = int(samplerate * 0.016)  # 25ms帧长
    frame_step = int(samplerate * 0.008)  # 10ms帧移

    # 创建汉明窗
    hamming_window = scipy.signal.windows.hamming(frame_length)
    hamming_window = np.reshape(hamming_window, [frame_length, 1])

    sample_num = data.shape[0] - frame_length + 1

    HH = calculate_circular_array_steering_vector(180)

    # 手动分帧和加窗
    frames = []
    out1 = np.zeros(int(fft_size/2), dtype=float)

    for i in range(0, sample_num, frame_step):
        frame = data[i:i + frame_length, :]
        windowed_frame = frame * hamming_window
        fft_frame = np.fft.fft(windowed_frame, axis=0)
        fft_frame1 = np.transpose(fft_frame[:freq_bin, :])

        Rxx_frame_real = np.matmul(fft_frame1, np.conjugate(fft_frame1).transpose()) / 129 + 1e-6 * np.eye(6)

        fft_frame1 = mvdr(Rxx_frame_real, HH, fft_frame1)

        fft_frame11 = fft_frame1
        fft_frame21 = np.concatenate((fft_frame11, fft_frame11[1:-1][::-1].conj()))
        fft_frame21 = np.transpose(fft_frame21)
        ifft_frame1 = np.fft.ifft(fft_frame21)
        short_data1 = ifft_frame1[:int(fft_size/2)] + out1
        out1 = ifft_frame1[int(fft_size/2):]
        frames.extend(short_data1)

    frames1 = np.array(frames).reshape((-1)).real
    sf.write("output/simulate_role1_0_t60_0.2_role2_180_t60_0.2_out_mvdr_t180.wav", frames1, 16000)

main()

四、结果展示

4.1 0度方向

4.2 180度方向

五、总结

        从结果上看，我们做0度方向的波束时，role2方向的声音信号被抑制，同时role1方向的声音信号保持不变，做180度方向的波束时，role1方向的声音信号被抑制，同时role2方向的声音信号保持不变，达到了预期效果。