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

波束形成（BF）从算法仿真到工程源码实现-第九节-延迟相减波束形成（delay sub）

news 来源：原创 2025/4/27 10:08:48

一、概述

本节对delay sub算法进行仿真。更多资料和代码可以进入https://t.zsxq.com/qgmoN ，同时欢迎大家提出宝贵的建议，以共同探讨学习。

二、代码仿真

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 delay_sub(a, data):
    data1 = np.multiply(np.conjugate(a), data)
    data2 = np.sum(data1, axis=0) / 6

    result = np.zeros((freq_bin,), dtype=complex)
    for i in range(freq_bin):
        data_i = data1[:, i]
        data_ds = data2[i]
        for ch in range(5):
            result[i] += data_i[ch+1] - data_i[ch]
        result[i] /= 5

    return result




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, :])


        fft_frame1 = delay_sub(HH, fft_frame1)

        #1mic
        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_delaysub_t0.wav", frames1, 16000)

main()