简介：面向信号处理方向的学习者与研究方向，这套 MATLAB 声源定位程序集围绕 TDOA 时延估计展开，涵盖广义互相关（GCC）与空间谱估计两条主流技术路线，并配套简单房间回响模型实验，可用于算法复现、课堂实验与仿真验证。压缩包共 20 个文件、约 1.67MB，其中 10 个 m 脚本承载各章节主程序，以及互相关、谱估计、房间冲激响应等核心函数；5 个 mat 文件提供实测或仿真数据；2 个 wav 语音信号和 2 份 docx 房间声学模型说明共同构成了完整的实验输入与理论对照。目前已有 368 人学习下载，说明它在教学与自学场景中具备一定参考价值。目录按简单房间回响模型、广义互相关定位、空间谱估计三段递进编排；你可以直接运行这些程序，对比不同算法在混响条件下的定位表现，也可以把它们当作二次开发的代码骨架，省去从零搭建仿真链路的时间。

1. 双麦录音跑出来峰在 0 点：先搞清这套程序的分层

实验室里最常见的一幕：拿一份声源定位算法的代码包，双击主脚本，图出来了，谱峰也挺漂亮；换成自己录的两路 wav，峰值死活贴在 0 附近，或者角度一会儿 20° 一会儿 -20°。问题通常不在 matlab 定位算法本身，而在于没分清这套程序里的两条主线——广义互相关只做「时延估计」，输出一个采样点差值；空间谱估计做「方位搜索」，输出一个角度或坐标。两者输入接口、帧长要求、符号约定都不一样，混着调必然翻车。

这份程序集合把两条路线都收进来了：GCC / GCC-PHAT 一族负责两两通道的时延，MUSIC 一类子空间方法负责窄带高信噪比下的角度谱，SRP-PHAT 负责宽带混响下的网格搜索。适合做声学定位算法验证、阵型对比、课程设计复现的人，也适合已经有一个能跑的基线、想找第二套实现做交叉验证的从业者。下面按「时延怎么算准 → 谱怎么搜对 → 工程怎么跑通 → 哪里会炸」的顺序拆。

2. 广义互相关为什么先做 PHAT 加权：从互谱到峰值的三步

2.1 两路信号到互功率谱：时延被搬进了相位

先把物理关系摆正。远场平面波、两个间距为 d 的麦克风，声源方向与阵列法线夹角 θ，声程差是 d·sinθ，对应时延 τ = d·sinθ / c，c 取 343 m/s（常温常压下的常用值，温湿度变化大时要做修正）。采样后这个时延变成 D = τ·fs 个采样点，通常不是整数。

频域里这件事特别好写。设 x1 是参考通道，x2 比 x1 晚到 D 个采样，那么

X2(f) = X1(f) · exp(-j2πfD/fs)

把两路的互功率谱写出来：

G12(f) = X1(f) · conj(X2(f)) = |X1(f)|² · exp(+j2πfD/fs)

时延完全落在相位项里，且是随频率线性增长的斜率。对 G12 做逆 FFT，理想无噪无混响条件下会在 τ = D 处出现一个冲激。这就是广义互相关的骨架：互谱 → 加权 → 逆变换 → 找峰。注意符号，X1 .* conj(X2)这一支在峰值索引为正时，表示 x2 滞后于 x1，也就是声源先到 1 号麦。写反了不会报错，只会让你的角度整体镜像，这类错误在第 5 章还会再提一次。

2.2 PHAT 加权：把幅度谱扔掉，只留相位

直接对 G12 做逆变换，在混响房间里几乎不可用。原因是房间冲激响应把每一路的幅度谱都染上了深谷和尖峰，互谱的幅度被通道特性主导，峰被抹平。PHAT 加权干的事情很粗暴：除以自己的模。

G12_phat(f) = G12(f) / |G12(f)| = exp(+j2πfD/fs)

分子分母一除，幅度信息全丢，只剩单位相量。这一步相当于把所有频点拉到同一权重，让那些信噪比低但相位仍可用的高频段也参与投票，峰变得更尖。代价是噪声大的频点同样被放大，所以实际实现里分母要加一个小的正则项abs(G12) + eps，并且最好在加权前把明显没信号的频带（比如 100 Hz 以下、接近奈奎斯特的高频）置零。

这也解释了一个常见困惑：为什么 PHAT 在低信噪比下反而不如朴素互相关。因为在那种场景里，噪声主导的频点被等权放大，投票被污染。工程上的折中是 SCOT 加权（除以两路自谱几何平均的平方根）或者对 PHAT 结果做频带限制。程序包里一般把这几种加权都留了开关，选型时不要只看论文里的默认值，要按你的信噪比和混响时间挑。

还有两个参数必须一起定，否则后面的谱估计全是白算：

• 帧长 N：决定了可分辨的最大时延。互相关结果的有效搜索范围是 ±N/2 个采样点，对应最大可测时延 N/(2·fs) 秒。要保证N/(2·fs) > d/c，也就是帧长必须覆盖阵元间的最大声程差，否则真实峰被卷到边界外面去了。
• 分析频带：语音建议 300~3500 Hz，宽带噪声源可以放到 200~6000 Hz。频带越宽，时延分辨率越高（理论上约 1/带宽），但高频的空间混叠风险也越大，见 5.5。

2.3 一个能跑的最小验证脚本：先仿真，再上真实录音

下面这段是不依赖任何工具箱的 GCC-PHAT 最小闭环，用来验证你手上的包符号约定对不对。先跑仿真，能准确恢复出设定时延，再去碰真实数据。

fs = 16000; % 采样率，语音常用值 c = 343; % 声速 m/s d = 0.05; % 阵元间距 5 cm N = 8192; % 帧长，取 2 的幂便于 FFT theta = 30 * pi/180; % 真实入射角（相对阵列法线） tau_true = d * sin(theta) / c; % 真实时延，秒 D = round(tau_true * fs); % 折算成采样点 s = randn(1, N); % 白噪声源，宽带便于验证 x1 = s; % 参考通道 if D >= 0 x2 = [zeros(1, D), s(1:N-D)]; % x2 滞后 D 点 else x2 = [s(1-D:N), zeros(1, -D)]; % 负时延交换补齐方向 end X1 = fft(x1); X2 = fft(x2); G12 = X1 .* conj(X2); G12 = G12 ./ (abs(G12) + eps); % PHAT 加权，只留相位 r = real(ifft(G12)); % 广义互相关函数 [~, idx] = max(r); lag = idx - 1; % 索引 1 对应 0 时延 if lag > N/2 lag = lag - N; % 循环移位修正到 (-N/2, N/2] end tau_est = lag / fs; fprintf('真实时延 %.3f ms，估计 %.3f ms，误差 %.3f 采样点\n', ... tau_true*1e3, tau_est*1e3, lag - D);

逻辑链条是：构造两路只差一个整数时延的信号 → 频域求互谱 → PHAT 归一化 → 逆变换找峰 → 索引换算回时延。几个参数的作用需要交代清楚：eps是防止静音段分母为零，取值 1e-10 量级就够，太大会把弱信号频点的相位信息一起压掉；idx - 1是因为 MATLAB 索引从 1 开始；lag > N/2那段修正不能省，否则负时延会显示成 N-|lag| 这种巨大正数，看起来像估计彻底失效。

跑通之后把s换成你录的 wav，用audioread读入、按通道拆开、各自去均值，其余不动。如果这时峰跑偏了，问题九成在数据端而不是算法端，直接跳到第 5 章。

提示：仿真验证时把D设成奇数个采样点，能看到峰值落在相邻两个 bin 之间、幅度被摊薄的现象。这就是后面要讲抛物线插值的原因，别把这种正常的幅度下降当成 bug。

3. 空间谱估计怎么落地：MUSIC 与 SRP-PHAT 的两条路

3.1 阵列流型：把几何关系压成一个向量

时延是两两之间的标量，空间谱要的是整个阵列的联合信息。均匀线阵、M 个阵元、间距 d、远场窄带信号频率 f0，方向 θ 的流型向量（steering vector）写成：

a(θ) = [1, exp(-j2πf0·d·sinθ/c), ..., exp(-j2πf0·d·(M-1)·sinθ/c)]'

这个向量的每一维就是一个阵元相对参考阵元的相位差，跟 2.1 节的时延是同一件事，只是写在了频域。所有空间谱方法——MUSIC、MVDR、ESPRIT——都是在这个流型上做文章，区别只在于用什么准则去度量「观测数据与某个方向的匹配程度」。

落地时最容易错的是符号和顺序。exp(-j...)换成exp(+j...)，角度谱会整体镜像；阵元坐标写成从左到右还是从右到左，同样镜像。程序包里如果 MUSIC 和 GCC 用的约定不一致，两条路线的结果就永远对不上，而各自的图看起来都「有峰」，极具迷惑性。建议在包里放一个 5 行的自检脚本：给定 +30° 的仿真源，两条路线都必须给出 +30°，差一个符号就当场暴露。

阵元间距还要满足空间采样定理。窄带中心频率 f0 下，不产生栅瓣（空间混叠）的条件是 d ≤ λ/2 = c/(2·f0)。1.5 kHz 对应 λ ≈ 22.9 cm，d ≤ 11.4 cm，看起来宽松；但如果你的宽带信号上到 6 kHz，λ/2 就只有 2.9 cm 了。这是宽带阵列最现实的约束：为了覆盖语音高频，阵元必须做得很密，密了又会导致低频段孔径不足、分辨率下降。程序里的d不是一个可以随便填的数，它跟你的分析频带是绑死的。

3.2 MUSIC 的三个决定因素：快拍数、信噪比、网格步长

MUSIC 的核心是对协方差矩阵做特征分解，把空间分成信号子空间和噪声子空间，然后拿流型向量去噪声子空间上投影，投影越小说明该方向越像信号：

P(θ) = 1 / (a(θ)' · En · En' · a(θ))

En 是噪声子空间。信号数 K 必须预先知道（或者用 MDL / AIC 估），K 估错一位，谱峰就多一个或少一个。

c = 343; fs = 16000; M = 4; d = 0.04; % 4 元线阵，间距 4 cm f0 = 1500; % 窄带分析中心频率 snap = 256; % 快拍数 theta_src = [-25, 10]; % 两个源的入射角（度） SNR = 15; % 信噪比 dB m = (0:M-1)'; A = exp(-1j*2*pi*f0*d*m*sind(theta_src)/c); % M×K 流型矩阵 S = (randn(2,snap) + 1j*randn(2,snap)) / sqrt(2); % 复高斯源 X = A*S; X = X + 10^(-SNR/20) * (randn(M,snap) + 1j*randn(M,snap)) / sqrt(2); R = X * X' / snap; % 样本协方差 [V, Lam] = eig(R); [~, ord] = sort(diag(Lam), 'descend'); V = V(:, ord); En = V(:, 3:end); % 前 2 个为信号子空间 scan = -90:0.5:90; P = zeros(size(scan)); for k = 1:numel(scan) a = exp(-1j*2*pi*f0*d*m*sind(scan(k))/c); P(k) = 1 / abs(a' * (En*En') * a); % 谱函数 end [~, pk] = findpeaks(P, 'SortStr','descend', 'NPeaks', 2); fprintf('估计角度: %.1f°, %.1f°\n', scan(pk));

几个参数怎么定，是从业者最关心的部分。快拍数决定了协方差的估计精度，经验值是至少 2~3 倍的阵元数，实际用 200 以上比较稳，低于 100 时谱峰位置会随帧随机漂移，这就是 5.2 要讲的抖动。信噪比低于 0 dB 时 MUSIC 的分辨率急剧退化，因为噪声子空间和信号子空间的特征值开始交叠，此时不是调参能救的，得靠增加阵元数或换 SRP 这类依赖时延而非子空间的方法。网格步长0.5° 对应 361 个扫描点，计算量可以忽略；但如果谱峰很窄（阵元多、信噪比高），0.5° 的粗网格会因为峰落在格点之间而低估峰值高度，做精细化时先粗扫再在峰附近 ±2° 内细扫到 0.05°。

注意：findpeaks属于信号处理工具箱。目标机器上没装的话程序会直接报未定义函数，改成手写三点局部极大值判断即可，逻辑只有几行，别为了一个找峰函数去折腾环境。

还有一条容易被忽略：上面这段是窄带写法，假定所有能量都在 f0 附近。真实语音是宽带的，常见做法是对每个 FFT 频点分别做一次 MUSIC，再把各频点的谱按能量加权平均（宽带 MUSIC），或者干脆给每个频点单独出一张谱看峰是否一致。只取单一 f0 去处理宽带语音，结果会随着这一帧在该频点上有没有能量而剧烈波动，看起来像「算法时好时坏」。

3.3 SRP-PHAT：宽带混响场景下的暴力搜索与算力预算

SRP-PHAT 的思路非常工程化，不跟你谈子空间，直接把所有麦克风对的 PHAT 互相关算好，然后对每一个候选声源位置，把各对麦克风「理论上应该出现时延」处的互相关值加起来，和最大的位置就是估计结果：

P(x, y) = Σ_{i<j} r_ij( τ_ij(x, y) · fs )

其中 τ_ij(x,y) = (|p_i - (x,y)| - |p_j - (x,y)|) / c，是近场模型下的精确时延，不再用远场的 sinθ 近似。

% 预计算所有麦克风对的 PHAT 互相关 pairs = nchoosek(1:M, 2); Nf = size(Xf, 2); r = cell(size(pairs,1), 1); for p = 1:size(pairs,1) G = Xf(pairs(p,1), :) .* conj(Xf(pairs(p,2), :)); G = G ./ (abs(G) + eps); r{p} = real(ifft(G)); end % 网格搜索 [Xg, Yg] = meshgrid(-2:0.05:2, 0.2:0.05:2.5); P = zeros(size(Xg)); for p = 1:size(pairs,1) i = pairs(p,1); j = pairs(p,2); di = sqrt((Xg - px(i)).^2 + (Yg - py(i)).^2); dj = sqrt((Xg - px(j)).^2 + (Yg - py(j)).^2); lag = round((di - dj) / c * fs); idx = mod(lag, Nf) + 1; % 循环索引，需与帧长对齐 P = P + r{p}(idx); end [~, imax] = max(P(:));

px、py是各麦克风的实际坐标向量，必须用实测坐标，用理论值会引入系统性偏差；这是 3.1 里流型几何的另一种表达。lag用round取整是为了能查表，代价是约半个采样点的量化误差，对应 343/(2·16000) ≈ 1 cm 的距离误差，普通场景够用，要高精度就把相邻两个 bin 做线性插值。

算力是这条路线的真门槛。网格 81×47 ≈ 3800 个点，6 对麦克风，单帧就是 2.3 万次查表累加，MATLAB 里跑一帧几十毫秒，实时性勉强。想上实时，常见做法有三条：粗网格先定位、再在候选区域做二级细搜；把查表改成矩阵化，一次算完所有网格点（内存换时间）；或者把搜索空间从二维降到一维（远场只搜方位角），点数直接砍掉一个量级。第一帧就要想清楚你的预算是多少，别等跑完 200 帧实验才发现要一晚上。

3.4 两条路线怎么选：别拿窄带方法的结论去评判宽带场景

选型其实取决于三个问题，不是「哪个更先进」。

双麦、只是要一个方位角，GCC-PHAT 加几何换算就够，别上谱估计。四元以上的小阵列、要分辨两个相近的源、信噪比尚可，MUSIC 的分辨率优势明显。房间混响重、声源是宽带、要的是坐标而不仅是角度，SRP-PHAT 更稳，代价是算力。

实际项目里更常见的组合是：用 SRP-PHAT 出粗定位，再用 MUSIC 在局部做精化；或者用 GCC 给出每对的时延和置信度，置信度低的对直接剔除，剩下的做最小二乘解算。程序包把几种方法并列放，是为了让你能对比，不是让你在一条流水线里全用一遍。

4. 从算法文件到可复现实验：调用链、参数表与批量跑法

4.1 调用链：入口脚本、算法函数、数据接口三层

这类程序包到手，别急着改算法。先花十分钟确认三层结构在哪里，后面所有改动都依赖它。典型组织方式如下（具体文件名以你拿到的包为准，但层次划分大同小异）：

project/ ├── main_gcc.m # GCC 路线的入口脚本：读数据 → 算时延 → 出图 ├── main_music.m # 空间谱路线入口：读数据 → 谱估计 → 出图 ├── main_srp.m # 网格搜索入口 ├── algo/ │ ├── gcc_phat.m # 核心：两路信号 → 时延（含加权开关） │ ├── music_spectrum.m # 核心：数据矩阵 + 流型 → 谱 │ └── srp_grid.m # 核心：多通道 + 麦坐标 → 代价图 ├── utils/ │ ├── framing.m # 分帧、加窗、去均值 │ └── peak_pick.m # 找峰 + 抛物线插值 └── data/ # 录音或仿真数据

关键判断标准只有一个：算法函数的输入输出是不是纯数值。如果gcc_phat.m进去就要读文件、出来还顺手画图，那这个包几乎没法做批量实验，第一步就是把它拆干净——读数据留在入口脚本，算法函数只吃数组、只吐数值。这个动作花二十分钟，能省掉后面几十次改脚本的时间。

数据接口有三件事必须统一，包与包之间对不上多半是这里的锅：通道顺序（哪个是参考通道）、采样率（不同来源的录音混在一起会导致时延换算全错）、以及数据是行向量还是列向量（MATLAB 里fft默认按列，行向量进来会算出一堆莫名其妙的矩阵结果）。建议在入口脚本开头就做一次强制规整：

[x, fs] = audioread('rec.wav'); x = x - mean(x, 1); % 逐通道去均值，去掉直流偏置 if size(x,1) < size(x,2) x = x.'; % 统一成「每列一个通道」 end

去均值这一步不是可选项。直流偏置在频域是 0 Hz 处的一个巨大冲激，PHAT 归一化之后它照样占权重，会让互相关峰被拉向 0 时延，这正是 5.1 的第一个成因。

4.2 关键参数表：把这几个数固定下来再谈复现

调参之前先把参数和物理量对应上，这张表建议直接抄进你的实验记录：

最容易埋雷的是 d 和 N 之间的联动关系：间距加大能提高角度分辨率，但最大时延 d/c 也随之变大，帧长不够就直接卷绕。举个具体数字，d = 0.2 m 时最大时延约 0.58 ms，16 kHz 下是 9.3 个采样点，看着不大；但互相关的搜索范围是 ±N/2，N = 1024 时只有 ±32 个采样点，一旦混响导致峰偏移，很容易撞到边界。这也是为什么长帧在混响场景下更稳：不是因为算法变了，而是搜索窗口变宽了。

4.3 批量扫描与结果落盘：把「看起来对」变成「可复现地对」

单次实验的结论没有说服力。参数扫描是这类程序包最该被用起来的场景——换信噪比、换混响时间、换阵元间距，看误差曲线的拐点在哪。

snr_list = -5:5:20; err = zeros(size(snr_list)); for k = 1:numel(snr_list) tmp = zeros(1, 50); % 每个条件做 50 次蒙特卡洛 for trial = 1:50 [x1, x2, tau_true] = sim_pair(fs, d, c, theta, snr_list(k), N); tau_hat = gcc_phat(x1, x2, fs); % 纯数值接口，方便循环调用 tmp(trial) = abs(tau_hat - tau_true) * 1e6; % 误差换算成微秒 end err(k) = mean(tmp); end save('sweep_snr.mat', 'snr_list', 'err');

sim_pair是自己写的仿真函数，负责生成一对带噪、带指定时延的信号；gcc_phat就是 2.3 节那个函数封装成纯数值接口后的样子。这段脚本的价值在于：每个条件 50 次取均值，把单次实验的偶然性压下去，err随snr_list的变化曲线才能说明问题——通常能看到某个信噪比以下误差陡增，那个拐点就是你这份数据的可用边界。

结果落盘时把参数一起存进去，别只存误差数组。我一般会额外存一个结构体，记录 fs、d、c、N、频带、算法名、随机种子。后面画图或者别人复现时，没有这些参数的数字就是一堆没法解释的数。50 次 × 6 个条件的循环，GCC 路线几十秒就跑完了；换成 SRP 的网格搜索，同样的实验规模会到小时级，这时候更要在跑之前把条件数、每条件次数、网格大小乘一遍，估个总时长再开跑。

5. 避坑排查：五类翻车现场与对应解法

5.1 现象：互相关峰永远贴在 0 附近，怎么调都动不了

原因：按出现频率排，第一位是两路数据其实是同一路（录音设备选错了通道，或者读文件时只取了第一列复制两次），这时候真实时延为 0，峰当然在 0。第二位是直流偏置没去，0 Hz 分量在 PHAT 归一化后仍占据权重，把峰往零时延拽。第三位是eps取得太大（比如 1e-3），分母被这个常数主导，加权退化成普通互相关，混响下的宽峰会盖住真实峰。

解决：先做一次「同源检测」——把 x1 和 x2 的相关系数算出来，如果接近 1，说明两路根本是同一个信号，回去查录音配置。然后强制去均值 + 300 Hz 高通，把eps降到 1e-10 量级。最后打印互相关曲线（不只是找峰），看看是真没有峰还是峰被埋了，这一步能省掉大量瞎猜。

5.2 现象：MUSIC 谱峰位置随帧剧烈抖动，相邻两帧能差 15°

原因：快拍数不足导致样本协方差估计方差大，特征值和特征向量都不稳。另一个常被忽略的原因是相干源——房间多径让直达声和反射声高度相关，协方差矩阵秩亏，信号子空间维数塌缩，噪声子空间被污染，谱峰开始乱跑。第三个原因是信号数 K 估错，把噪声当成信号，多出来的「源」在两个真实峰之间游走。

解决：快拍数先加到 300 以上看是否收敛，不收敛就别在这个方向上继续调参。相干源用前后向空间平滑（均匀线阵下把 M 个阵元分成若干重叠子阵，子阵协方差取平均）解相干，代价是有效孔径缩小。K 用 MDL 准则估，别手填。另外对协方差做对角加载（R = R + σ·eye(M)，σ 取噪声特征值均值的量级）能明显改善低快拍下的稳定性，这是工程上最省事的一招。

5.3 现象：角度估出来是真实值的镜像，+30° 变成 -30°

原因：符号约定不一致。流型里写成exp(+j...)而 GCC 那边用conj的约定，两条路线的角度定义正好相反。也可能是物理上的：麦克风阵列的坐标定义（1 号在左还是右）跟你代入几何公式时的假设反了。这类错误最阴险，因为谱图漂亮、峰也尖锐，只是位置镜像，不拿已知角度去验证根本发现不了。

解决：做一个符号自检脚本，放一个已知 +30° 的仿真源进去，GCC 给出的时延符号、MUSIC 给出的峰位置、SRP 给出的坐标偏移三者必须一致。不一致就逐个改符号，改到三个都对上为止。这个脚本写好之后放进包里长期保留，以后每次换阵列、换代码分支都先跑一遍。

5.4 现象：混响房间内单帧定位点满图乱跳，取平均之后又偏到墙角

原因：单帧估计把最强的那个互相关峰当成直达声，混响强的时候最强峰经常来自一次反射，于是定位点跳到镜像位置。直接对坐标取算术平均更糟——双峰分布下均值落在两个峰中间的空地上，物理上根本不存在声源。

解决：不要平均坐标，要做投票。把连续若干帧的估计结果累计成直方图（角度直方图或坐标网格），取票数最高的 bin，这是 SRP 类方法天然在做的事。同时给每帧加置信度筛选：互相关峰均比（峰值 / 次峰或峰值 / 均值）低于阈值的帧直接丢弃。实测中这两步组合能把乱跳压掉大半，比调任何算法参数都有效。

5.5 现象：换了麦克风间距，整套结果全崩，连之前能对上的角度也错了

原因：空间混叠。间距超过最高分析频率对应的半波长（d > c/(2·f_max)）之后，不同方向会对应相同的相位差，谱上出现多个等高峰（栅瓣），算法随机挑一个。另一个原因是帧长没跟着改，d 变大后最大时延 d/c 超出了 N/(2·fs) 的搜索窗口，峰被循环卷绕到错误位置。

解决：先按最高频率算出 d 的上限，再决定阵型；如果必须用大间距（为了低频分辨率），就把分析频带上限压下来，或者换非均匀阵列打散栅瓣。帧长和 d 要同步改：N ≥ 2·fs·d/c是硬约束，留一倍余量更稳。改完参数重新跑一次 2.3 节的仿真验证，别直接上真实数据。

6. 进阶：用已知角度做一次闭环校准，把「看起来对」变成「量得准」

真实系统里，误差很少是算法本身带来的，更多是几何和声速的系统性偏差：阵元间距标称 5 cm，实际装出来是 4.8 cm；声速按 343 算，实测环境可能是 338。这类偏差在角度上表现为「小角度还好、大角度偏得厉害」，因为它正比于 sinθ。

校准办法很朴素：把声源固定在若干个已知角度上（比如 -60° 到 +60° 每 10° 一个点），每个角度录一段，用 GCC-PHAT 估时延，然后跟理论时延做一次线性回归。

theta_list = -60:10:60; % 已知声源角度 tau_est = zeros(size(theta_list)); for k = 1:numel(theta_list) [x1, x2] = load_pair(theta_list(k)); % 该角度下的一对录音 tau_est(k) = gcc_phat(x1, x2, fs); % 纯数值接口 end tau_th = d * sind(theta_list) / c; % 理论时延 p = polyfit(tau_th, tau_est, 1); % tau_est ≈ p(1)*tau_th + p(2) c_eff = c / p(1); % 有效声速（含间距标定误差） res = tau_est - polyval(p, tau_th); rms_res = sqrt(mean(res.^2)); % 残差 RMS，单位秒 fprintf('等效声速 %.2f m/s，固定偏置 %.2f us，残差 RMS %.2f us\n', ... c_eff, p(2)*1e6, rms_res*1e6);

怎么读这三个数，是这套流程的关键。斜率 p(1)偏离 1，说明 d/c 这个比值有系统性误差，用c_eff回代到所有后续计算里，能把系统性偏差压掉；注意单靠时延只能定出 d/c 的比值，间距和声速谁错要靠卷尺和温湿度计另外确认。截距 p(2)是两路的固定延时，最典型的来源是两块采集设备的启动不同步或采样时钟漂移——这个值如果稳定在几十微秒量级，说明必须做通道同步，不然任何算法都救不回来。残差 RMS是这套系统的噪声底：它决定了你能分辨的最小角度变化，残差 20 μs 对应约 1.5 cm 的声程差，5 cm 间距下大约是 17°，也就是小于这个角度差的两个源在你的系统上根本分不开，再换算法也没用。

我第一次做这套校准时，把回归跑完发现斜率 1.13，当时以为是算法有问题，折腾了两天，最后拿卷尺一量，麦克风实际装在 4.4 cm 的位置，标称值一直写的是 5 cm。从那以后我每次碰新的阵列，都强制先走一遍这个闭环——量间距、跑角度扫描、看三个数、把有效声速写进配置文件，再谈后面的算法调优。省下来的不是调试时间，是不停怀疑算法的那些夜晚。

希望帮到你。

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