核心代碼如下：

%------延時求和波束形成-------
close all;
clear all;
clc;

%%%%%%基陣的有關參數
c=1500;                 %m/s
L=1.3;                    %基陣長
N=32      %基元個數
d=L/(N-1);       %陣元間距
a=9;  %參考基元的在基陣中的位置
xi=[1:1:N]*d;

f=40000;                  %信號頻率
bl=c/f;                 %波長
r0=60;%參考基元到聲源距離
ang0=40;      %改變ang0值可改變束控方向
angle0=ang0*pi/180;                 %角度
fs=f*10;              %采樣頻率
T=0.0005;             %脈沖寬度
w=2*pi*f;
b=30                 %扇面開角范圍-30度至+30度
Ts=1/fs;             %采樣時間間隔
t=[0:Ts:T];          %采樣樣點個數
angle=[-b:0.05:b]*pi/180; %波束個數
R=30;

tao=zeros(N,1); %每個基元的時延
signl=zeros(N,length(t));%每個基元輸出的時間序列
signal=zeros(N,length(t));
for k1=1:N
    tao(k1,1)=(-xi(k1)*sin(angle0))./c;%*cos(angle0)^2
    signl(k1,:)=exp(j*w*t).*exp(j*w*tao(k1,1));%自然指向入射角
    signal(k1,:)=signl(k1,:);
end;

s=zeros(1,length(angle));
y=zeros(length(angle),length(t));
    
for k2=1:length(angle)%轉過波束角
    for k3=1:N
        tao=(-xi(k3)*sin(angle(k2)))./c;
        dfy_signl=signal(k3,:)*exp(-j*w*tao);%用轉過波束做相位修正;*exp(-j*fy(ii))
        y(k2,:)=dfy_signl+y(k2,:);  % 18 路相加在該波束方向的時域輸出
    end
    s(1,k2)=sum(abs(y(k2,:)).^2)/length(t);%在各個波束方向的總輸出
end
figure(1)
plot(angle*180/pi,10*log10(s/max(s)));