一、 项目概览

目标1：脱离 MATLAB，使用 Eigen 库实现 CBF、MVDR 和 MUSIC 算法。

目标2：使用 Python 训练 DOA（波达方向估计）模型，导出 ONNX，并在 C++ 中使用 ONNX Runtime 进行推理。

技术栈：

  • 开发环境：VS Code + CMake

  • 数学库：Eigen 3

  • AI 推理引擎：Microsoft ONNX Runtime

  • 训练框架：PyTorch

  • 平台：Windows (MinGW) / Linux (WSL)

二、 系统建模

1.1 物理模型

假设有一个均匀线性阵列 (ULA)，包含 $M$ 个阵元，阵元间距 $d$ 为半波长。空间中有一个远场窄带信号 $S(t)$ 从角度 $\theta$ 入射。

接收到的信号模型可以表示为：

$$\mathbf{X}(t) = \mathbf{A}(\theta) \mathbf{S}(t) + \mathbf{N}(t)$$

其中最核心的概念是 导向矢量 (Steering Vector) $\mathbf{a}(\theta)$。它是信号在空间上的“指纹”，描述了信号到达不同阵元时的相位延迟：

$$\mathbf{a}(\theta) = \begin{bmatrix} 1 \\ e^{-j \frac{2\pi}{\lambda} d \sin(\theta)} \\ \vdots \\ e^{-j \frac{2\pi}{\lambda} d (M-1) \sin(\theta)} \end{bmatrix}$$

1.2 MATLAB 数据生成代码

为了验证 C++ 算法的正确性，我们首先在 MATLAB 中生成“标准答案”。

% === MATLAB Data Generation Script ===
clear; clc;

% 参数设置
M = 10;             % 阵元数
K = 100;            % 快拍数 (Snapshots)
theta = 30;         % 目标角度 (度)
SNR = 10;           % 信噪比 (dB)
d_lambda = 0.5;     % 间距/波长

% 1. 构建导向矢量
theta_rad = theta * pi / 180;
index = 0:M-1;
a = exp(-1j * 2 * pi * d_lambda * index' * sin(theta_rad));

% 2. 生成信号 S 和 噪声 N
S = (randn(1, K) + 1j*randn(1, K)) / sqrt(2); % 复高斯信号
N = (randn(M, K) + 1j*randn(M, K)) / sqrt(2); % 复高斯白噪声

% 3. 合成接收数据 X
Ps = 1;             % 信号功率
Pn = Ps / (10^(SNR/10)); % 根据SNR计算噪声功率
X = a * S + sqrt(Pn) * N;

% 4. 导出为 txt (供 C++ 读取)
% 格式：每一行是一个快拍，实部 虚部 交替
fid = fopen('cbf_input.txt', 'w');
for k = 1:K
    for m = 1:M
        fprintf(fid, '%.6f %.6f ', real(X(m,k)), imag(X(m,k)));
    end
    fprintf(fid, '\n');
end
fclose(fid);
disp('数据已生成：cbf_input.txt');

实现三种经典算法：

  1. CBF (常规波束形成)：通过相位补偿进行能量扫描，分辨率较低。

  2. MVDR (最小方差无畸变响应)：利用协方差矩阵的逆 $R^{-1}$ 压制干扰，分辨率较高。

  3. MUSIC (多重信号分类)：基于特征空间分解（子空间方法），利用信号子空间与噪声子空间的正交性，实现超高分辨率测向。

三、 核心算法原理与 C++ 实现

在 C++ 中，我使用 Eigen 3 库来处理复数矩阵运算。相比手动写循环，Eigen 的效率和可读性极高。

2.0 数据预处理：协方差矩阵

所有算法的第一步都是计算 样本协方差矩阵 $R$。它利用时间平均来提取空间相关性，压制随机噪声。

$$\mathbf{R} \approx \frac{1}{K} \sum_{t=1}^{K} \mathbf{X}(t) \mathbf{X}(t)^H = \frac{1}{K} \mathbf{X} \mathbf{X}^H$$

C++ 实现 (Beamformer.cpp):

void Beamformer::computeCovariance() {
    // X_.adjoint() 即共轭转置 (Hermitian)
    R_ = (X_ * X_.adjoint()) / static_cast<double>(n_snapshots_);
}

2.1 常规波束形成 (CBF)

原理：

CBF 是最直观的方法。既然信号从 $\theta$ 来，我就把阵列对准 $\theta$ 听。通过对各个阵元进行相位补偿然后求和，使目标方向信号同相叠加增强。

空间谱公式：

$$P_{CBF}(\theta) = \mathbf{a}(\theta)^H \mathbf{R} \mathbf{a}(\theta)$$

C++ 实现：

double Beamformer::runCBF(double theta_deg) {
    VectorXcd a = getSteeringVector(theta_deg);
    // 纯矩阵运算，直接翻译公式
    MatrixXcd P = a.adjoint() * R_ * a;
    return std::abs(P(0,0));
}

特点：鲁棒性好，但主瓣宽（分辨率低），旁瓣高（抗干扰差）。

2.2 最小方差无畸变响应 (MVDR)

原理：

MVDR (或 Capon) 更加聪明。它试图寻找一个最优权重 $\mathbf{w}$，使得：

  1. 对目标方向 $\theta$ 的增益为 1（无畸变）。

  2. 总输出功率最小（最小方差）。

    这就意味着算法会自动在干扰方向形成“零陷”，从而压制干扰。

空间谱公式：

$$P_{MVDR}(\theta) = \frac{1}{\mathbf{a}(\theta)^H \mathbf{R}^{-1} \mathbf{a}(\theta)}$$

C++ 实现：

double Beamformer::runMVDR(double theta_deg) {
    VectorXcd a = getSteeringVector(theta_deg);
    // 利用 Eigen 求逆
    MatrixXcd R_inv = R_.inverse();
    MatrixXcd denom = a.adjoint() * R_inv * a;
    return 1.0 / std::abs(denom(0,0));
}

特点：分辨率优于 CBF，能有效抑制强干扰，但对阵列误差敏感，且需要求逆（计算量大）。

2.3 多重信号分类 (MUSIC)

原理：

MUSIC 是子空间算法的代表，它不再扫描能量，而是利用数学上的正交性。

对 $R$ 进行特征值分解 (EVD)，可以将空间分为两个正交的子空间：

  1. 信号子空间：大特征值对应的特征向量。

  2. 噪声子空间 ($U_n$)：小特征值对应的特征向量。

理论上，真实信号的导向矢量 $\mathbf{a}(\theta)$ 与噪声子空间 $U_n$ 是正交的，即 $\mathbf{a}^H U_n = 0$。

空间谱公式：

$$P_{MUSIC}(\theta) = \frac{1}{||\mathbf{a}(\theta)^H \mathbf{U}_n||^2}$$

当 $\theta$ 对准目标时，分母趋近于 0，谱峰趋近于无穷大。

C++ 实现：

double Beamformer::runMUSIC(double theta_deg, int n_sources) {
    // 1. 特征值分解 (SelfAdjoint 针对厄米特矩阵优化)
    Eigen::SelfAdjointEigenSolver<MatrixXcd> eig_solver(R_);
    MatrixXcd evecs = eig_solver.eigenvectors();

    // 2. 截取噪声子空间 (前 M-K 列)
    int n_noise = n_sensors_ - n_sources;
    MatrixXcd Un = evecs.block(0, 0, n_sensors_, n_noise);

    // 3. 计算正交性
    VectorXcd a = getSteeringVector(theta_deg);
    MatrixXcd denom = a.adjoint() * Un * Un.adjoint() * a;
    return 1.0 / std::abs(denom(0,0));
}

特点：突破瑞利限，具备超高分辨率。但计算复杂度极高（$O(N^3)$）。

四、 深度学习实现

1. PyTorch 训练

构建了一个简单的多层感知机 (MLP)，输入是协方差矩阵的实部和虚部拼接向量，输出是角度。

# 模型定义 (简单粗暴)
class DOANet(nn.Module):
    def __init__(self, input_dim=200):
        super(DOANet, self).__init__()
        self.net = nn.Sequential(
            nn.Linear(input_dim, 128),
            nn.ReLU(),
            nn.Linear(128, 64),
            nn.ReLU(),
            nn.Linear(64, 1) # 直接回归角度值
        )

2. 导出 ONNX

这是连接 Python 和 C++ 的关键一步。注意： 如果你的 PyTorch 版本较新，导出时可能会遇到 Opset Version 冲突，建议指定 opset_version=18。

torch.onnx.export(model, dummy_input, "doa_model.onnx", opset_version=18, ...)

五、 C++ 部署 ONNX Runtime

需要在 C++ 中加载 .onnx 模型文件，并把 Eigen 的矩阵数据喂给模型。

1. CMake 配置

需要在 CMakeLists.txt 中链接 ONNX Runtime 的库文件。

set(ORT_ROOT "${CMAKE_SOURCE_DIR}/3rdparty/onnxruntime")
include_directories("${ORT_ROOT}/include")
link_directories("${ORT_ROOT}/lib")
target_link_libraries(SmartRadar onnxruntime)

2. C++ 推理代码

核心逻辑：Eigen Matrix -> std::vector -> Ort::Value -> Run -> Result。

double Beamformer::runSmartDOA() {
    // 1. 初始化环境
    Ort::Env env(ORT_LOGGING_LEVEL_WARNING, "SmartRadar");
    Ort::SessionOptions session_options;
    const wchar_t* model_path = L"doa_model.onnx"; // 注意宽字符路径
    Ort::Session session(env, model_path, session_options);

    // 2. 数据预处理：把复数协方差矩阵拉平
    std::vector<float> input_tensor_values;
    for(int i=0; i<R_.size(); ++i) input_tensor_values.push_back((float)R_(i).real());
    for(int i=0; i<R_.size(); ++i) input_tensor_values.push_back((float)R_(i).imag());

    // 3. 执行推理
    auto output_tensors = session.Run(..., &input_tensor, ...);

    // 4. 获取结果
    float* floatarr = output_tensors[0].GetTensorMutableData<float>();
    return (double)floatarr[0];
}

六、 运行结果

在 Windows 和 WSL (Ubuntu) 环境下编译运行，针对一个设置在 30° 的目标，结果如下：

=== Modern C++ Array Processing Framework ===
[Info] Data loaded. Shape:10x100
[Info] Covariance matrix R computed.

CBF Power at 30 degrees: 104.055   (传统算法，波束较宽)
MVDR Power at 30 degrees: 0.931144 (抑制干扰，峰值锐利)
MUSIC Power at 30 degrees: 84.7295 (超高分辨率)

--------------------------------
[AI] Loading ONNX model...
AI Predicted Angle: 30.9711 deg    (端侧推理)

可以看到，AI 模型的预测结果与真实值误差仅在 1° 左右，且无需进行复杂的角度扫描，计算效率提升。

七、 总结

遇到的问题：

1. MATLAB 与 C++ 的数据对齐

  • 问题现象：C++ 读取的数据全是乱码、0，或者算出来的协方差矩阵不对。

  • 原因分析：

    • 内存布局不同：MATLAB 默认是列优先 (Column-major)，而 C++ 通常按行读取。

    • 复数格式：文本文件中实部和虚部的排列方式，C++ 读取时需要严格对应（real >> imag）。

  • 解决方法：

    • 统一标准：在 MATLAB 生成数据时，显式地按“一行一个快拍（Snapshot）”的格式写入 .txt。

    • 代码修正：C++ 读取循环改为外层 k (时间)，内层 m (阵元)。

    • 验证技巧：先读取前 5 行打印出来，肉眼对比 .txt 文件，确保一致。

2. 矩阵运算的精度与接口

  • 问题现象：MVDR 求逆失败，或者 MUSIC 算出来的谱峰位置不对。

  • 关键点：

    • 求逆：必须使用 .inverse()，且矩阵必须满秩（所以仿真时必须加噪声，不能太纯净）。

    • 共轭转置：在数学公式里是 $A^H$。在 Eigen 库里，要用 .adjoint()，千万不能只用 .transpose()（那是转置不共轭）。

    • 特征值排序：MATLAB 的 eig 返回的结果通常是无序或升序，Eigen 的 SelfAdjointEigenSolver 默认是升序。写 MUSIC 取噪声子空间时，要取前 $M-N$ 个（对应小特征值）。

3. PyTorch 版本与 ONNX Opset 冲突

  • 问题现象：运行 train_doa.py 导出模型时报错 RuntimeError: ... Opset version ... 或提示安装 onnxscript。

  • 原因分析：你的 PyTorch 版本较新（2.x），默认支持的高版本 ONNX 算子（Opset 18），而代码里写死了旧版本（Opset 11）。

  • 解决方法：

    • 安装缺失依赖：pip install onnxscript。

    • 修改代码：在 torch.onnx.export 中将 opset_version 改为 18。

4. 输入数据的维度匹配

  • 问题现象：C++ 推理时报错，提示 Input Shape 不匹配。

  • 关键点：

    • PyTorch 训练时：输入通常是 Batch x Dim (例如 64 x 200)。

    • C++ 推理时：通常只推一条数据，所以 Input Shape 必须是 1 x 200。

    • 数据展平：协方差矩阵是二维复数 ($10 \times 10$)，必须按训练时的逻辑（先实部后虚部，或者交替）拉平成一维向量 (std::vector<float>)。

打通了从 理论算法 到 工程落地 的简单流程。

# matlab # c++ # 深度学习 # python
Logo
九章云极普惠算力

更多推荐

cover

阿里云 Coding Plan Lite 下架，各家算力吃紧，上哪买还能支持GLM-5和5.1的coding plan？_2026-04-15

avatar 九章云极普惠算力
cover

看遍了市面上的coding plan，我发现还是这个好用

avatar 九章云极普惠算力

Kandinsky-5.0-I2V-Lite-5s惊艳效果展示：赛博朋克街景→霓虹闪烁+雨滴滑落动态视频

本文介绍了如何在星图GPU平台上自动化部署Kandinsky-5.0-I2V-Lite-5s镜像，实现高效图生视频转换。该工具能将静态赛博朋克街景图片快速转换为动态视频，添加霓虹闪烁、雨滴滑落等效果，适用于短视频制作、广告创意等场景，显著提升内容创作效率。

avatar 九章云极普惠算力