基于被动时反的 ac-MVDR 水下声源定位python代码

以下是基于被动时反的 ac-MVDR 水下声源定位的 Python 代码：

import numpy as np
import scipy.signal as signal

def ac_mvdr_doa_estimation(data, freq_bins, array_geometry, sound_speed, num_src=None):
    """
    基于被动时反的 ac-MVDR 水下声源定位
    :param data: 录音数据，shape=(num_channels, num_samples)
    :param freq_bins: 频率点，shape=(num_freq_bins,)
    :param array_geometry: 阵列几何，shape=(num_channels, 3)，第一维是x坐标，第二维是y坐标，第三维是z坐标
    :param sound_speed: 声速
    :param num_src: 声源个数，默认为None，表示不知道声源个数，需要进行自动估计
    :return: DOA估计结果，shape=(num_src,)
    """
    num_channels = array_geometry.shape[0]

    # 计算每个传感器之间的距离
    d = np.zeros((num_channels, num_channels))
    for i in range(num_channels):
        for j in range(num_channels):
            d[i, j] = np.linalg.norm(array_geometry[i] - array_geometry[j])

    # 计算传感器间相对延迟
    tau = d / sound_speed
    tau -= np.min(tau)

    # STFT
    window = signal.hamming(256)
    hop_size = 128
    f, t, spec = signal.stft(data, window=window, nperseg=256, noverlap=128)
    spec = np.transpose(spec, (2, 0, 1))  # shape=(num_freq_bins, num_channels, num_frames)

    # 声源个数估计
    if num_src is None:
        num_frames = spec.shape[2]
        P = np.zeros((num_freq_bins, num_channels, num_channels), dtype=np.complex64)
        for i in range(num_freq_bins):
            Rxx = np.dot(spec[i], np.conj(spec[i]).T) / num_frames
            P[i] = np.linalg.inv(Rxx)
        P_sum = np.sum(P, axis=0)
        eig_values, _ = np.linalg.eig(P_sum)
        lambda_min = np.min(np.real(eig_values))
        num_src = num_channels - np.sum(np.real(eig_values) < lambda_min / 10)

    # 计算空间谱
    spec_2 = np.abs(spec) ** 2
    Ps = np.zeros((num_freq_bins, num_channels, num_channels), dtype=np.complex64)
    for i in range(num_freq_bins):
        Rxx = np.dot(spec[i], np.conj(spec[i]).T) / spec.shape[2]
        Ps[i] = np.linalg.inv(Rxx)

    # 计算空间协方差矩阵的逆
    P_sum = np.sum(Ps, axis=0)
    inv_P_sum = np.linalg.inv(P_sum)

    # 计算指纹矩阵
    X = np.zeros((num_freq_bins, num_channels, num_src), dtype=np.complex64)
    for i in range(num_freq_bins):
        U, s, V = np.linalg.svd(Ps[i])
        U = U[:, :num_src]
        X[i] = np.dot(np.conj(U).T, inv_P_sum)

    # 计算DOA
    theta_range = np.arange(-90, 90, 0.5)
    num_theta = len(theta_range)
    P_music = np.zeros((num_freq_bins, num_theta))
    for i in range(num_freq_bins):
        Xf = X[i]
        for j in range(num_theta):
            steering_vec = np.exp(-1j * 2 * np.pi * freq_bins[i] * tau *
                                  np.sin(theta_range[j] / 180 * np.pi))
            P_music[i, j] = 1 / (np.dot(np.conj(steering_vec).T, np.dot(Xf, Xf.conj().T.dot(steering_vec))))

    # DOA估计
    doa_estimates = np.zeros(num_src)
    for i in range(num_src):
        theta_idx = np.argmax(np.sum(P_music[:, :num_theta // 2], axis=1))
        doa_estimates[i] = theta_range[theta_idx]
        P_music[:, theta_idx - 2:theta_idx + 3] = 0  # 降低峰值，避免重复估计同一声源

    return doa_estimates

以上代码使用了一些依赖库，如numpy、scipy等，需要提前安装。代码中的注释可以帮助你更好地理解每个步骤的计算过程。在使用时，你需要准备好录音数据、频率点、阵列几何和声速等参数，然后调用`ac_mvdr_doa_estimation`函数即可得到DOA估计结果。