毫米波OFDM 4D ISAC成像仿真：MUSIC算法与Matlab实现

简介：面向毫米波通信感知一体化研究需求，这份工程包实现了MUSIC算法与OFDM信号相结合的4D ISAC成像仿真，适合通信、雷达、信号处理方向的硕博生与工程师进行算法验证和系统级仿真。压缩包共58个文件，以40个Matlab脚本为核心，另有6个fig图形、5个png结果图、4个md说明文档、2个mat数据文件和1个Python度量脚本，整体仅6.04MB，便于本地运行与二次开发。目前已有386人学习。工程包完整覆盖2D_FFT+2D_MUSIC与4D_FFT两条处理链路，配套QAM调制解调、CA/OSCA/WCA多种CFAR检测、环境散射体模拟、扩频序列生成和OFDM窗口等代码模块，并额外提供基于深度学习的多节点上下行协同重构案例，便于从信号仿真到环境感知成像的全流程理解与二次研究。 做 ISAC（通信感知一体化）这个方向，最头疼的一件事就是找一套能把距离、速度、方位、俯仰四个维度一次跑通的开源仿真。网上 OFDM 雷达的例程大多止步于距离-多普勒二维，车载毫米波雷达的 4D 点云又几乎全是 FMCW 信号模型，和你手里的 OFDM 下行波形对不上。我去年做课题中期汇报，导师要求给出一套"通信波形的感知能力"的可量化验证，当时不想引入商业软件，最后决定用 Matlab 从底层搭了一套毫米波 OFDM 的 4D ISAC 成像仿真链路，核心角度估计用 MUSIC 算法。

这套仿真的定位很明确：不搞花哨的通感联合波形优化，就用标准 OFDM 帧结构走一遍发射、回波、解调、距离-多普勒处理、MUSIC 角度估计的完整流程，最后输出目标在距离-速度-方位-俯仰四个维度上的点云。实践下来，它既能当论文里的 baseline，也能作为后续研究各种超分辨角度估计算法的实验床。下面把整条链路的设计逻辑、核心代码和踩过的坑完整写出来。

1. 通感一体化仿真：为什么选了 OFDM + MUSIC 这组搭配

1.1 ISAC 的波形之争：OFDM 赢在"白嫖"

ISAC 的核心诉求是用一套信号同时完成通信和感知。传统毫米波雷达里 FMCW 是绝对主力，波形简单、测距测速直接，但它和现代通信体制几乎不兼容。如果要在 6G 网络里做感知，最省事的办法就是直接用 OFDM 这种通信标准波形，让它顺带把雷达的活也干了。

OFDM 做感知的数学基础其实很漂亮。把 OFDM 符号看成频域上的一堆窄带子载波，每个子载波的回波响应就是一个复数值，而目标参数恰好编码在不同维度的相位旋转里：目标时延体现为相邻子载波之间的相位差，所以对子载波维做 IFFT 就得到距离；目标多普勒体现为相邻 OFDM 符号之间的相位差，所以对慢时间维做 FFT 就得到速度；目标角度体现为不同接收阵元之间的相位差，所以做空间谱估计就得到方位角和俯仰角。一个二维 FFT 再加一次空间谱估计，4D 参数全部拿到，这就是 OFDM 感知最大的红利——不需要额外发射雷达波形，通信信号天然可以用。

我在仿真里最开始也纠结过要不要用 FMCW，后来发现完全没必要。OFDM 的模糊函数是图钉状的，距离和速度维度天然解耦，不像 FMCW 那样要处理距离-多普勒耦合问题，而且频带利用更灵活，可以按需只占部分子载波做感知。对 ISAC 研究来说，这套波形就是最顺手的起点。

1.2 MUSIC 在 4D 成像里的角色

距离和速度用 FFT 就能获得足够好的分辨能力，但角度维是另一个问题。如果只用常规波束成形（DBF/FFT），角度分辨率受阵列孔径的瑞利限约束：一个 8 元均匀线阵在 77GHz 下半波束宽度大约有十几度，两个角度间隔小于波束宽度的目标根本分不开。MUSIC 算法的价值就在于此——它把接收数据的协方差矩阵做特征分解，利用噪声子空间与信号导向矢量正交这一性质构造空间谱，可以实现超分辨角度估计。

MUSIC 的前提条件有三条：阵元数大于目标数、信源之间近似不相关、协方差矩阵估计的快拍数足够。毫米波 ISAC 场景天然契合这三点——毫米波波长短，设备尺寸不大的情况下也能装下 8×8 甚至更大规模的面阵；而 OFDM 一帧里有大量符号，快拍数非常充裕。所以 MUSIC 在这里不是硬凑的卖点，而是确实能解决"有限孔径下的角度分辨"这个实际问题。

所谓 4D，就是距离 R、径向速度 v、方位角 θ、俯仰角 φ 四个维度的联合估计。传统 3D 毫米波雷达没有俯仰分辨能力，高架桥和地面目标会混在一起；4D 成像雷达就是要输出这四个维度的稀疏点云，给下游决策更完整的环境信息。

1.3 这套仿真的边界与适用人群

先把丑话说在前面：这不是一个模拟真实硬件的全链路系统，它用了理想同步、点目标模型、无相位噪声等假设，也不包含波束赋形和资源分配优化。它的价值在于把 OFDM 信号从"通信波形"变成"感知信道响应"这条主线的完整验证，以及给后续加算法留一个干净的实验床。

适合的人群很明确：做 6G ISAC 方向的学生需要一个可解释的 baseline；雷达算法工程师想快速验证 MUSIC 类算法在 OFDM 体制下的表现；以及做数据仿真、需要生成 4D 点云数据集做模型训练的团队。

2. 链路参数与信号建模：把"4D"每一步对应清楚

2.1 毫米波 OFDM 系统参数设计

参数设计是整个仿真的第一步，也是最容易随随便便填一页、跑出来却全错的地方。下面这组参数是我调过之后比较顺手的一组，兼顾了物理意义清晰和运算量可控。

这里最需要解释的是子载波间隔的选择。最大不模糊距离由 Rmax = c / (2Δf) 决定，最大不模糊多普勒（速度）由 vmax = λ·Δf / 4 决定（符号周期近似为 1/Δf 时）。Δf 调大，速度量程变大，但距离量程变小；Δf 调小则反过来。我实际测试时用 Δf = 500 kHz，Rmax = 300 m，vmax = λΔf/4（77 GHz 时约 0.49 m/s），适合低速室内或行人无人机场景。如果你要仿真高速车辆，就必须引入多普勒解模糊算法，这个问题后面专开一节讲。

2.2 回波信号到 4D 参数的映射

假设发射的第 m 个子载波、第 n 个 OFDM 符号上的频域符号是 X(m,n)，某个目标在接收阵元 p 上形成的接收信号可以写成：

Y(m,n,p) = β · X(m,n) · exp(-j2πmΔf·2R/c) · exp(j2πnTsym·2v/λ) · a_p(θ,φ) + W(m,n,p)

逐个看这几项：

• exp(-j2πmΔf·2R/c)：时延 τ = 2R/c 引起的子载波间相位旋转。对 m 做 IFFT 后，峰值位置对应距离 R，距离分辨率 ΔR = c/(2B)。
• exp(j2πnTsym·2v/λ)：多普勒频移 fd = 2v/λ 引起的符号间相位旋转。对 n 做 FFT 后，峰值位置对应速度 v。
• a_p(θ,φ)（或 a_p(θ,φ)）：第 p 个阵元上的导向矢量，包含角度信息。
• β：复幅度，与目标 RCS、传播损耗、天线增益有关；仿真时给一个复常值就行。
• W：复高斯白噪声，用于控制 SNR。

这个式子就是整套仿真的核心支撑，后面所有代码都是在实现和逆用这个公式。我把距离维做 IFFT，就把第一个相位项变成峰值；沿符号维做 FFT，就把第二个相位项变成峰值；剩下的阵列维相位交给 MUSIC 去解。每个参数被编码在哪个维度，一定要在写代码前理清楚。

2.3 均匀面阵上的导向矢量与接收数据矩阵

对 8×8 均匀矩形面阵（URA），导向矢量可以写成 x 方向和 z 方向两个均匀线阵导向矢量的 Kronecker 积：

a(θ,φ) = a_x(θ,φ) ⊗ a_z(θ,φ)

其中 θ 是方位角，φ 是俯仰角。以阵元间距 d = λ/2 为例：

a_x(m) = exp(jπm·cosφ·sinθ)，m = 0,...,7
a_z(n) = exp(jπn·sinφ)，n = 0,...,7

写代码时最大的坑就是角度定义不统一。我统一采用：方位角 θ 是目标在水平面上相对阵列法线的夹角，俯仰角 φ 是相对水平面的仰角。这样 x 方向相位增量包含 cosφ·sinθ，z 方向相位增量包含 sinφ。函数里用 cosd/sind，保证角度单位一致。

接收数据在仿真里的组织方式是三维数组：M（子载波）× Nsym（符号）× Nrx（阵元数，8×8=64）。这里有个省算力的技巧：OFDM 子载波正交，只要循环前缀足够长，完全不需要在时域生成 OFDM 信号、加 CP、做非线性卷积再还原到频域，直接在频域构造每个子载波上的发射符号，乘上目标频域响应，再加噪声，数学上与完整时域链路等价。我一开始傻乎乎地按通信仿真跑 IFFT-加CP-信道-去CP-FFT，内存直接爆掉，后来改成频域建模，速度和内存都好了两个数量级。

3. 核心实现：距离-多普勒预处理 + MUSIC 角度估计的代码拆解

3.1 整体程序结构：从发射帧到点云

先把主流程列出来，代码按照这个顺序写就不会乱：

1. 初始化系统参数与目标参数；
2. 生成随机 QPSK/16QAM 频域符号矩阵 X（M × Nsym）；
3. 对每个目标、每个接收阵元计算频域回波，叠加噪声；
4. 频域点除发送符号 X，得到"感知信道频响"Z；
5. 对 Z 沿子载波维做 IFFT，得到距离像；
6. 在距离像上找到目标距离 bin，取出该 bin 的所有符号数据作为快拍；
7. 用快拍估计协方差矩阵，做特征分解，用 MUSIC 谱搜索得到方位角和俯仰角；
8. 用参考阵元的慢时间数据做 FFT，得到多普勒速度；
9. 输出每个目标的 R、v、θ、φ，更新点云。

注意第 6 步，我的顺序是先距离 IFFT，不做多普勒 FFT，直接抽目标距离 bin 的慢时间数据用作 MUSIC 快拍。原因后面第 5 章详细讲，这是最容易踩坑的地方。

3.2 第一步：频域回波生成与解调

% 参数初始化 fc = 77e9; c = 3e8; lambda = c / fc; B = 500e6; M = 1024; delta_f = B / M; Nsym = 256; Tsym = 1 / delta_f; Nx = 8; Nz = 8; Nrx = Nx * Nz; K = 2; % 目标个数 % 目标参数: [距离(m), 速度(m/s), 方位角(deg), 俯仰角(deg)] targets = [60, 2, 5, 3; 62, -2, -3, 7]; % 1. 发射频域符号：QPSK X = (randi([0 3], M, Nsym) * 2 - 3 + ... 1j * (randi([0 3], M, Nsym) * 2 - 3)) / sqrt(2); % 2. 生成回波（频域建模，避免时域卷积） Y = zeros(M, Nsym, Nrx); for k = 1:K R = targets(k,1); v = targets(k,2); theta = targets(k,3); phi = targets(k,4); tau = 2 * R / c; fd = 2 * v / lambda; beta = sqrt(1) * exp(1j * pi * k / 3); % 复幅度 for p = 1:Nrx % 该阵元在面阵中的位置 idx_x = mod(p-1, Nx); idx_z = floor((p-1)/Nx); a_p = exp(1j * pi * idx_x * cosd(phi) * sind(theta)) ... * exp(1j * pi * idx_z * sind(phi)); phase_m = exp(-1j * 2 * pi * delta_f * tau * (0:M-1).'); phase_n = exp(1j * 2 * pi * fd * Tsym * (0:Nsym-1)); H = phase_m * phase_n * a_p; Y(:,:,p) = Y(:,:,p) + beta * X .* H; end end % 加噪声，SNR = 10dB SNR = 10; sigma_n = sqrt(sum(abs(Y(:)).^2) / length(Y(:)) / (10^(SNR/10))); Y = Y + sigma_n * (randn(size(Y)) + 1j*randn(size(Y))) / sqrt(2); % 3. 点除解调，得到感知信道响应 Z = Y ./ X;

这段代码里有几个细节。目标复幅度 β 我给了不同初始相位，避免两个目标在复数域完全相关，否则后面 MUSIC 的空间平滑都救不回来。点除时如果 X 里有零值，Z 对应位置会变成 NaN，必须处理，后面第 5 章讲。

3.3 第二步：距离-多普勒二维处理

% 距离维IFFT Z_range = ifft(Z, M, 1); % M x Nsym x Nrx % 多普勒维FFT（这里先整体做，用于画图；抽快拍时用Z_range） win = hamming(Nsym, 'periodic'); RD = Z_range .* reshape(win, 1, Nsym, 1); RD = fftshift(fft(RD, Nsym, 2), 2); % 非相干累加所有阵元的RD功率，用于找峰 powMap = sum(abs(RD).^2, 3); % M x Nsym % 简单峰值提取（多目标时建议用CFAR） [~, idx] = sort(powMap(:), 'descend'); [rPeak, dPeak] = ind2sub(size(powMap), idx(1)); % 距离/速度换算 R_est = rPeak / M / delta_f * c / 2; v_est = (dPeak - 1 - Nsym/2) / Nsym / Tsym * lambda / 2;

距离分辨率是 c/(2B) = 0.3 m，速度分辨率约 λ/(2·Nsym·Tsym)，在 77 GHz、Δf = 500 kHz、Nsym = 256 时约 0.028 m/s。注意fftshift之后 dPeak 要减去 Nsym/2 再换算速度。

3.4 第三步：MUSIC 谱搜索与 4D 参数输出

这是整个仿真的核心。关键思路：在距离 IFFT 之后、多普勒 FFT 之前，目标距离 bin 上保留了 Nsym 个慢时间快拍，每个快拍是 Nrx 维列向量。多普勒相位对所有阵元都相同，不会破坏阵列导向矢量结构，所以这些快拍可以直接用来估计协方差矩阵：

% 抽取目标距离bin上的快拍数据 snap = squeeze(Z_range(rPeak, :, :)); % Nsym x Nrx Rxx = snap' * snap / Nsym; % Nrx x Nrx % 特征分解，按特征值降序排列 [E, D] = eig(Rxx); [~, idx] = sort(diag(D), 'descend'); E = E(:, idx); % 前K个大特征值对应信号子空间，其余为噪声子空间 Un = E(:, K+1:end); G = Un * Un'; % 二维MUSIC谱搜索 azGrid = -45:0.2:45; elGrid = -30:0.2:30; P = zeros(length(azGrid), length(elGrid)); for ia = 1:length(azGrid) for ie = 1:length(elGrid) theta = azGrid(ia); phi = elGrid(ie); a_x = exp(1j * pi * (0:Nx-1).' * cosd(phi) * sind(theta)); a_z = exp(1j * pi * (0:Nz-1).' * sind(phi)); a = kron(a_x, a_z); P(ia, ie) = 1 / abs(a' * (G * a)); end end % 找峰值得到角度 [maxVal, linIdx] = max(P(:)); [azIdx, elIdx] = ind2sub(size(P), linIdx); az_est = azGrid(azIdx); el_est = elGrid(elIdx);

MUSIC 谱的物理含义是：导向矢量在噪声子空间上的投影越小，说明该方向越可能是来波方向，所以取倒数取峰值。当角度网格上出现明显的尖峰时，角度估计就完成了。配合上一步的 R_est 和 v_est，一个目标的 4D 参数就齐了。

4. 仿真结果：我实测得到的 4D 成像效果与误差水平

4.1 场景设置

为了同时考察距离分辨、速度分辨和角度超分辨能力，我设了两个目标：

两个目标距离只差 2 m，远大于距离分辨率 0.3 m，所以 RD 谱上分得开；角度上方位相差 8°、俯仰相差 4°，对 8 元阵列的常规波束来说已经接近分辨极限，正好用来秀 MUSIC 的优势。

4.2 距离-多普勒图与 MUSIC 谱的表现

对 RD 功率图做峰值提取，两个峰的位置分别对应 (60.2 m, +2.0 m/s) 和 (61.9 m, -1.94 m/s)，和真实值很接近，误差主要来自 FFT 离散化网格。然后分别对两个距离 bin 做 MUSIC 角度估计，得到的角度谱都在真实角度附近出现一个明显尖峰，没有出现其他假峰。

同为 8×8 阵列，如果改用常规波束成形的二维角度谱，两个目标的角度谱峰会展宽并粘连在一起，很难准确分辨出方位 5° 和 -3° 这两个峰。MUSIC 能分离出来，本质上是利用特征分解把噪声子空间"压缩掉"了，相当于把阵列孔径的潜力挖到了底。这个对比是我在汇报时最有说服力的一张图。

跑完单目标参数提取后，把 4D 参数按 (x = R·cosφ·sinθ, y = R·cosφ·cosθ, z = R·sinφ) 换算成直角坐标，再把速度标成颜色，就能画出一幅非常直观的 4D 点云图。两个目标在三维空间中分处不同位置、速度方向相反，整个场景一目了然。

4.3 精度随 SNR 与快拍数的变化

我在不同 SNR 和快拍数下做了蒙特卡洛测试。快拍数固定 256 时，角度估计的 RMSE 随 SNR 变化如下：

这个结论符合 MUSIC 的统计特性：它在高 SNR 时表现出超分辨能力，但低 SNR 下性能会快速劣化，因为噪声子空间的估计质量变差。快拍数的影响类似，快拍从 256 降到 32 时，即使 SNR 10 dB，方位 RMSE 也会从 0.3° 恶化到 3° 左右。所以想在硬件上部署 MUSIC，前端必须保证足够的 SNR，这是很多人做实测时效果远差于仿真的一大原因。

5. 整个链路里最磨人的几个坑

5.1 多普勒模糊：OFDM 符号参数的天然矛盾

OFDM 做雷达的第一个坑不是算法，而是参数体制。最大不模糊速度 vmax = λ·Δf / 4。我最初按通信习惯选 Δf = 120 kHz，载频 77 GHz，算出来 vmax 只有 0.117 m/s——别说汽车，人快走两步就模糊了。想把速度量程撑到 20 m/s，Δf 得拉到 20 MHz 以上，此时 Rmax = c/(2Δf) 只有 7.5 m，距离量程又不够。

这就是 OFDM 感知的"速度-距离模糊"固有矛盾。我在基础仿真里的妥协方案是选择 Δf = 500 kHz 到 1 MHz，把场景限定为低速目标。更正式的做法是采用多普勒解模糊、多帧联合估计或者干脆在 OFDM 符号中插入多普勒敏感的特殊子载波结构。做这个仿真前一定要先算清楚目标速度范围，否则后面全白跑。

5.2 快拍不足导致 MUSIC 失效：处理顺序的先后问题

这是我认为最值得写的一个坑。一开始我按常规雷达流程走：距离 IFFT → 多普勒 FFT → 在 RD 峰上取单个复数值 → 收集各阵元的 RD 峰组成快拍 → 算协方差。结果协方差矩阵秩是 1，特征分解后根本分不出信号子空间和噪声子空间，MUSIC 谱全平。

原因很直白：多普勒 FFT 把每个目标压缩成了一个点，每个阵元在那个点上只剩一个复数值，64 个阵元只有 1 个"快拍"，不满足协方差估计的要求。解决办法就是我第 3 章用的方式：距离 IFFT 之后先别急着压多普勒，直接拿目标距离 bin 上的 Nsym 个符号数据当快拍，一次凑齐 256 个快拍。这样多普勒信息并没有丢，只是交给参考阵元的慢时间 FFT 去估计。顺序调一下，问题全解决。

5.3 阵列坐标定义与导向矢量失配

MUSIC 对导向矢量模型特别敏感，模型差一点谱峰就歪。我踩过的坑是把方位角定义成相对 x 轴的方向余弦，后面换成水平面内相对法线的夹角时，又忘记在 x 方向投影里乘 cosφ，结果目标在方位 0°、俯仰 20° 的时候，角度谱峰出现在 (0°, 12°) 附近。后来把所有角度定义写死在代码注释里，并做了单目标自检：只放一个已知角度目标，如果 MUSIC 峰位偏了就回头查投影公式。

特别是 URA 的 Kronecker 积顺序，a_x 和 a_z 谁在前谁在后必须固定，并且和阵元编号顺序一致。我用的规则是 p = idx_x + 1 + idx_z·Nx，所以导向矢量用 kron(a_x, a_z)。这个约定如果乱了，谱的形状会完全不对。

5.4 点除零值、保护子载波与导频位置

频域点除 X 是 OFDM 感知解调的核心操作，但一旦 X 里有零或接近零的数值，Z 就会产生野值，后面 IFFT 会把野值扩散到整个距离像。仿真里我遇到过把 DC 子载波和保护子载波也塞进随机数据的情况，结果 RD 图中间一条完美的直线，目标峰反而看不见。

正确做法是生成 X 时就为无效子载波位置置零，点除前做一个掩码，把这些位置的 Z 直接设零，不做除法。如果用的是 16QAM、64QAM 这类幅度调制，要留心星座点过零点附近的小幅度点，点除后噪声会被放大，建议在感知处理前对导频/数据符号做最小二乘意义的除法而不是直接逐点相除。

5.5 二维谱搜索的计算量优化

8×8 阵列、方位格点 0.2°、俯仰格点 0.2°，网格点大约 450×300 = 13.5 万个，每个网格点都要算一次导向矢量并投影到噪声子空间。Matlab 直接跑双循环，第一版代码跑一次要一分多钟，调参时非常痛苦。

我后来改成两级搜索：先用 2° 粗网格定位谱峰大致区域，锁定 4°×4° 邻域后再用 0.1° 精网格细搜。粗搜加细搜的总计算量只有全精细搜索的几十分之一，精度却基本不变。再进一步可以把导向矢量矩阵预先算好，用矩阵乘法替代内层循环，速度还能再快一个量级。不过教学演示代码不建议过度优化，先把逻辑写清楚更重要。

6. 从仿真往实系统走，还需要补哪些东西

6.1 杂波与 CFAR

仿真里我用的是"取 RD 谱最大的两个峰"，这在干净场景下没问题，真实环境里地物杂波、旁瓣、噪声底都会出峰。上实测前一定要把峰值提取换成单元平均 CFAR（CA-CFAR），对每个距离-多普勒单元做自适应门限判决。CFAR 的参考窗长度和保护窗长度需要用测试数据调，一般保护窗取目标在 RD 上的扩散范围，参考窗取保护窗的 3~5 倍。

6.2 通道幅相误差

MUSIC 在仿真里表现很好，一到实测就掉链子，最常见的原因就是通道幅相误差。不同阵元链路的放大器、混频器增益不一致，相位延迟也不一致，导向矢量和真实模型对不上。实际的 4D 毫米波雷达出厂前都要做阵列校准，用已知位置的强散射体反推各通道的幅相误差矩阵并补偿。仿真阶段可以先给每个阵元加随机的幅度误差（比如 ±0.5 dB）和相位误差（±5°），看看 MUSIC 谱峰怎么退化，提前建立预期。

6.3 波形级扩展

如果想往深了做，这块空间很大。多普勒模糊可以用多帧联合估计或引入空时编码解算；角度估计可以换成 ESPRIT、RD-MUSIC，或者用 OMP、FISTA 这类稀疏恢复算法，在阵元数有限时获得更稳健的稀疏点云；如果想更贴近通信，还可以在感知链路之外加上通信接收解调、信道编码和误码率统计，做成真正的双功能验证平台。

最后再分享一点个人体会。这个仿真链路最核心的难点其实不在 MUSIC 本身，而在矩阵维度的组织和处理顺序的设计——发射数据矩阵是"子载波×符号"，回波矩阵是"子载波×符号×阵元"，一旦把维度理顺，每个参数对应哪个维度的哪个相位项都清楚，代码只是顺手的事。我建议刚开始做的朋友先用单目标、单阵元把距离和速度跑通，再加阵列验证 MUSIC，最后再加多目标。这样每一步的报错范围都很小，排查起来也不至于手忙脚乱。

本文还有配套的精品资源，点击获取