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

波束形成（BF）从算法仿真到工程源码实现-第六节-广义旁瓣消除算法(GSC)

news 来源：原创 2025/5/18 0:50:10

一、概述

本节我们讨论广义旁瓣消除算法(GSC)，包括原理分析及代码实现。更多资料和代码可以进入https://t.zsxq.com/qgmoN ，同时欢迎大家提出宝贵的建议，以共同探讨学习。

二、原理分析

广义旁瓣消除(GSC)算法

GSC算法是与LCMV算法等效的，其权矢量被分解为自适应部分和非自适应部分，其中自适应部分正交于约束子空间，而非自适应部分位于约束子空间内，其权矢量可以表示为

$w = w_{q}-Bw_{a}$

其中， $w_{q}=(CC^{H})^{-1}Cf$ ， $w_{a}=(B^{H}R_{x}B)^{-1}B^{H}R_{x}w_{q}$ 。 $B$ 为阻塞矩阵，正交于约束矩阵 $C$ ，其作用是为了阻止期望信号进入辅助支路。关于 $B$ 可以通过求 $C$ 的补空间来确定

$B = I - C(C^{H}C)^{-1}C^{H}$

主支路的输出 $y = w_{q}^{H}x$ ，阻塞矩阵投影后的输出为 $z = B^{H}x$ ,那么自适应的权矢量可以表示为

$w_a = R_z^{-1}P_z$

其中， $R_z = B^{H}R_xB$ 是 $z$ 的协方差矩阵， $P_z = B^{H}R_xw_q$ 是 $z$ 和 $y$ 的互协方差矩阵。GSC是LCMV的等效，其将后者的有约束的优化问题变成了无约束的优化问题，当 $z$ 中含有较少期望信号时，GSC还能正常工作，反之，其性能会大幅度下降。1999年，O.Hoshuyama等人采用约束自适应滤波的方法代替原来的对齐相减，以及采用当期望语音存在时只更新阻塞矩阵，而当期望语音不存在的时候只更新自适应抵消器的系数来减小期望语音信号的泄露。

三、代码仿真

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 gsc(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))
        wq_i = np.matmul(C_i, np.conjugate(C_i).transpose())
        wq_i = np.linalg.pinv(wq_i)
        wq_i = np.matmul(wq_i, C_i)
        wq_i = np.matmul(wq_i, d)

        B_i = np.matmul(np.conjugate(C_i).transpose(), C_i)
        B_i = np.linalg.pinv(B_i)
        B_i = np.matmul(C_i, B_i)
        B_i = np.matmul(B_i, np.conjugate(C_i).transpose())
        B_i = np.eye(6) - B_i

        Rz_i = np.matmul(np.conjugate(B_i).transpose(), np.reshape(Rxx[:,:], (6, 6)))
        Rz_i = np.matmul(Rz_i, B_i)

        Pz_i = np.matmul(np.conjugate(B_i).transpose(), np.reshape(Rxx[:,:], (6, 6)))
        Pz_i = np.matmul(Pz_i, wq_i)

        wa_i = np.matmul(np.linalg.pinv(Rz_i), Pz_i)

        w = wq_i - np.matmul(B_i, wa_i)
        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 = gsc(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_gsc_t0_i180.wav", frames1, 16000)

main()