声源定位算法4----CLEAN-SC（2）

最新推荐文章于 2026-08-21 12:37:40 发布

原创最新推荐文章于 2026-08-21 12:37:40 发布 · 920 阅读

22 ·

本内容遵循CC 4.0 BY-SA版权协议

收录于

声源定位算法

上篇博客声源定位算法4----CLEAN-SC（1）-CSDN博客中，推导了Clean-SC算法，这篇文章我将推导该算法工程实现，同时使用matlab仿真实现该算法。

CLEAN-SC：标准公式链 → 工程实现对照

1) 频域互谱矩阵（CSM）

理论：在每个频点 f 上，构造阵列复声压向量 p(f)，互谱矩阵

工程实现：用频带内 FFT 片段做样本平均得到 R，并把它当作 C 使用：

不同点：理论是期望 E{ ⋅ }，工程用有限快拍/频带平均。

2) DAS/CBF 的声源图（dirty map）

对扫描点 r，定义导向/聚焦向量：

理论：DAS/CBF 的 dirty map 为

工程实现：取幅值：

不同点：理论里应为实数非负，工程上数值误差/归一化/噪声会导致小的负值或复数残留，代码需要取绝对值保证非负可视化。

3) CLEAN-SC 的“残差 CSM”迭代框架

理论：初始化残差矩阵

第 k次迭代先用残差矩阵形成新的 dirty map：

工程实现：在进入 CLEAN-SC 前只计算了一次，然后每次迭代都用 PSF 扣除法更新，而不是每次重新用回算整张图。这样可以节省很多计算时间。

这是一个重要工程差异：

理论形式：每次迭代都用新的重新计算（代价很大）
工程形式：用“从 dirty map 中扣 PSF”的方式更新（更快）

这两者在理想条件下是等价思路，工程里常用后者。

4) 寻找最大谱峰并记录

理论：选取当前最大峰位置

并以进行累加：

工程实现：与理论相同，每轮找 Pmax（即） 和最大谱峰位置 (kx0,ky0)，代码示意如下：

Q(kx0,ky0)=Q(kx0,ky0)+Pmax*phi;

5) 相干场估计与 rank-1 剥离

这是 CLEAN-SC 的核心。

理论：令峰值点导向矢量

CLEAN-SC 的相干场向量常写成：

然后从残差 CSM 中剥离掉该相干分量（rank-1）：

工程实现：与理论相同，代码示意如下：

h = D*wmax./Pmax;
D = D - Pmax.*(h*h');

6) 用 PSF 从 dirty map 中扣除

理论：从残差矩阵出发，峰值分量对任意扫描点 r 的贡献可以写成：

然后更新 dirty map：

工程实现：把所有网格的导向矢量缓存成矩阵 W，一次性计算所有点：

不同点：理论上也可以每次从重新计算；这里用 PSF 扣除更新，从而避免重复计算整张图，这是 CLEAN 类算法最常见的工程优化。

7) 迭代停止条件

常用停止条件是满足其一即可：

达到最大迭代次数
当前峰值低于阈值
残差能量低于阈值

工程里为了数值稳定，经常加入“非负约束”和“能量上限/单调性保护”，这些不是理论必要条件，但很常用。

8) 工程输出 : 采用“CLEAN + 非负残差”组合

工程上把非负残差图补回去：

工程对应：

Pclean_temp = abs(Q)./phi + max(Pcbf,0);

与“严格 PSF restore”最大的差别：

严格 restore：
把 Q 再与 PSF 做卷积，得到物理一致的 restored map。
工程输出：
直接把 “干净声图” 加上“非负残差图”，好处是更稳更像声图、对扩展源更友好；代价是残差里仍可能包含部分旁瓣/未完全剥离的脏能量。

9) 宽带输出：逐频点求后累加

理论宽带融合：

也可以加频率权重 wk，比如 A-weighting 或能量归一化：

CLEAN-SC 在每个频点上以 CSM 为输入，先生成 CBF 脏图并定位最大谱峰；随后用相干场构造 rank-1 分量并从残差 CSM 中剥离，同时用快速更新脏图；最终工程上常用输出稳定的声图，宽带结果通过对频点累加获得。

仿真：

采用matlab实现宽频信号的CLEAN-SC算法。

算法模型采用近场，麦克风阵列采用72通道的多臂螺旋（与前面CBF的一致）。

麦克风阵列

仿真参数设置（与CBF几乎一致）：

%% 参数设置
c = 343;     % 声速 (m/s)
N = 2048;    % 所需信号长度
Fs = 64e3;   % 采样频率 (Hz)
fc1 = 3500;  % 第一个声源频率 (Hz)
fc2 = 4500; % 第二个声源频率 (Hz)
SNR = 10;    % 信噪比 (dB)
t = 0:1/Fs:1; % 生成 1 秒长的时间序列

%% 声源位置
source_location1 = [-0.45, 0.6, 2];  % 第一个声源位置
source_location2 = [0.15, -0.2, 2]; % 第二个声源位置

根据仿真参数设置，计算得到两个仿真信号以及最重要的时延信息，时延计算代码如下（网格自定义）：

%% 计算距离差矩阵，模拟网格
Max_x = 1.3;
Max_y = 0.75;
pixel = 0.04;
Site_matrix = -Max_x:pixel:Max_x;
Site_matriy = -Max_y:pixel:Max_y;
Length_matrix = length(Site_matrix);
Length_matriy = length(Site_matriy);

tau = zeros(Length_matrix, Length_matriy, BSN);
for ky = 1:Length_matriy
    site_ky = Site_matriy(ky);
    for kx = 1:Length_matrix
        site_kx = Site_matrix(kx);
        % 探测点到各麦克风的时延
        DS = [site_kx, site_ky, 2];   % 探测点位置
        for kki = 1:BSN
            ri = sqrt((BS(1, kki) - DS(1))^2 + (BS(2, kki) - DS(2))^2 + (BS(3, kki) - DS(3))^2) - sqrt(DS(1)^2 + DS(2)^2 + DS(3)^2);
            tau(kx, ky, kki) = ri / c;
        end
    end
end

得到时延信息，进行成像：

%% 功率谱成像 —— CLEAN-SC（逐频点）
tic

X_ki = zeros(BSN, N);
for ki = 1:BSN
    X_ki(ki, :) = fft(s(ki, :)); % FFT
end

fre_step  = 100;
fre_range = 3000:fre_step:5200;

% ========= CLEAN-SC 参数========
phi = 0.6;      
gen = 25;       

Pclean = zeros(Length_matrix, Length_matriy);

for kf = 1:length(fre_range)

    % 当前频率
    fc = fre_range(kf);
    w  = 2 * pi * fc;

    % ===== 频段协方差矩阵 =====
    Tna = round(fc * N / Fs) - round(fre_step * N / Fs) / 2;
    Tnb = round(fc * N / Fs) + round(fre_step * N / Fs) / 2;
    Tna = round(max(1, Tna));
    Tnb = round(min(N, Tnb));
    X0  = X_ki(:, Tna:Tnb);
    R   = X0 * X0' / (Tnb - Tna);

    % ===== 1) 计算 CBF 脏图 + 缓存导向矢量 W =====
    Kgrid = Length_matrix * Length_matriy;
    W     = zeros(BSN, Kgrid);                  % 缓存导向矢量
    Pcbf  = zeros(Length_matrix, Length_matriy);

    for kx = 1:Length_matrix
        for ky = 1:Length_matriy
            K = (kx-1)*Length_matriy + ky;      
            b = exp(-1i * w * squeeze(tau(kx, ky, :)));
            W(:, K) = b;
            Pcbf(kx, ky) = abs(b' * R * b);      % CBF dirty map
        end
    end

    % ===== 2) CLEAN-SC 迭代 =====
    D = R;                                      % 迭代用的“残差”CSM
    Q = zeros(Length_matrix, Length_matriy);   

    prevSum = sum(sum(abs(Pcbf)));
    allowSum = prevSum * 2;                     
    kgen = 1;

    while (kgen <= gen) && all(Pcbf(:) >= 0) && (prevSum < allowSum)

        % --- 寻找最大谱峰 ---
        [Pmax, idx] = max(Pcbf(:));
        [kx0, ky0]  = ind2sub(size(Pcbf), idx);

        Q(kx0, ky0) = Q(kx0, ky0) + Pmax * phi;

        % --- 最大谱峰处导向矢量 ---
        Kmax = (kx0-1)*Length_matriy + ky0;
        wmax = W(:, Kmax);

        h = D * wmax ./ Pmax;

        D = D - Pmax .* (h * h');

        % --- 计算并扣除“PSF 贡献” Ppsf
        Ppsf_vec = Pmax .* abs(W' * h).^2;               
        Ppsf = reshape(Ppsf_vec, [Length_matriy, Length_matrix])';
        Pcbf = Pcbf - Ppsf;

        % --- 更新停止判据用的能量 ---
        newSum = sum(sum(abs(Pcbf)));
        if newSum >= prevSum
            break;                                          % 防止数值导致能量不降
        end
        prevSum = newSum;

        kgen = kgen + 1;
    end

    % ===== 3) 该频点的 CLEAN-SC 声图=====
    Pclean_temp = abs(Q) ./ phi + max(Pcbf, 0);

    % ===== 4) 宽带累加 =====
    Pclean = Pclean + Pclean_temp;

end

toc

[x_location, y_location] = find(Pclean == max(max(Pclean)));

仿真结果：