延时求和波束形成:从麦克风阵列到毫米波雷达

延时求和波束形成:从麦克风阵列到毫米波雷达 简介延时求和波束形成是雷达与麦克风阵列信号处理中的基础技术通过调整各阵元接收信号的时间延迟可在期望方向形成主瓣并抑制其他方向干扰。该MATLAB代码资源面向信号处理初学者和雷达、声学阵列方向的开发者提供最简可运行的实现示例适合快速理解延迟计算、加权求和与波束图形成等核心流程。压缩包共2个文件以.m源码为主并含原始打包存档整体仅2KB轻量便于快速仿真与修改。目前已有504人学习下载适合作为波束形成课程实验或入门练习的辅助材料。借助该源码读者可清晰看到从阵元延迟、加权合成到波束输出的完整代码路径并能在此基础上进一步开展目标定向、干扰抑制、信噪比提升以及多目标扫描等扩展实验同时算法也适用于麦克风阵列的声学场景可帮助消除背景噪声、提升语音采集质量。正如描述所指延时求和波束形成能够增强特定方向回波、提升系统信噪比该源码为理解这些实际效果提供了直观的练习入口。1. 延时求和波束从麦克风阵列到毫米波雷达的共性信号模型拿到一份多通道录音比如 ES8311 双麦克风评估板的 I2S 输出直接做延时求和做完发现方向性还不如单麦——这是很多初涉波束形成的工程师都撞过的墙。原因不在求和而在时延估计延时求和波束形成Delay-and-Sum Beamforming本质上是用一组已知的几何延迟把阵列各通道对齐到同一等相位面再叠加获得阵列增益。它同时是雷达阵列如 TI AWR2243 的虚拟阵列和麦克风阵列里最基础、最容易落地的波束策略。本文从数学推导到 MATLAB 仿真覆盖角度扫描、分数时延实现、声学前端硬件参数以及和雷达距离方程、4D 毫米波雷达目标检测的衔接点适合需要快速验证“这个角度能不能分开两个目标”的算法工程师和硬件工程师。2. 线性阵列的延迟模型与 time_delay_beamforming.m 实现2.1 均匀线性阵列的几何时延公式延时求和波束的输入是一组空域采样信号。以 N 元均匀线阵为例阵元间距为 d信号来自远场入射方向与阵列法线夹角为 θ。此时相邻阵元接收到的信号存在一个固定的传播时间差τ d * sin(θ) / c其中 c 是传播速度。射频雷达场景下 c 取光速 3e8 m/s声学场景下 c 取声速约 343 m/s。这个公式只对远场成立目标距离 R 必须远大于阵列孔径 D (N-1)*d一般要求 R 2D²/λλ 是信号波长。近场场景例如桌面麦克风阵列拾取 20cm 外的说话人走的是球面波模型每个阵元的时延要按各自到声源的距离单独算不能直接用 sin(θ) 近似这点在后续麦克风阵列章节会再讨论。时延计算的精度直接决定波束指向的准确性。对窄带信号而言时延 τ 等价于相位差 2πf*τ所以窄带情况下延迟求和可以退化为相位加权求和也就是相控阵里常用的“相移波束形成”。但对宽带信号语音、线性调频雷达回波相位近似只在窄带假设下成立必须真正做时间轴上的延迟对齐。这也是 delay-and-sum 名称里“delay”存在的意义——它不是在频域乘一个复数权重而是在时域把信号搬移。MATLAB 仿真中建议先按实信号做时域延迟验证通过后再换成复基带信号做相移等效两个结果在窄带条件下应当一致。2.2 time_delay_beamforming.m 核心逻辑拆解源码包里的time_delay_beamforming.m实现了最基本的延迟对齐-求和过程。我一般会把它拆成两个函数一个负责计算各通道的整数采样点延迟一个负责叠加输出。下面是最小可运行版本逻辑与源码包保持一致function [beamformed, delay_samples] delay_sum_beamforming(x, fs, d_mic, theta_deg) % x: [N_channel, N_sample] 多通道输入信号, N_channel为阵元数 % fs: 采样率(Hz) % d_mic: 阵元间距(m), 均匀线阵 % theta_deg: 期望波束指向角度(deg) N_channel size(x, 1); N_sample size(x, 2); c 343; % 声速, 雷达场景改成3e8 theta theta_deg * pi / 180; % 计算每个通道相对参考阵元的时延(秒) % 参考阵元取第一个, 其余通道按几何关系计算 delay_s (0:N_channel-1). * d_mic * sin(theta) / c; % 转换为采样点, 四舍五入取整 % 注: 这里做了round, 只保留整数采样点延迟 delay_samples round(delay_s * fs); % 对齐并叠加 beamformed zeros(1, N_sample); for ch 1:N_channel shift delay_samples(ch); if shift 0 % 正延迟: 该通道需要往后平移 aligned [zeros(1, shift), x(ch, 1:end-shift)]; else % 负延迟: 该通道需要往前截断 aligned [x(ch, -shift1:end), zeros(1, -shift)]; end beamformed beamformed aligned; end beamformed beamformed / N_channel; % 归一化, 保持信号幅度一致 end代码里有三个关键参数采样率 fs 决定时延量化精度阵元间距 d_mic 决定阵列有效频率范围指向角度 theta_deg 决定波束扫描方向。需要注意的是round操作会引入最大半个采样周期的量化误差——如果 fs 16kHz误差约 31μs对应语音信号高频段3kHz 以上会产生明显的相位偏差。所以工程实现里通常不会直接使用整数延迟而是配合分数时延滤波器这个在最后一章单独展开。调用逻辑也很直接先把所有通道读入矩阵 x然后用一个循环对角度扫描比如从 -60° 扫到 60° 步进 1°每个角度调用一次上面的函数输出能量的峰值对应的角度就是目标方向。这是最原始但最稳的 DOA 估计方式计算量随角度扫描密度线性增长适合阵元数少4 元以内的场合。2.3 波束图验证与栅瓣抑制条件写完延时求和函数后第一件事是画波束图不是直接丢实测数据。窄带情况下延时求和等效于对各通道施加相位旋转阵列响应可以写成B(theta) sum_{n1..N} exp(j * 2*pi*f * (n-1) * d * sin(theta) / c)把期望指向角度和目标角度代入得到方向增益。用 MATLAB 画出来能直观看到主瓣宽度、旁瓣水平和栅瓣位置。下面的脚本生成波束图f 8e3; % 信号频率 8kHz fs 48e3; % 采样率 48kHz d c / f / 2; % 半波长间距, 避免栅瓣 N 8; % 8 元阵 theta_scan -90:0.5:90; theta_target 30; % 期望指向 30 度 % 计算每个扫描角度对应的阵列响应 B zeros(size(theta_scan)); for idx 1:length(theta_scan) % 各通道相对参考通道的相位差 phase_diff 2 * pi * f * (0:N-1). * d * (sind(theta_scan(idx)) - sind(theta_target)) / c; B(idx) abs(sum(exp(1j * phase_diff))); end B B / max(B); % 归一化到 0dB figure; plot(theta_scan, 20*log10(B eps)); xlabel(扫描角度 (deg)); ylabel(归一化增益 (dB)); grid on; ylim([-40 5]);运行后会看到主瓣在 30° 附近第一旁瓣约 -13dB。这是均匀加权延时求和的固有特征旁瓣高度与阵元数无关只与窗函数有关。如果你想压低旁瓣把简单求和改成加窗求和比如 Hamming 窗或 Taylor 窗代价是主瓣变宽、角度分辨率下降。一个重要边界条件阵元间距 d 不能超过半个波长。如果 d λ/2扫描到非期望角度会出现栅瓣——此时目标方向的信号会被多个角度同时响应DOA 估计会出歧义。雷达场景尤其要小心毫米波雷达的阵列设计通常直接按 λ/2 固定布阵而麦克风阵列为兼顾低频段阵元间距可比半波长更大但要确认工作频段不触发栅瓣。3. 麦克风阵列实战声速、采样率与 ES8311 双麦前端3.1 声学与射频时延的尺度差异麦克风阵列和雷达阵列用的是同一套延时求和理论但工程参数差异巨大。第一个差异是传播速度声速 343 m/s 只有光速的百万分之一所以声学场景里同样 4cm 的阵元间距时延约 116μs雷达场景下 4cm 间距对应 133ps——差了六个数量级。这决定了实现路径完全不同麦克风阵列可以在时域直接做延迟线雷达阵列尤其是窄带雷达必须转换成相位加权因为子采样级别的时延在射频硬件上没法直接实现。第二个差异是信号带宽。语音信号带宽约 4kHz中心频率与带宽之比相对带宽接近 1是典型的宽带信号不能简单用一个相位旋转替代时延雷达信号通常是窄带线性调频或单脉冲相对带宽通常小于 10%所以相移波束形成在雷达里是主流。理解了这一点就能明白为什么同一个 “delay-and-sum” 名字在声学文章里看到的是抽头延迟线在雷达文章里看到的是复数乘积相加。3.2 到达时间差估计与语音对齐麦克风阵列做延时求和时最常见的坑是采样率跟不上时延精度。ES8311 这类音频编解码器最常用的采样率是 48kHz 或 16kHz48kHz 下每个采样周期约 20.8μs对应声程差约 7.1mm。如果两个麦克风间距 4cm那么正面来声入射角 90°延迟约 5.6 个采样点30° 入射则是 2.8 个采样点——后者直接四舍五入会产生很大的指向误差。角度不同对应的延迟采样点如下表目标角度°时延μs48kHz 延迟采样点16kHz 延迟采样点00001534.41.65四舍五入 20.55四舍五入 13058.02.78四舍五入 30.93四舍五入 14581.83.93四舍五入 41.31四舍五入 1从上表可见16kHz 采样率下 15° 和 30° 的延迟都被量化成 1 个采样点这意味着波束指向在这两个角度之间几乎没有区分能力。所以要做高精度指向要么提高采样率代价是 I2S 带宽和功耗要么在时延对齐时使用分数延迟插值。实际工程上我一般先用互相关估计到达时间差TDOA再用估计结果做精细对齐而不是纯靠几何计算。下面给一段常用的 TDOA 估计代码这段代码在调试双麦克风阵列时非常实用[y1, fs] audioread(mic1.wav); [y2, ~] audioread(mic2.wav); % 截取相同长度 len min(length(y1), length(y2)); y1 y1(1:len); y2 y2(1:len); % 互相关估计延迟 [corr, lags] xcorr(y1, y2, coeff); [~, idx] max(abs(corr)); lag lags(idx); % 正数表示信号2滞后于信号1 fprintf(估计时延: %.2f ms (%.2f 采样点)\n, lag/fs*1000, lag); % 可视化 figure; plot(lags/fs*1000, corr); xlabel(延迟 (ms)); ylabel(互相关系数); grid on;这里xcorr返回的峰值位置对应两个通道的采样点延迟。如果峰值不明显先检查信号是否包含足够的宽带成分——纯单音信号互相关会出现多个等高峰值无法判断唯一延迟。语音信号通常是宽带信号峰值一般很干净。另外要注意互相关的延迟分辨率只有 1/fs和直接几何计算的精度相当想要亚采样精度需要用抛物线插值对峰值附近做拟合或者用频域相位法。3.3 双麦克风阵列与 ES8311 的硬件连接细节ES8311 是 Cirrus Logic 的 24bit 低功耗音频 codec常见的双麦克风评估板上两颗 ES8311 各自通过 I2S 输出一路麦克风数据作为 I2S 输入。延时求和算法对两路数据的时间同步要求非常高如果两条 I2S 总线之间存在几百微秒的启动延迟差几何时延计算就没有意义了——所以硬件设计上两颗 codec 必须共用一个 MCLK 和 LRCLK从根上避免采样时钟偏差。软件上采集时也要先做同步校准放一个固定位置的白噪声源用互相关测出两通道固定延迟这个延迟在后续算法中作为一个常量补偿。麦克风阵列的阵元间距不是随便定的。语音信号上限频率约 8kHz考虑女性声和 s 音半波长约 2.1cm所以正常语音增强的双麦间距取 1.5~2.5cm 比较常见如果间距过大高频段会出现空间混叠也就是栅瓣。反过来间距太小时延差太小低频段相位差辨识度下降。做会议麦克风远场拾音时间距设 4cm 左右兼顾中频段的波束宽度和硬件体积。这些参数在time_delay_beamforming.m里对应d_mic变量换硬件时只改这一处和fs就能跑通算法流程。4. 雷达方向距离方程、4D 毫米波雷达与 AWR2243 数据衔接4.1 雷达距离方程中的阵列增益雷达系统用延时求和波束形成收益直接体现在雷达距离方程里。理想点目标回波信噪比可以写成SNR Pt * Gt * Gr * λ² * σ / ((4π)³ * R⁴ * kTBF)其中 Gt 和 Gr 是发射和接收增益。在相控阵或 MIMO 雷达中接收阵列通过延时求和获得 N 元阵列增益G 正比于阵元数 N。更直白的说法把 N 个通道相干叠加信号幅度增加 N 倍噪声功率增加 N 倍非相干噪声信噪比增加 N 倍换算成 dB 就是 10log10(N)。4 元阵带来约 6dB 增益16 元阵带来 12dB探测距离按四次方根关系对应提升约 1.4 倍和 1.8 倍。前提是“相干叠加”真的成立。通道之间存在幅度不一致、相位不一致或时延误差时增益会打折扣公式变成 N·ηη 是相干损失系数。工程上要求通道间相位误差控制在 10° 以内对应时延误差约 0.03 个周期。AWR2243 这类毫米波雷达芯片出厂做了通道校准但温度漂移和 PCB 加工公差仍会在高频段引入相位偏差实际处理时一般先采集已知角度目标的回波做校准矩阵。4.2 延时求和与 4D 毫米波雷达的角度维形成4D 毫米波雷达如 TI AWR2243 级联方案能同时输出距离、多普勒、方位角、俯仰角四维信息。距离维通过 FFT 处理快时间采样得到多普勒维通过对慢时间采样做 FFT 得到方位角和俯仰角维则是靠天线阵列的波束形成。以水平方位角为例如果接收天线是水平排布的均匀线阵那么对每个距离-多普勒单元沿着天线维做延时求和窄带下等效于 DFT 加权峰值对应的角度就是目标方位角。多个目标在同一距离-多普勒单元时延时求和波束会受限于波束宽度无法区分角度间隔小于半波束宽度的两个目标。此时常见做法是换用超分辨算法但代价是计算量上升且需要信源数先验。工程上先用延时求和做初测确认目标数目和大致方位再在局部角度范围内用 MUSIC 或压缩感知类方法精化是雷达信号处理仿真中比较稳妥的流程。4.3 AWR2243 雷达数据读取后的第一级波束处理TI 的 AWR2243 采集到的 raw ADC 数据格式通常是一个四维数组[采样点, chirp, 发射天线, 接收天线]。解出原始数据后后续处理链路的第一步是做距离-多普勒变换第二步才是跨接收天线的波束形成。代码逻辑示意如下// 伪代码: 距离-多普勒后处理 角度维波束形成 for (int chirp 0; chirp num_chirps; chirp) { // 1. 对每个chirp的ADC采样做距离FFT (快时间维) range_fft[chirp][rx] fft(adc_data[chirp][rx], range_fft_size); // 2. 对同一距离bin的多chirp数据做多普勒FFT (慢时间维) for (int r 0; r range_fft_size; r) { doppler_fft[r][rx] fft(range_fft[:, r][rx]); } } // 3. 对每个距离-多普勒单元, 沿RX天线维做延时求和 for (int r 0; r range_bins; r) { for (int d 0; d doppler_bins; d) { for (int theta 0; theta num_scan_angles; theta) { // 每个角度下, 计算各天线相位补偿并累加 complex_sum 0; for (int rx 0; rx num_rx; rx) { phase 2*PI * (rx * rx_spacing * sin(theta_rad)) / wavelength; complex_sum doppler_fft[r][d][rx] * exp(-j * phase); } angle_spectrum[r][d][theta] abs(complex_sum) * abs(complex_sum); } } }注意这里用exp(-j * phase)做相位补偿本质是窄带延时求和的频域形式。设置扫描角度时需要在意两点角度分辨率大约为 λ/(N·d·cosθ)在法线方向θ0最优偏离法线后分辨率下降视角范围受阵元方向图限制一般不超过 ±60°。num_scan_angles取 128 还是 256看角度分辨率和实时性要求的折中我一般先在 MATLAB 里用 1° 步进离线验证一遍再在 C674x DSP 里降采样到 2° 步进控制算力。5. 分数时延、通道失配校准与定向验证技巧延时求和波束形成最容易忽略的一个精度瓶颈是整数采样点近似。前文提到四舍五入引入最大半采样周期误差这个误差在窄带场景表现为一个额外相位在宽带语音场景则直接表现为波束指向偏移和旁瓣抬高。工程上补救方案有两种。第一种方案是频域相位旋转。把信号做 FFT乘以exp(-j*2π*f*τ)再逆变换回时域τ 可以是任意小数实现亚采样精度对齐。以 MATLAB 为例function y fractional_delay(x, tau, fs) % tau: 需要补偿的时延(秒), 可为小数 % fs: 采样率 N length(x); f (0:N-1) * fs / N; % 频率轴 phase exp(-1j * 2 * pi * f * tau); % 对应频域相位旋转 y real(ifft(fft(x) .* phase)); end这个方法的精度受 FFT 分辨率和频谱泄漏影响适合对整段信号做固定延迟补偿。第二种方案是时域 sinc 插值FIR 滤波器的抽头系数用截断的 sinc 函数做插值生成阶数在 9~17 之间较常用能实现逐点的小数延迟适合实时流式处理。我一般先用频域旋转验证算法逻辑再切到 FIR 方式做嵌入式实现。通道失配校准时一个实用的方法是采集一段正前方声源或雷达角反信号理想情况下所有通道时延相同输出应与单通道一致。实际结果如果出现零点说明存在幅度不一致或相位偏差幅度不一致通过每个通道除以自身 RMS 解决相位偏差用互相关估计出固定延迟差后补偿到对应通道。最后一个定向验证技巧延时求和输出的波束能量峰值位置不一定正好落在目标角度上。原因可能是目标信号在阵列孔径范围内不是平面波近场效应或者反射路径叠加导致等效波前畸变。验证时把目标从期望位置沿等距圆弧移动几个点位记录每个角度下波束输出电平把实测电平与仿真的波束图对比。两者主瓣位置偏差超过一个波束宽度时优先怀疑阵元位置和采样率配置而不是算法本身——这通常能节省半天到一天的调试时间。本文还有配套的精品资源点击获取