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

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

news 来源：原创 2025/9/6 19:42:44

一、概述

本节我们讨论广义旁瓣消除算法(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 pltfft_size = 256
freq_bin = 129def 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_sizeif frequency == 0:phase_delay = 0steering_vector[n, f] = np.exp(1j * phase_delay)else:lambda_val = c / frequencytheta_mic = -2 * np.pi * n / N + 2 * np.pitheta_signal = np.pi * angle / 180phase_delay = 2 * np.pi * np.cos(theta_signal - theta_mic) * r / lambda_valsteering_vector[n, f] = np.exp(1j*phase_delay)return steering_vectordef 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_sizeif frequency == 0:phase_delay = 0steering_vector[n, f] = np.exp(1j * phase_delay)else:lambda_val = c / frequencytheta_mic = 2 * np.pi * n / Ntheta_signal = np.pi * angle / 180phase_delay = 2 * np.pi * np.cos(theta_signal - theta_mic) * r / lambda_valsteering_vector[n, f] = np.exp(1j * phase_delay)return steering_vectordef 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_iRz_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) / 6return data2def 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)#lcmvC = 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_windowfft_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_frame1fft_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)] + out1out1 = ifft_frame1[int(fft_size/2):]frames.extend(short_data1)frames1 = np.array(frames).reshape((-1)).realsf.write("output/simulate_role1_0_t60_0.2_role2_180_t60_0.2_out_gsc_t0_i180.wav", frames1, 16000)main()

四、仿真结果

4.1 0度方向为期望信号，180度为干扰方向

4.2 180度方向为期望信号，0度方向为干扰方向

五、总结

对比发现GSC和LCMV的结果一致，可以验证GSC是LCMV的等效形式。自适应滤波器形式可以参考athena中的代码。