当前位置：首页 > news >正文

波束形成（BF）从算法仿真到工程源码实现-第五节-线性约束最小方差波束形成算法(LCMV)

news 来源：原创 2025/6/7 3:40:33

一、概述

本节我们讨论线性约束最小方差波束形成算法(Linearly constrained minimum variance,LCMV)波束形成算法，包括原理分析及代码实现。更多资料和代码可以进入https://t.zsxq.com/qgmoN ，同时欢迎大家提出宝贵的建议，以共同探讨学习。

二、原理分析

1972年，L.Frost提出了线性约束最小方差(Linearly Coastrained Minimum Variance,LCMV)波束形成器。LCMV波束形成在效果上实际是MVDR波束形成的扩展形式，它将后者中期望信号不受影响的这一约束扩展为一组约束，即为目标方向无失真同时对其它噪声干扰方向陷零。随后L.Frost基于约束最小均方自适应滤波器提出了LCMV算法的自适应结构。

原理:在满足一组约束的同时，使波束形成输出(干扰信号、噪声)的功率最小化。即

$\begin{cases} \min_{w} w^{H}R_{x}w\\ \text{s.t. } C^{H}w = f \end{cases}$

其中， $f=\begin{bmatrix}1\\0\\\vdots\\0\end{bmatrix}^T$ 为Nx1的约束值向量， $C = \begin{bmatrix}a(\theta_{01})&a(\theta_{02})&\cdots&a(\theta_{0N})\end{bmatrix}$ 为MxN维的约束矩阵， $\theta_{0n}$ ,n=1,2 ,... , N 为可能的期望信号方向。 $a(\theta_{0n})$ 为对应的导向矢量。 $R_{x} = E\left[x(t)x^{H}(t)\right]$ 表示输出协方差矩阵。

通过拉格朗日乘数法，可以求解得到最终的权系数为

$w = R_{x}^{-1}C\left(C^{H}R_{x}^{-1}C\right)^{-1}f$

当LCMV方法的约束条件取 $w^{H} a(\theta)=1$ 时，演变为最小方差无失真响应波束形成器(MVDR，minimum variance distortionless-response)。也就是说MVDR算法是LCMV算法的一个特例。其原理是在阵列输出信号能量保持不变的约束条件下，通过调节权重系数使阵列信号输出总功率(相关功率与非相关功率之和)达到最小，由于目标信号的强度得以保持，而噪声的方差被最小化，可以说MVDR使阵列输出信号的信噪比(SNR)达到最大。

三、代码仿真

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 lcmv(C, d, Rxx, data):

    beamformer = np.zeros((6, freq_bin), dtype=complex)
    for i in range(freq_bin):
        C_i = np.reshape(C[i, :, :], (6, 2))
        w_mid = np.linalg.pinv(np.matmul(np.matmul(np.conjugate(C_i).transpose(), np.linalg.pinv(Rxx)), C_i)) # 2*6 * 6*6 * 6*2 = 2*2
        w_front = np.matmul(np.linalg.pinv(Rxx), C_i) # 6*6 * 6*2 = 6*2
        w = np.matmul(w_front, w_mid)  # 6*2 * 2*2 = 6*2
        w = np.matmul(w, d) # 6*2 * 2*1 = 6*1
        beamformer[:, i] = w.reshape(6,)
    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

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

    #lcmv
    C = np.zeros((freq_bin, 6, 2), dtype=complex)
    d = np.reshape(np.array([1, 0]), (1, 2)).transpose()

    desire = calculate_circular_array_steering_vector(0)
    interf = calculate_circular_array_steering_vector(180)
    C[:, :, 0] = np.transpose(desire)
    C[:, :, 1] = np.transpose(interf)

    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 = lcmv(C, d, Rxx_frame_real, 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_lcmv_t0_i180.wav", frames1, 16000)

main()