hblxjm.m

目标连续声纳回波模型仿真 应用于水下目标回波分析与仿真

close all
clear all
clc

f0=35000;
w=2*pi*f0;
fs=192000;
r=0.191;
d=0.225;
l=7.500;
c=1480.000;
%-------------------------------------
N=20;
M=20;
%=====================================
fan=2*pi/M;                         %圆角划分
h=0.5*cos(22.5*pi/180)/N;   
S1=2*pi*r*(2*r/N)/M;
S2=2*pi*r*0.735/N/M;
for i=1:N;
    for j=1:M
        S3(i,j)=pi*r*(2*(N-i)+1)*0.5/N/N/M;
    end
end
S4=2*r^2/N;
%--------------------------
tt=1/fs:1/fs:6/f0;            %单频信号
s=1000*sin(2*pi*f0*tt); 
%--------------------------
n=length(s);
k=10;                %数据长度设定！！！    
m=(2*k+1)*n;
mn=5;                %取包络
t=1/fs:1/fs:m/fs;
t=t*1000;
x=[zeros(1,k*n),s(1:n),zeros(1,k*n)];
%--------------------------------------------------------------
%构造影射距阵
%----------------------------- jiao ---------------------------
for ij=1:4 
%     cta=0;                                                  %cta=0度
    cta=pi*ij/4;
for i=1:N
    for j=1:M
        jiao1(i,j)=(r-i*2*r/N+r/N)/r;                       %1，球面
        jiao2(i,j)=sin(cta)*cos(j*fan-pi/2-fan/2);          %2，柱面
        jjiao(i)=atan((N-i+0.5)*r/N/(i-0.5)*h);             %3，锥面
        jjiao(i)=abs(pi-jjiao(i)-cta);
        if jjiao(i)>=pi/2 
            jiao3(i,j)=0;
        elseif  jjiao(i)<pi/2
            jiao3(i,j)=cos(jjiao(i))*cos(j*fan-pi/2-fan/2);
        end
 %------------------------------ r --------------------------------
        if cta<=pi/2
            rr1=-r-d*cos(cta)+i*2*r/N-r/N;
        else
            rr1=-(r/sin(cta)+d*cos(cta)-r*cos(cta)/tan(cta))+i*2*r/N-r/N;
        end
        r1(i,j)=l+rr1;                         %球面正对入射距离
        tao1(i,j)=2*rr1/c;
        rr2=-(r/sin(cta)+(d-i*0.735/N+0.735/2/N)*cos(cta)-r*cos(cta)/tan(cta))+r*sin(cta)*(1-cos(j*fan-pi/2-fan/2));
        r2(i,j)=l+rr2;                         %柱面正对入射距离
        tao2(i,j)=2*rr2/c;
        rl=0.5-i*0.5/N+0.25/N;                 %锥面正对入射距离
        rr3=0.940-rl*cos(22.5*pi/180)-rl*sin(22.5*pi/180)*tan(cta);
        rr3=rr3*cos(cta)+rl*sin(22.5*pi/180)*sin(cta)*(1-cos(j*fan-pi/2-fan/2));
        r3(i,j)=l+rr3;            
        tao3(i,j)=2*rr3/c;
        rr4=(0.510+i*2*r/N)*cos(cta);          %平面正对入射距离
        r4(i)=l+rr4;            
        tao4(i)=2*rr4/c;
    end
end
xx=[x,x,x];
X=zeros(1,m);
%--------------------------------------------------------------
%cta=0,球面回波仿真
%--------------------------------------------------------------
A1=6*S1*jiao1./r1.^2/2/pi;                  
m1=n*(2*k+1)-tao1*fs;
mm1=m1+m-1;
if cta==0                                       %cta=0,球面回波仿真
    for i=1:N/2 
       for j=1:M
           X=A1(i,j)*xx(m1(i,j):mm1(i,j))+X;
       end
    end
elseif cta==pi/4                                %cta=pi/4,球面回波仿真
      for i=1:N/4  
         for j=1:M
             X=A1(i,j)*xx(m1(i,j):mm1(i,j))+X;
         end
      end
      for i=N/4+1:N/2                                
         for j=1:M/2
             X=A1(i,j)*xx(m1(i,j):mm1(i,j))+X;
         end
      end
elseif cta==pi/2                                %cta=pi/2,球面回波仿真
      for i=1:N/2  
         for j=1:M/2
             X=A1(i,j)*xx(m1(i,j):mm1(i,j))+X;
         end
      end
elseif cta==pi*3/4                              %cta=pi*3/4,球面回波仿真
      for i=1:N/4                                
         for j=1:M/2
             X=A1(i+N/4,j)*xx(m1(i,j):mm1(i,j))+X;
         end
      end
end
%-------------------------------------------------------------
% figure(1)                                         %cta=0度
% subplot(421),plot(t,X)
% % axis ([0 0.00319 -2e4 2e4]);
% title('0度入射时仿真回波信号')
% xlabel('时间/ms')
% ylabel('幅度/mV')
%--------------------------------------------------------------
%柱面回波仿真
%--------------------------------------------------------------
A2=6*S2*jiao2./r2.^2/2/pi;
m2=n*(2*k+1)-tao2*fs;
mm2=m2+m-1;
for i=1:N                                    %cta=0,柱面回波仿真
    for j=1:M/2
         X=A2(i,j)*xx(m2(i,j):mm2(i,j))+X;
    end
end
% figure(2)
% subplot(421),plot(t,X)
% % axis ([0 0.00319 -2e4 2e4]);
% title('45度入射时仿真回波信号')
% xlabel('时间/ms')
% ylabel('幅度/mV')
%--------------------------------------------------------------
%锥面回波仿真
%--------------------------------------------------------------
A3=6*S3.*jiao3./r3.^2/2/pi;                       %cta=0,锥面回波仿真
m3=n*(2*k+1)-tao3*fs;
mm3=m3+m-1;
if cta<=pi*112.5/180 
   for i=1:N                                    %0<cta & cta<pi*112.5/180,锥面回波仿真
      for j=1:M/2
         X=A3(i,j)*xx(m3(i,j):mm3(i,j))+X;
      end
   end
else
    for i=1:N                                    %pi*112.5/180<cta & cta<pi,锥面回波仿真
      for j=1:M
         X=abs(A3(i,j))*xx(m3(i,j):mm3(i,j))+X;
      end
   end
end
% figure(3)
% subplot(421),plot(t,X)
% % axis ([0 0.00319 -2e4 2e4]);
% title('90度入射时仿真回波信号')
% xlabel('时间/ms')
% ylabel('幅度/mV')
%--------------------------------------------------------------
%平面回波仿真
%--------------------------------------------------------------
A4=6*S4*sin(cta)./r4.^2/2/pi;                       %cta=0,平面回波仿真
m4=n*(2*k+1)-tao4*fs;
mm4=m4+m-1;
for i=1:N                                    %cta=0,柱面回波仿真
   X=A4(i)*xx(m4(i):mm4(i))+X;
end
xg(i,:)=xcorr(X,s);
xg(i,:)=abs(xg(i,:));
xg(i,:)=xg(i,:)/max(xg(i,:));
%+++++++++++++++++++++++++++++++++++++++++++++++
for iiii=1:(length(xg(i,:))-1)/mn
    XG(i,iiii)=max(xg(i,((iiii-1)*mn+1):iiii*mn));
end
figure(ij)
subplot(421),plot(t,X)
% axis ([0 0.00319 -2e4 2e4]);
title('45度入射时仿真回波信号')
xlabel('时间/ms')
ylabel('幅度/mV')
end