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

以下是基于被动时反的 ac-MVDR 水下声源定位的 Matlab 代码，仅供参考：

clear all;
close all;

% 读取数据
load('data.mat');

% 声源位置
s_pos = [1500, 1000]; 

% 传感器位置
r1_pos = [0, 0];
r2_pos = [0, 3000];
r3_pos = [4000, 3000];
r4_pos = [4000, 0];

% 采样频率
fs = 8000;

% 音速
c = 1500;

% 信号长度
N = length(data);

% 声源到各传感器的距离
d1 = norm(s_pos - r1_pos);
d2 = norm(s_pos - r2_pos);
d3 = norm(s_pos - r3_pos);
d4 = norm(s_pos - r4_pos);

% 接收信号
r1 = data(:, 1);
r2 = data(:, 2);
r3 = data(:, 3);
r4 = data(:, 4);

% 平移对齐
r1 = r1(1: end - round(d1 / c * fs));
r2 = r2(round(d2 / c * fs) + 1: end);
r3 = r3(round(d3 / c * fs) + 1: end);
r4 = r4(round(d4 / c * fs) + 1: end - round(d1 / c * fs));

% 信号长度
N = length(r1);

% STFT参数
L = 256;
M = round(L / 4);
K = L - M;
NFFT = L;
nframes = floor((N - L) / M) + 1;

% STFT
X1 = zeros(NFFT, nframes);
X2 = zeros(NFFT, nframes);
X3 = zeros(NFFT, nframes);
X4 = zeros(NFFT, nframes);

for i = 1: nframes
    n1 = (i - 1) * M + 1;
    n2 = n1 + L - 1;
    
    X1(:, i) = fft(r1(n1: n2), NFFT);
    X2(:, i) = fft(r2(n1: n2), NFFT);
    X3(:, i) = fft(r3(n1: n2), NFFT);
    X4(:, i) = fft(r4(n1: n2), NFFT);
end

% 平均谱密度
S11 = mean(abs(X1) .^ 2, 2);
S22 = mean(abs(X2) .^ 2, 2);
S33 = mean(abs(X3) .^ 2, 2);
S44 = mean(abs(X4) .^ 2, 2);

% 协方差矩阵
R = zeros(4, 4, NFFT);

for k = 1: NFFT
    R(:,:,k) = [X1(k, :); X2(k, :); X3(k, :); X4(k, :)] * ...
        [X1(k, :); X2(k, :); X3(k, :); X4(k, :)]' / nframes;
end

% ac-MVDR
theta = -90: 1: 90;
p = zeros(length(theta), 1);
for i = 1: length(theta)
    a1 = exp(-1i * 2 * pi * fs * (0: NFFT - 1)' * d1 / c * sin(theta(i) / 180 * pi));
    a2 = exp(-1i * 2 * pi * fs * (0: NFFT - 1)' * d2 / c * sin(theta(i) / 180 * pi));
    a3 = exp(-1i * 2 * pi * fs * (0: NFFT - 1)' * d3 / c * sin(theta(i) / 180 * pi));
    a4 = exp(-1i * 2 * pi * fs * (0: NFFT - 1)' * d4 / c * sin(theta(i) / 180 * pi));
    a = [a1, a2, a3, a4];
    p(i) = 1 / (a' * inv(R(:,:,1) + R(:,:,2) + R(:,:,3) + R(:,:,4)) * a);
end

% 绘制谱图
figure;
plot(theta, 10 * log10(p / max(p)), 'linewidth', 2);
xlabel('角度/°');
ylabel('幅度/dB');
grid on;

需要注意的是，以上代码仅提供了一个基本的实现框架，实际应用中还需要根据具体情况进行一些参数调整和优化。