bf_phaseshift.m

相移波束形成仿真程序,用24个阵元仿真,形成181个波束,适用于水声信号处理里的波束形成技术

%**************************************************************************
%                          相移波束形成算法
%**************************************************************************
%=============================变量初始化====================================
f0=2e5;       %发射信号频率
fs=10*f0;     %采样频率
c=1500;       %声速
M=24;         %阵元个数
t=5e-4;       %发射脉冲宽度
%==========================================================================
theta=-10*pi/180;    %信号方向
d=0.5*c/f0;          %阵间距
N=fs*t;              %仿真点数
w0=2*pi*f0;          %发射信号角频率
tao=d*sin(theta)/c;  %相邻阵元之间时延
fai=w0*tao;          %相邻阵元之间相位差
%==========================契比雪夫加权系数=================================
A=10.^(30/20);                         %由分贝数改为数值。
n=M-1;
x=0.5*((2*A)^(1/n)+(1/(2*A))^(1/n));   %计算x。
a=1-1/(x^2);                           %计算参数a。
%==========================================================================
w(1)=1;                                 %计算权值。
for k=2:(M/2)
    w(k)=0;
    for s=0:(k-2)
     w(k)=w(k)+(M-1)/(M-k)*nchoosek(k-2,s)*nchoosek(M-k,s+1)*a^(s+1);
    end
end
for i=(M/2+1):M                         %对称得出另一半。
    w(i)=w(M-i+1);
end
%=============================信号的产生===================================
for i=1:M            %造带固定相位差的信号（相对与参考阵元）
    for n=1:N
        X(i,n)=w(i)*exp(w0*n/fs)*exp(j*(i-1)*fai);   %加权
        Y(i,n)=exp(w0*n/fs)*exp(j*(i-1)*fai);        %不加权
    end
end
%**************************************************************************
%                            做常规相移波束形成
%**************************************************************************
beam=-90:1:90;                      %波束指向角
beta=w0*d*sin(beam*pi/180)/c;       %由波束指向角得到的相移
for i=1:181                         %造一组相移信号
    for n=1:M
        aw(i,n)=exp(-j*(n-1)*beta(i));
    end
end
%==========================================================================
S1=aw*X;                             %波束形成(加权)
S2=aw*Y;                             %波束形成(不加权)
for i=1:181
     R1(i)=sum(abs(S1(i,:)))/N;      %检波(加权)
     R2(i)=sum(abs(S2(i,:)))/N;      %检波(不加权)
end
xax=-90:1:90;
plot(xax,10*log(R1/max(R1)),'r:');   %输出波束图(加权)
hold on 
plot(xax,10*log(R2/max(R2)));        %输出波束图(不加权)
grid on;
%==========================================================================