""" MUSIC 算法核心实现 基于博客：https://blog.csdn.net/weixin_/article/details/

实现步骤：

1. 信号生成与预处理（采样、加噪）
2. 阵列流形矩阵构建
3. 协方差矩阵计算
4. 特征分解提取噪声子空间
5. 空间谱估计与谱峰搜索 """

import numpy as np from numpy.linalg import eig, svd

def generate_signal(signal_angles, snapshots, snr, num_elements, fs=1000):

""" 生成多源信号并添加高斯白噪声 参数: signal_angles: 信号源到达角度列表 (度) snapshots: 快拍数 snr: 信噪比 (dB) num_elements: 阵元数量 fs: 采样频率 返回: X: 接收信号矩阵 (num_elements x snapshots) """ num_signals = len(signal_angles) # 生成时间向量 t = np.arange(snapshots) / fs # 生成基带信号 (复指数信号) signals = np.zeros((num_signals, snapshots), dtype=complex) for i, angle in enumerate(signal_angles): # 每个信号源使用不同的频率以避免相干 freq = 100 + i * 50 # Hz signals[i, :] = np.exp(1j * 2 * np.pi * freq * t) # 构建阵列流形矩阵 A = build_array_manifold(signal_angles, num_elements) # 无噪声接收信号 X_clean = A @ signals # 计算信号功率 signal_power = np.mean(np.abs(X_clean)2) # 根据 SNR 计算噪声功率 # SNR(dB) = 10 * log10(signal_power / noise_power) snr_linear = 10 (snr / 10) noise_power = signal_power / snr_linear # 生成复高斯白噪声 noise = np.sqrt(noise_power / 2) * (np.random.randn(num_elements, snapshots) + 1j * np.random.randn(num_elements, snapshots)) # 接收信号 = 信号 + 噪声 X = X_clean + noise return X 

def build_array_manifold(signal_angles, num_elements, d_lambda=0.5):

""" 构建均匀线阵的阵列流形矩阵 参数: signal_angles: 信号源角度列表 (度) num_elements: 阵元数量 d_lambda: 阵元间距与波长的比值 (通常为 0.5) 返回: A: 阵列流形矩阵 (num_elements x num_signals) """ num_signals = len(signal_angles) A = np.zeros((num_elements, num_signals), dtype=complex) # 角度转换为弧度 angles_rad = np.radians(signal_angles) for i, theta in enumerate(angles_rad): # 方向向量 a(theta) # a(theta) = [1, exp(-j*pi*cos(theta)), exp(-j*2*pi*cos(theta)), ...]^T phase_shift = -np.pi * np.arange(num_elements) * d_lambda * np.cos(theta) A[:, i] = np.exp(1j * phase_shift) return A 

def compute_covariance_matrix(X, method=‘direct’):

""" 计算接收信号的协方差矩阵 参数: X: 接收信号矩阵 (num_elements x snapshots) method: 计算方法 ('direct' 或 'forward_backward') 返回: R: 协方差矩阵 (num_elements x num_elements) """ num_elements, snapshots = X.shape if method == 'direct': # 直接样本协方差估计 R = (X @ X.conj().T) / snapshots elif method == 'forward_backward': # 前向 - 后向空间平滑 # 前向协方差 R_forward = (X @ X.conj().T) / snapshots # 后向协方差 J = np.fliplr(np.eye(num_elements)) # 交换矩阵 X_backward = J @ np.conj(X) R_backward = (X_backward @ X_backward.conj().T) / snapshots # 平均 R = (R_forward + R_backward) / 2 return R 

def extract_noise_subspace(R, num_signals):

""" 通过特征分解提取噪声子空间 参数: R: 协方差矩阵 (num_elements x num_elements) num_signals: 信号源数量 返回: Vn: 噪声子空间 (num_elements x (num_elements - num_signals)) Vs: 信号子空间 (num_elements x num_signals) eigenvalues: 特征值 (降序排列) """ num_elements = R.shape[0] # 特征分解 eigenvalues, eigenvectors = eig(R) # 特征值降序排列 idx = np.argsort(np.real(eigenvalues))[::-1] eigenvalues = np.real(eigenvalues[idx]) eigenvectors = eigenvectors[:, idx] # 分离信号子空间和噪声子空间 # 大特征值对应信号子空间，小特征值对应噪声子空间 Vs = eigenvectors[:, :num_signals] Vn = eigenvectors[:, num_signals:] return Vn, Vs, eigenvalues 

def music_spectrum(Vn, scan_angles, num_elements, d_lambda=0.5):

""" 计算 MUSIC 空间谱 参数: Vn: 噪声子空间 scan_angles: 扫描角度范围 (度) num_elements: 阵元数量 d_lambda: 阵元间距与波长比值 返回: spectrum: 空间谱值 angles: 扫描角度 """ angles = np.array(scan_angles) num_scan_points = len(angles) spectrum = np.zeros(num_scan_points) for i, theta_deg in enumerate(angles): theta_rad = np.radians(theta_deg) # 构造测试方向向量 phase_shift = -np.pi * np.arange(num_elements) * d_lambda * np.cos(theta_rad) a_theta = np.exp(1j * phase_shift).reshape(-1, 1) # MUSIC 谱函数：P(θ) = 1 / (a^H * Vn * Vn^H * a) # 其中 Vn * Vn^H 是噪声子空间的投影矩阵 Pn = Vn @ Vn.conj().T # 噪声子空间投影矩阵 denominator = np.abs(a_theta.conj().T @ Pn @ a_theta)[0, 0] # 避免除零 if denominator < 1e-10: denominator = 1e-10 spectrum[i] = 1.0 / denominator # 归一化到 dB spectrum_db = 10 * np.log10(spectrum / np.max(spectrum) + 1e-10) return spectrum_db, angles 

def peak_search(spectrum, angles, num_peaks=None, threshold_db=-20):

""" 搜索空间谱峰值，估计 DOA 参数: spectrum: 空间谱值 (dB) angles: 对应的角度 num_peaks: 要检测的峰值数量 (None 表示自动检测) threshold_db: 峰值检测阈值 (相对于最大值的 dB 数) 返回: estimated_angles: 估计的到达角度 peak_values: 峰值强度 """ # 寻找局部极大值 peaks = [] peak_values = [] max_spectrum = np.max(spectrum) threshold = max_spectrum + threshold_db for i in range(1, len(spectrum) - 1): if (spectrum[i] > spectrum[i-1] and spectrum[i] > spectrum[i+1] and spectrum[i] > threshold): peaks.append(angles[i]) peak_values.append(spectrum[i]) # 如果没有找到峰值，降低阈值 if len(peaks) == 0: threshold = max_spectrum - 10 for i in range(1, len(spectrum) - 1): if spectrum[i] > spectrum[i-1] and spectrum[i] > spectrum[i+1]: peaks.append(angles[i]) peak_values.append(spectrum[i]) # 按峰值强度排序 if len(peaks) > 0: sorted_idx = np.argsort(peak_values)[::-1] peaks = [peaks[i] for i in sorted_idx] peak_values = [peak_values[i] for i in sorted_idx] # 限制峰值数量 if num_peaks is not None and len(peaks) > num_peaks: peaks = peaks[:num_peaks] peak_values = peak_values[:num_peaks] return peaks, peak_values 

def estimate_num_signals(eigenvalues, num_elements):

""" 基于特征值分布估计信号源数量 使用 AIC 或 MDL 准则 参数: eigenvalues: 特征值 (降序排列) num_elements: 阵元数量 返回: num_signals: 估计的信号源数量 """ # 简单方法：寻找特征值的"拐点" # 信号特征值较大且下降快，噪声特征值较小且接近 # 计算特征值差分 diff = np.diff(eigenvalues) # 寻找最大的下降 if len(diff) > 0: # 归一化差分