基于被动时反的 ac-MVDR 水下声源定位matlab代码

以下是基于被动时反的 ac-MVDR 水下声源定位的 MATLAB 代码示例：

% 清空工作空间
clear all;

% 加载声源和阵列数据
load('source_data.mat');
load('array_data.mat');

% 声源和阵列参数
fs = 48000;        % 采样频率
c = 1500;          % 声速
N = 8;             % 阵列元素数量
d = 0.5;           % 阵列元素间距
theta = 30;        % 水平方向声源角度
phi = 45;          % 垂直方向声源角度

% 构造波束形成矩阵
theta_range = -90:90;
phi_range = -90:90;
P = zeros(length(theta_range), length(phi_range));
for i = 1:length(theta_range)
    for j = 1:length(phi_range)
        a = steering_vector(N, d, fs, c, theta_range(i), phi_range(j));
        P(i,j) = a' * inv(R) * a;
    end
end

% 计算波束形成权重向量
desired_doa = [theta, phi];
a_desired = steering_vector(N, d, fs, c, desired_doa(1), desired_doa(2));
w = inv(R) * a_desired / (a_desired' * inv(R) * a_desired);

% 计算 ac-MVDR 权重向量
noise_doa = [-30, 60];    % 噪声方向
a_noise1 = steering_vector(N, d, fs, c, noise_doa(1), 0);
a_noise2 = steering_vector(N, d, fs, c, 0, noise_doa(2));
A_noise = [a_noise1, a_noise2];
V = A_noise' * inv(R) * A_noise;
w_ac_mvdr = inv(A_noise' * inv(R) * A_noise) * A_noise' * inv(R) * a_desired / (a_desired' * inv(R) * a_desired);

% 计算各个方向的 DOA 估计值
Rmvdr = zeros(N, N);
for i = 1:length(theta_range)
    for j = 1:length(phi_range)
        a = steering_vector(N, d, fs, c, theta_range(i), phi_range(j));
        Rmvdr = Rmvdr + P(i,j) * (a * w) * (a * w)';
    end
end
[~, doa_est] = max(diag(abs(Rmvdr)));
theta_est = theta_range(doa_est);

Racmvdr = zeros(N, N);
for i = 1:length(theta_range)
    for j = 1:length(phi_range)
        a = steering_vector(N, d, fs, c, theta_range(i), phi_range(j));
        Racmvdr = Racmvdr + P(i,j) * (a * w_ac_mvdr) * (a * w_ac_mvdr)';
    end
end
[~, doa_acmvdr_est] = max(diag(abs(Racmvdr)));
theta_acmvdr_est = theta_range(doa_acmvdr_est);

% 显示结果
disp(['水平方向声源真实角度：', num2str(theta)]);
disp(['水平方向声源估计角度（MVDR）：', num2str(theta_est)]);
disp(['水平方向声源估计角度（ac-MVDR）：', num2str(theta_acmvdr_est)]);

% 画出波束形成图像
figure();
imagesc(theta_range, phi_range, 10*log10(P));
xlabel('水平方向角度（度）');
ylabel('垂直方向角度（度）');
title(['波束形成图像（DOA：', num2str(theta), ',', num2str(phi), '）']);
colorbar;

% 声源信号波束形成
s_beamform = a_desired' * array_data;

% MVDR 估计
s_mvdr = w' * array_data;

% ac-MVDR 估计
s_acmvdr = w_ac_mvdr' * array_data;

% 显示结果
figure();
subplot(3,1,1);
plot(s_beamform);
title('波束形成');
subplot(3,1,2);
plot(s_mvdr);
title('MVDR 估计');
subplot(3,1,3);
plot(s_acmvdr);
title('ac-MVDR 估计');

其中 `source_data.mat` 是包含声源信号数据的 `.mat` 文件，`array_data.mat` 是包含阵列信号数据的 `.mat` 文件。`steering_vector` 函数用于计算阵列的空间滤波器系数，`R` 是阵列信号协方差矩阵，可以通过阵列信号数据计算得到。代码中分别计算了 MVDR 和 ac-MVDR 的波束形成权重向量，并根据波束形成权重向量计算了各个方向的 DOA 估计值。最后将声源信号经过波束形成、MVDR 和 ac-MVDR 估计，可以得到三个不同的水平方向声源信号估计结果，可以通过比较估计结果和真实声源信号进行定位误差分析。