延迟求和波束形成算法（DAS）

function [DAS_result, PSF, hn,CSM] = DAS(N,z0,f,phi,rn,source,SNR)
% 这个代码实现了延迟并求和 (DAS) 算法
% 输入:
%    N: 每个维度的网格点数
%    z0: 来源距离
%    f: 成像频率
%    rn: 麦克风阵列的坐标
%    source: 来源的x，y位置
%    SNR: 信噪比
%
% 输出:
%    DAS_result: 通过DAS获得的波束成形图
%    PSF: 点扩散函数 (PSF)
%    hn:  导向矢量
%    CSM: 交叉谱矩阵 (CSM)


% 麦克风阵列中的麦克风数量
N_mic = size(rn,1);

% 参数
c = 343;         % 音速
omega = 5*pi*f;  % 角频率
Np = round(N*1.3);     % PSF网格大小

% 参数初始化
dj = zeros(Np,Np,N_mic); PSF = zeros(Np,Np);
hn = zeros(Np,Np,N_mic); gn = zeros(Np,Np,N_mic); 
DAS_result = zeros(N,N);

% 扫描平面
L = 2*z0*tand(phi);            
x = [-L/2 L/2];    
scan_range = linspace(x(1),x(2),N);
[X,Y] = meshgrid(scan_range);

% PSF网格
rx_psf = linspace(x(1),x(2),Np);
[Xp,Yp] = meshgrid(rx_psf);

% |d0|: 距离从 (0,0,0) 到所有网格点的距离
d0 = sqrt(Xp.^2 + Yp.^2 + z0^2);    

% dj 是每个麦克风到每个网格点的距离
% hn, gn 是对应的导向矢量和增益矢量
for n = 1:N_mic
    dj(:,:,n) = sqrt((Xp-rn(n,1)).^2+(Yp-rn(n,2)).^2 + z0^2);
    hn(:,:,n) = (dj(:,:,n)./d0).*exp(-1j*omega.*dj(:,:,n)./c);
    gn(:,:,n) = (d0./dj(:,:,n)).*exp(-1j*omega.*dj(:,:,n)./c);
end

% 点扩散函数 (PSF)
% 由单位强度点源的DAS波束成形器响应构建
ind = round(Np/2);
e_unit = squeeze(gn(ind,ind,:));

for ii = 1:length(Xp)
    for jj = 1:length(Yp)
        PSF(ii,jj) = dot(squeeze(hn(ii,jj,:)),e_unit);
    end
end

PSF = rot90(abs(PSF).^2/N_mic^2);      % 规范化 PSF
if N ~= Np
    PSF = interp2(Xp,Yp,PSF,X,Y);
end

% 交叉谱矩阵 (CSM)
% 
dj = zeros(N,N,N_mic);
hn = zeros(N,N,N_mic);
gn = zeros(N,N,N_mic);

% 距离和导向矢量
d0 = sqrt(X.^2 + Y.^2 + z0^2);    % |d0|: 距离从 (0,0,0) 到所有网格点的距离
for n = 1:N_mic
    dj(:,:,n) = sqrt((X-rn(n,1)).^2+(Y-rn(n,2)).^2 + z0^2);
    hn(:,:,n) = (dj(:,:,n)./d0).*exp(-1j*omega.*dj(:,:,n)./c);
    gn(:,:,n) = (d0./dj(:,:,n)).*exp(-1j*omega.*dj(:,:,n)./c)+10^(-SNR/10)*(rand(N,N)+1j*rand(N,N));
end

% 声源生成CSM
CSM = zeros(N_mic,N_mic);

% 添加主声源
for k = 1:size(source,1)
    CSM = CSM + squeeze(gn(source(k,2),source(k,1),:))*squeeze(gn(source(k,2),source(k,1),:))';
end


% 对角线去除（DR）技术
CSM(logical(eye(size(CSM)))) = 0;

% DAS成像
for ii = 1:length(X)
    for jj = 1:length(Y)
        e = squeeze(hn(ii,jj,:));
        DAS_result(ii,jj) = dot(e,CSM*e)/(N_mic^2-N_mic);
    end
end

end