基于MATLAB的波速形成仿真

Posted 2022-01-03 fpga&matlab

close all
clear all

%%波束形成
f=[4000,5000,6000];%信号频率，单位Hz
w=[1,1,1,1,1;%每一行表示一种加权
    1,2,4,2,1;
    4,2,1,2,4];
fs=4*f;%采样频率取信号频率的4倍
T=0.1;%信号长度,单位s
N=T*fs;%采样点数
c=1500;%声速,m/s
d=0.15;%阵元间距，单位m
jay=sqrt(-1);%复单位
theta=-90:1:90;%观测角度范围，单位度
tao=d*sin(pi*theta/180)/c;%相邻阵元时延
xx=zeros(3,2400);%准备产生三种频率下的参考信号

for i=1:3
    xx(i,1:N(i))=exp(jay*(2*pi*f(i)/fs(i)*(0:N(i)-1)));%三种频率下的参考信号
end

for j=1:3%三种频率循环
    for i=1:length(tao)%不同角度循环
        for k=(j-1)*5+1:(j-1)*5+5%第一种频率下的五个阵元的信号
            x(k,:)=xx(j,:)*exp(-jay*2*pi*f(j)*(k-1)*tao(i));
        end
        r=x((j-1)*5+1:(j-1)*5+5,:)*x((j-1)*5+1:(j-1)*5+5,:)'/N(j);%求某一角度下的波束图，用公式B=w*(∑x*x')*w'/N.
        B((j-1)*3+1,i)=w(1,:)*r*w(1,:)';%第j个频率下的第一种加权的波束图输出
        B((j-1)*3+2,i)=w(2,:)*r*w(2,:)';%第j个频率下的第二种加权的波束图输出
        B((j-1)*3+3,i)=w(3,:)*r*w(3,:)';%第j个频率下的第三种加权的波束图输出
    end
end
%c产生的举证B表示B=[B11;B12;B13;B21;B22;B23;B31;B32;B33];第一个下角标表示第几个频率，第二个下角标表示第几种加权而得到的波束
for k=1:9%分别对9个波束图进行归一化和取对数
    B(k,:)=real(B(k,:))/max(real(B(k,:)));%归一化
    B(k,:)=20*log10(B(k,:));%取对数
end
figure(1);%频率f=4000，加权不同
subplot(3,1,1);
plot(theta,B(1,:));axis([-90 90 -60 0]);title('频率f=4000，加权[1,1,1,1,1]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,2);
plot(theta,B(2,:));axis([-90 90 -60 0]);title('频率f=4000，加权[1,2,4,2,1]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,3);
plot(theta,B(3,:));axis([-90 90 -60 0]);title('频率f=4000，加权[4,2,1,2,4]');xlabel('角度/°');ylabel('波束图/dB');

figure(2);%频率f=5000，加权不同
subplot(3,1,1);
plot(theta,B(4,:));axis([-90 90 -60 0]);title('频率f=5000，加权[1,1,1,1,1]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,2);
plot(theta,B(5,:));axis([-90 90 -60 0]);title('频率f=5000，加权[1,2,4,2,1]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,3);
plot(theta,B(6,:));axis([-90 90 -60 0]);title('频率f=5000，加权[4,2,1,2,4]');xlabel('角度/°');ylabel('波束图/dB');

figure(3);%频率f=6000，加权不同
subplot(3,1,1);
plot(theta,B(7,:));axis([-90 90 -60 0]);title('频率f=6000，加权[1,1,1,1,1]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,2);
plot(theta,B(8,:));axis([-90 90 -60 0]);title('频率f=6000，加权[1,2,4,2,1]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,3);
plot(theta,B(9,:));axis([-90 90 -60 0]);title('频率f=6000，加权[4,2,1,2,4]');xlabel('角度/°');ylabel('波束图/dB');

figure(4);%加权[1,1,1,1,1],频率不同
subplot(3,1,1);
plot(theta,B(1,:));axis([-90 90 -60 0]);title('频率f=4000，加权[1,1,1,1,1]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,2);
plot(theta,B(4,:));axis([-90 90 -60 0]);title('频率f=5000，加权[1,1,1,1,1]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,3);
plot(theta,B(7,:));axis([-90 90 -60 0]);title('频率f=6000，加权[1,1,1,1,1]');xlabel('角度/°');ylabel('波束图/dB');

figure(5);%加权[1,2,4,2,1],频率不同
subplot(3,1,1);
plot(theta,B(2,:));axis([-90 90 -60 0]);title('频率f=4000，加权[1,2,4,2,1]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,2);
plot(theta,B(5,:));axis([-90 90 -60 0]);title('频率f=5000，加权[1,2,4,2,1]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,3);
plot(theta,B(8,:));axis([-90 90 -60 0]);title('频率f=6000，加权[1,2,4,2,1]');xlabel('角度/°');ylabel('波束图/dB');

figure(6);%加权[4,2,1,2,4],频率不同
subplot(3,1,1);
plot(theta,B(3,:));axis([-90 90 -60 0]);title('频率f=4000，加权[4,2,1,2,4]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,2);
plot(theta,B(6,:));axis([-90 90 -60 0]);title('频率f=5000，加权[4,2,1,2,4]');xlabel('角度/°');ylabel('波束图/dB');
subplot(3,1,3);
plot(theta,B(9,:));axis([-90 90 -60 0]);title('频率f=6000，加权[4,2,1,2,4]');xlabel('角度/°');ylabel('波束图/dB');

%计算主瓣宽度D和旁瓣级P
for k=1:9
    [i,j]=find(abs(B(k,:)+3)<0.5);%误差限设在0.5dB。
    D(k)=(j(length(j))-j(1))*1;%每一个点之间是一度
end    %根据计算D的结果可以发现，同一种频率下三种加权分辨率最高的是第三种，其次是第一种。同一加权下，频率越高，分辨率越高
disp('9个波束图的参数依次为，其中分别表示:');
disp('f=4000，w=[1,1,1,1,1];f=4000，w=[1,2,4,2,1];f=4000，w=[4,2,1,2,4]');
disp('f=5000，w=[1,1,1,1,1];f=5000，w=[1,2,4,2,1];f=5000，w=[4,2,1,2,4]');
disp('f=6000，w=[1,1,1,1,1];f=6000，w=[1,2,4,2,1];f=6000，w=[4,2,1,2,4]');
disp('波束宽度为：');
ss=sprintf('%d ',D);
disp(ss);
NN=length(theta);%每一个波束图离散的点数
for k=1:9
    for i=(NN-1)/2+1:NN-1%从主波束开始往后找到第一个次级大的位置
        if B(k,i)>=B(k,i+1)&&B(k,i)>=B(k,i-1)
            CC(k)=B(k,i);%记录第一个次级大的分贝数
        end
    end
end
CC(2)=B(2,1);CC(5)=B(5,1);%由于在[-90,90]度的范围内，第二种和第五种情况下没有找到次级大，故用它在-90°的地方分贝数近似。
disp('旁絆级为：');
s=sprintf('%f ',CC);
disp(s);