多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

固定波束形成:麦克风阵列语音增强的稳定基础与工程实践

固定波束形成:麦克风阵列语音增强的稳定基础与工程实践 简介一份面向语音增强与麦克风阵列信号处理学习者与开发者的固定波束形成MATLAB实现。代码以BeamF.m为核心完整呈现从阵列信号预处理、波束形成参数配置、各麦克风通道权重计算到目标方向语音信号合成的处理流程并含增强效果评估环节便于观察信噪比提升情况。该算法主要用于噪声环境下指定方向语音的拾取与干扰抑制可应用于智能语音交互、会议系统、机器人听觉等场景。对于正在理解延迟求和等固定波束形成算法或想快速搭建麦克风阵列增强原型的人来说能提供直观的代码参考。压缩包仅含1个m文件体积仅约1KB轻量但核心步骤齐全适合课程设计、实验验证或二次扩展。已有328人学习浏览属于短小精悍的入门级示例阅读时可重点关注权重系数计算与信号加权叠加部分这有助于后续掌握自适应波束形成等技术。1. 固定波束形成为什么先做固定波束再谈自适应在语音识别和通话降噪的项目里我经常看到这样的现象开发人员一开始就上自适应波束形成结果在真实噪声场中权重乱跳去混响效果还不如不启用的状态。其实先花半天把固定波束形成Fixed Beamforming跑通很多问题就解决了一大半。固定波束形成的思路很简单用一组固定的加权系数把麦克风阵列接收到的多路信号合成一路让目标方向来的语音尽量同相叠加非目标方向的噪声相消从而提升信噪比。它不需要统计估计没有迭代收敛CPU开销低非常适合麦克风阵列语音增强的硬件实现。这篇文章以BeamF.m这段MATLAB代码为线索拆开讲讲固定波束形成的原理、线阵模型、代码实现和在实际双麦克风阵列上的部署要点。2. 线阵模型与波束图固定波束形成的数学基础2.1 均匀线阵的导向矢量固定波束形成的核心是“相位对齐”。以最简单的M元均匀线阵为例麦克风间距为d声波以平面波形式从θ方向入射θ为与阵法线方向的夹角。第m个麦克风相对于参考麦克风接收到的时延差为τ_m (m-1)·d·sinθ/c其中c是声速工程上常取343 m/s。在频域该时延对应相位偏移e^(-j2πf τ_m)。因此从θ方向来的信号在M个麦克风上形成的响应向量就是通常说的导向矢量steering vectora(θ,f) [1, e^(-j2πfdsinθ/c), ..., e^(-j2πf(M-1)dsinθ/c)]^T理解了导向矢量固定波束形成器的设计就转化为寻找一组权重向量w使得w^H·a(θ_target,f)尽量大同时让非目标方向的响应尽量小。经典的延迟求和Delay-and-Sum波束形成器直接取w a(θ_target,f)/M每个麦克风信号在频域乘上共轭相位补偿后再相加目标方向信号相干叠加幅度放大M倍不相关的传感器噪声叠加后只放大√M倍因此阵列增益约为10log10(M) dB。对双麦克风来说是3dB对四元阵列是6dB这也是固定波束在麦克风阵列语音增强中稳定有效的原因。2.2 用MATLAB计算波束图在拆BeamF.m之前先写一个最小可运行的脚本用于直观观察波束图。波束图beam pattern就是|w^H·a(θ,f)|²关于角度的曲线它清楚地展示了主瓣方向和旁瓣水平。% beam_pattern_demo.m fs 16000; % 采样率 Hz M 4; % 麦克风数量 d 0.04; % 麦克风间距 4cm c 343; % 声速 m/s freq 1000; % 观察频率 Hz theta -90:0.5:90; % 角度扫描范围 theta_t 0; % 目标方向法线方向 % 构造目标方向导向矢量 a_t exp(-1j*2*pi*freq*d*sin(theta_t*pi/180)*(0:M-1)./c); w a_t / M; % 延迟求和权归一化保证目标增益为1 % 扫描各方向响应 resp zeros(size(theta)); for i 1:length(theta) a_scan exp(-1j*2*pi*freq*d*sin(theta(i)*pi/180)*(0:M-1)./c); resp(i) abs(w * a_scan); end figure; plot(theta, 20*log10(resp eps)); xlabel(入射角度 (degree)); ylabel(幅度响应 (dB)); grid on; title(4元均匀线阵波束图 1kHz);这段代码先生成目标方向的导向矢量a_t取其共轭转置再除以M作为权矢量w对应延迟求和波束形成。扫描循环里对每个入射角度计算导向矢量a_scan用w*a_scan得到该方向的增益最后转成dB画图。0度方向响应为0 dB其他方向低于0 dB凹陷的位置和深度由阵列孔径和频率决定。这里有一个容易忽略的参数约束当麦克风间距d大于最高频率对应波长的一半时波束图会出现与主瓣同样幅度的栅瓣导致噪声从栅瓣方向泄入。所以d的选择不能只看物理安装还要结合算法的最高工作频率。下表汇总了线阵常用参数对波束的影响后面配置BeamF.m时会直接用到。参数典型值对波束的影响注意点M麦克风数量2~8M越大主瓣越窄旁瓣数量增多通道不一致会破坏波束d麦克风间距0.02~0.05 m间距越大低频指向性越好高频可能出现栅瓣d 须小于 λ_min/2θ_target目标方向0°~±45°主瓣指向偏移旁瓣分布不对称超过±60°时波束畸变明显采样率 fs16k~48k决定可处理频率上限和帧长与ADC时钟绑定不能随意设置2.3 频率依赖与宽带波束形成上面画的波束图只在单一频率上有效。语音是宽带信号从300 Hz到3400 Hz甚至更高不同频率的波束宽度差异很大。举例来说1 kHz时一个4cm间距的四元阵主瓣宽度约30度到4kHz时主瓣缩到不到10度而低频段指向性很弱。如果只做一个频率的相位补偿整个频带上的噪声抑制效果会很不均匀。所以工程上固定波束形成普遍采用频域实现将多通道信号分帧加窗后做FFT在每个频点独立计算导向矢量和权重加权求和后再IFFT并重叠相加合成时域信号。这样既能正确处理分数时延又能在某些频段根据需要压低旁瓣。BeamF.m正是按照这个思路编写的。下一章直接进入函数代码看权重到底是怎么算出来的。3. BeamF.m代码解析权重计算与信号合成实现3.1 BeamF.m的函数结构与参数约定BeamF.m的典型入口是这样一个函数输入多通道信号矩阵、采样率、麦克风坐标、目标方向和频率范围输出增强后的单通道信号。为了能直接运行并对照调试我把它整理成下面的完整实现。这份代码兼顾了可读性和工程性适合在此基础上修改。function [y, W_all] BeamF(x, fs, mic_pos, theta_target, freq_range) % BeamF.m - 固定波束形成频域延迟求和实现 % 输入: % x 多通道信号, [M x samples] % fs 采样率, Hz % mic_pos 麦克风坐标, [M x 2] 或 [M x 3]单位米 % theta_target 目标声源方向度相对阵列法线 % freq_range 处理的频率范围, [f_low, f_high] Hz % 输出: % y 增强后的单通道信号, [1 x samples] % W_all 各频点权向量, [nfft/21 x M] [M, ns] size(x); c 343; x x - mean(x, 2); % 去直流 % 分帧参数32ms帧长50%重叠 frame_len round(fs * 0.032); hop round(frame_len / 2); nfft 2^nextpow2(frame_len); win hann(frame_len, periodic); num_frames floor((ns - frame_len) / hop) 1; theta_r theta_target * pi / 180; freqs (0:nfft/2) * (fs / nfft); low_bin find(freqs freq_range(1), 1); high_bin find(freqs freq_range(2), 1, last); if isempty(low_bin), low_bin 2; end if isempty(high_bin), high_bin nfft/21; end low_bin max(2, low_bin); high_bin min(nfft/21, high_bin); % 离线计算每个频点的权向量 W_all zeros(nfft/21, M); for k low_bin:high_bin f freqs(k); % 目标方向由 (sin(theta_r), cos(theta_r)) 决定 tau (mic_pos(:,1) * sin(theta_r) mic_pos(:,2) * cos(theta_r)) / c; a exp(1j * 2 * pi * f * tau); % 导向矢量 [M x 1] W_all(k, :) conj(a). / M; % 延迟求和权归一化 end % 频域逐帧处理 y_all zeros(1, num_frames * hop frame_len); win_sum zeros(1, num_frames * hop frame_len); for idx 1:num_frames start_pos (idx - 1) * hop 1; seg x(:, start_pos:start_posframe_len-1); seg seg .* repmat(win, M, 1); % 加入分析窗 X fft(seg, nfft, 2); % M x nfft Y zeros(1, nfft); for k low_bin:high_bin Y(k) W_all(k, :) * X(:, k); % 频点加权融合 end % 补齐共轭对称部分用于实数IFFT for k 2:nfft/2 Y(nfft - k 2) conj(Y(k)); end y_frame real(ifft(Y, nfft)); y_frame y_frame(1:frame_len); y_all(start_pos:start_posframe_len-1) ... y_all(start_pos:start_posframe_len-1) y_frame .* win; win_sum(start_pos:start_posframe_len-1) ... win_sum(start_pos:start_posframe_len-1) win.^2; end % 重叠相加增益归一化 win_sum(win_sum 0.001) 1; y y_all(1:ns) ./ win_sum(1:ns); y y / max(abs(y)) * 0.95; % 归一化输出幅度 end代码里的关键点要逐一说明。首先是low_bin和high_bin的查找find(freqs freq_range(1), 1)返回第一个大于等于下限频率的索引find(freqs freq_range(2), 1, last)返回最后一个小于等于上限的索引。这样可以把不关心的频段掏空避免极低频的直流漂移和超高频的栅瓣进入输出。频率默认从bin 2开始是因为bin 1是直流波束形成对直流没有意义。权重计算部分tau是每个麦克风相对于坐标原点的到达时间差mic_pos(:,1)*sin(theta_r) mic_pos(:,2)*cos(theta_r)是声源方向单位向量与麦克风位置的点积。注意这个公式对二维坐标成立当目标方向是法线方向(theta0)时只有y坐标参与延时当theta90度时只有x坐标参与。这种写法比直接两个麦克风做差更通用也方便扩展到三维坐标。W_all(k,:) conj(a)./M做的是共轭相位旋转乘到频域信号上正好把目标方向的相位对齐每个麦克风贡献1/M幅度实现无损相加。3.2 频域处理的几个工程细节分帧参数上我习惯用32ms帧长、50%重叠加Hann窗。这个参数组合在语音质量和延迟之间比较平衡。16k采样时帧长512点FFT也是512点低频分辨率达到31.25Hz对于波束权的频率变化足够细。window在分析时乘在时域信号上重叠相加时再乘一次所以功率归一化用win.^2来累计。这里很多人容易出错如果使用wsum wsum win输出信号会带有周期性调制听感像“哇哇”的鼓声因为重叠区域的幅度不是恒定增益。频点循环里Y(nfft-k2) conj(Y(k))这一步填充负频率部分。MATLAB的fft结果中正频率是索引2到nfft/2负频率部分是共轭对称。不补这一句ifft出来的是复数信号直接取实部会导致波形严重失真。另外一个细节Y(1)和Y(nfft/21)保持为零虽然会丢掉直流和奈奎斯特频率但对语音影响极小还能顺便消除传感器偏置带来的直流分量。3.3 参数怎么调频带、间距和通道数下表给出BeamF.m在高频上限出现栅瓣时的临界间距参考值可以按这个表快速判断当前阵列配置是否安全。最高处理频率 f_high半波长 λ/2声速343m/s麦克风间距上限3000 Hz5.72 cm建议d ≤ 5 cm4000 Hz4.29 cm建议d ≤ 4 cm8000 Hz2.14 cm建议d ≤ 2 cm如果你的麦克风间距超出上限最直接的办法是在BeamF.m中把freq_range(2)压到临界频率以下。比如d4cm时最高只能处理到4.3kHz电话语音完全够用要做宽频音乐录音就得缩小麦克风间距或改用自适应波束。M从4降到2方向性会明显变松但低频灵敏度损失不大。这解释了为什么很多智能音箱只用双麦克风也要坚持把两个麦克风拉开到4cm以上就是为了保证3kHz语音有足够的空间采样间隔。4. 双麦克风阵列部署es8311采集链路上的参数校准4.1 双麦克风阵列的几何配置固定波束形成在真实产品里最常见的载体是双麦克风阵列配合一颗低功耗音频编解码器典型的组合就是双麦克风阵列和es8311音频编解码器电路。两颗麦克风放在设备的同一块面板上间距取15mm到50mm方向性在前方形成一个“扇区”用于捕捉人的语音、抑制侧面和后面的噪声。BeamF.m里把麦克风坐标定义为mic_pos因此真实硬件的位置就映射成两个点。假设设备正面朝向用户双麦沿水平线排列则坐标可以写成[-0.02, 0; 0.02, 0]这时theta_target0对应正前方。在es8311这一侧通常通过I2S接口输出双通道音频数据到主控。es8311内部集成了两路ADC和MICBIAS麦克风被偏置到1.0V左右。采集回来的数据先进入算法前处理再做BeamF。因为es8311是立体声ADC两通道的采样是同一时钟相位一致性天然好这一点比用两颗独立codec的电路更可靠。这也是固定波束形成在es8311平台上容易落地的原因。4.2 从MATLAB权重到嵌入式C代码嵌入式处理器比如常见的DSP或带FPU的MCU无法直接跑MATLAB的矩阵计算。标准做法是在PC上运行BeamF.m把W_all导出为一组常量数组烧录到设备里。在线处理时只需要做FFT、复数乘加和IFFT。以16kHz采样512点FFT双麦为例每个频点需要两次复数乘加共256个频点主频在100MHz以下的MCU也能实时跑完。// beam_f_fixed.c // 双麦克风固定波束频域处理示意 // offline_weight[k][0] 表示麦克风0的复权offline_weight[k][1] 表示麦克风1 extern const float offline_weight[257][4]; // 实0,虚0,实1,虚1 void beam_f_apply(complex_t *X0, complex_t *X1, complex_t *Y, int nfft) { int k; for (k 0; k nfft/2; k) { Y[k].real offline_weight[k][0]*X0[k].real - offline_weight[k][1]*X0[k].imag offline_weight[k][2]*X1[k].real - offline_weight[k][3]*X1[k].imag; Y[k].imag offline_weight[k][0]*X0[k].imag offline_weight[k][1]*X0[k].real offline_weight[k][2]*X1[k].imag offline_weight[k][3]*X1[k].real; } // 剩余负频率部分在IFFT前补齐由DSP库完成 }这段C代码的每一行都对应BeamF.m里的Y(k) W_all(k,:) * X(:,k)。注意offline_weight的排列顺序是[实部0, 虚部0, 实部1, 虚部1]方便DSP的SIMD指令一次加载四个浮点。负频率部分一般由DSP库的逆FFT自动处理前提是只填充正频率或者算出负频率填入两种方式在代码风格上要提前约定好。从MATLAB导出权重时我一般直接打印成C数组fid fopen(weights.c,w); for k 1:size(W_all,1) fprintf(fid, {%f,%f,%f,%f},\n, ... real(W_all(k,1)), imag(W_all(k,1)), ... real(W_all(k,2)), imag(W_all(k,2))); end fclose(fid);这段生成代码不必多说关键是保持W_all的行数与FFT的bin数一致。如果嵌入式使用非标准FFT点数记得在BeamF.m里把nfft改成一致。4.3 实际校准的三个关键点固定波束形成对通道失配非常敏感。双麦克风阵列在 es8311 电路上看起来是两路完全相同的放大链路但实际元件容差、麦克风灵敏度差、PCB走线长度差都会让“固定”的相位补偿失效。所以我习惯在出厂前做以下三项校准校准项测试方法合格阈值通道增益失配播放1kHz正弦采集两通道RMS差值 0.5 dB通道相位差播放1kHz正弦用互相关估计时延差值 0.1 ms本底噪声安静环境静音采集两通道噪声级差 3 dB校准不是改硬件而是在算法输入端做数字补偿。增益失配可以在BeamF之前给幅度较小的通道乘一个系数相位差可以用一个分数延时滤波器把两路对齐。es8311的ADC直流偏移也要检查如果直流偏移不一致会直接进入FFT的直流bin虽然BeamF把它置零了但可能造成相邻bin泄漏所以最好在分帧前先做高通滤波或去均值。4.4 es8311电路上的相位一致性设计最后说电路层面的细节。双麦克风阵列的PCB走线从麦克风焊盘到es8311的MICIN引脚长度差应当尽量小。以16kHz为例声波在空气中传播1cm需要约29us在PCB走线中信号传播速度大约是光速的一半但走线长度差引起的电气延时会限制波束形成器的极限精度。好的设计是把两个麦克风到codec的走线做等长至少在5mm以内这样相位误差小于0.5度16kHz对波束图的影响可以忽略。es8311的电路还有两个常见坑一是MICBIAS电阻选大了导致麦克风偏置不稳表现为低频噪声抖动二是耦合电容容值太小会让低频响应塌陷影响300Hz以下的语音。固定波束形成对低频本来指向性就弱再把低频砍掉语音会变得单薄。一般用2.2uF或4.7uF的隔直电容低频截止点控制在20Hz以下。5. 用波束图快速验证固定波束的有效性5.1 验证波束主瓣是否对准目标拿到BeamF.m输出的W_all不要急着接到实车上先在MATLAB里画一幅多频波束图确认主瓣方向和旁瓣高度。下面这段代码取某个频点的权向量画出该频率下全角度的响应。% verify_beam.m k find(freqs 1000, 1); % 取1kHz附近的频点 w W_all(k, :).; % 当前频点的权向量 theta_scan -90:0.5:90; resp zeros(size(theta_scan)); for i 1:length(theta_scan) tau mic_pos(:,1)*sind(theta_scan(i)) mic_pos(:,2)*cosd(theta_scan(i)); a exp(1j*2*pi*freqs(k)*tau/c); resp(i) abs(w * a); end plot(theta_scan, 20*log10(resp eps));把这段代码跑完合格的标准是0度方向的幅度比其他方向至少高6dB主瓣两侧没有超过-10dB的旁瓣。如果画出来的图形在主瓣之外出现一个和主瓣等高的“鬼峰”说明当前频率已经超过栅瓣临界值需要降低freq_range(2)或减小麦克风间距。5.2 快速定位两类典型问题第一类问题是波束主瓣偏了。比如目标声源明明在正前方0dB峰值却落在了40度位置。绝大多数情况是theta_target的单位或mic_pos的坐标轴映射写错比如cosd和sind用反了或者坐标原点不在阵列中心。第二类问题是输出语音“嗡声”明显这通常是因为低频段权重保留了过大的旁瓣或者去直流没做好。检查方法很简单把频带从100Hz到200Hz单独跑一遍波束图如果这个频段响应接近全向属于正常如果连目标方向都被压制赶紧检查高通的转折频率是不是设到了300Hz以上。样机调试时建议先把波束图打在屏幕上让声源轮流从几个角度发声看到响应变化再入手调参数。本文还有配套的精品资源点击获取
返回列表