浅水声信道模型的建立（5）----自由场，没有散射，3个目标研究

clc;close all;clear all;

%参考论文：声呐目标回波的亮点模型

%自由场下的回波

%发射信号：（z=0,r=10）

%目标：% r1=5000;z1=15;%1号亮点% r2=5010;z2=15;%2号亮点% r3=5050;z3=15;%3号亮点

%接收信号： rr=0;zz=20:30；

%1号目标假设为直径a=4的球体

%参数==================================

c = 1500;            %声速             Unit：m/s

SNR = 60;            %信噪比           Unit：dB

f0=1000;             %信号频率         Unit：Hz

fs=10000;            %采样率          Unit：Hz     

T=5/f0;              %信号脉宽         Unit：s

ts=0:1/fs:T-1/fs;    %一个脉冲包络的持续时间序列

Ns=length(ts);  

startt=0.5;endt=8;   %混响持续时间

w=2*pi*f0;

k=w/c;

%===============================

Sample_time=endt;

Ts = 1/fs;      %采样时间间隔

sample_num = fix(Sample_time*fs);       %采样总点数 FIX：让x向0靠近取整 

nTs = (0:sample_num-1)/fs;              %离散的采样时刻    s

ample_num1 = fix(T*fs);              %信号的采样点数 

nTs1 = (0:sample_num1-1)/fs;            %信号的离散的采样时刻 

receive_signal = [zeros(size(nTs))];      %用于存储水听器接收的信号 

%===============================

%发射信号（z0=0,r0=10）

z0=0;

r0=10;

p=sin(2*pi*f0*ts);       %发射信号

p_fft=fft(p);

%目标1============================================

%目标1位置：r1=5000，z1=15;

%目标1为球形，采用球形的传递函数的幅度反射因子

%求回波

%解决三个参数问题：幅度反射因子Am; 时延taom,相位跳变faim

R1=2;  %目标1的半径为2

r1=5000;z1=15;

distance1=sqrt((r1-r0)^2+(z1-z0)^2);  %声源与目标1等效中心的距离

Am1=0.5*sqrt(R1^2)/sqrt(1+R1/r1)^2;   %球的传递函数幅度反射因子

taom1=(distance1-R1)/c;

faim1=pi;

pb1_fft=exp(1i*k*r1)/r1*Am1*exp(1i*w*taom1)*exp(1i*faim1).*p_fft;   %目标1回波

pb1=ifft(pb1_fft); real_time =(distance1-R1)/c;                                                 %信号的实际到达时刻 3.3333s

signal_start_time = fix(real_time*fs)+1;                                     %信号第一个采样点的时刻

phase = (signal_start_time*Ts-real_time)/Ts*2*pi;                            %得到信号的第一个采样点的相位

receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =pb1;    %模拟采样

pb1 = receive_signal;

%目标1---------------------------------------------------

%散射信号：

sanshe1=conv(p,pb1);

sanshe1_fft=fft(sanshe1);

N=length(sanshe1);

t=startt:1/fs:startt+(N-1)/fs;

%目标1--------------------------------------------------

%阵列接收 rr=0,zz=21:30;

rr=0;

for kkk=1:10    

zz=20+kkk;    

distance1_temp=sqrt((zz-z1)^2+(rr-r1)^2);    %阵元i与目标1的距离    

Am1=0.5*sqrt(R1^2)/sqrt(1+R1/r1)^2;   %球的传递函数幅度反射因子    

taom1=(distance1_temp-R1)/c;    

faim1=pi;    

pb1__fft=exp(1i*k*r1)/r1*Am1*exp(1i*w*taom1)*exp(1i*faim1).*p_fft;   %目标1回波    

pb1__=ifft(pb1__fft);    %1x80049    

%%%%%%%%%%    

real_time =(distance1_temp-R1)/c;                                                    

signal_start_time = fix(real_time*fs)+1;                                        

phase = (signal_start_time*Ts-real_time)/Ts*2*pi;                               

receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =pb1__;     %signal_start_time:33321,sample_num1=50    

pb1_(kkk,:) = receive_signal;  %阵列接收到的信号 end %得到的阵列接收到的信号与散射信号相关。

%这里有一个问题，我在

for kkk=1:10    

total1(kkk,:)=conv(pb1_(kkk,:),sanshe1);      %10x160048    

total1_fft=fft(total1);

end

 %目标1===================================================================

c = 1500;            %声速             Unit：m/s

D=1;                %阵元的距离

M = 10;             % 阵元数

m=[0:M-1];          %阵元序号

d=0.2;           %系统实孔径0.2

T=5/f0;              %信号脉宽         Unit：s %0.0031

k=d/(c*T);

doa1=atan(5/100);                                  %doa为线阵法线方向与声源的夹角

theta = linspace(-90,90,360);                        % 扫描方位角   1X360矩阵 %要构造一个160048的矩阵放drive的数值

drive_=zeros(M,length(total1));

for kkk = 1:length(theta)          %1:360    

tow=sin(theta(kkk)*pi/180)*D/c;          

for kk=1:M                           %M个阵元        

Dtheta = exp(-1i*2*pi*f0*tow*(kk-1)); %补偿    %每一个角度，与每一个阵元形成一个信号补偿        

drive(kk,1) = Dtheta;    %10x1矩阵   %在该扫描角度下，阵列接收的信号补偿            

for temp=1:length(total1)                

drive_(kk,temp)=drive(kk,1);              

end                    

end                  

tempbeamout1(kkk) = sum(abs(sum(total1.*drive_))); %输出   %10x160048矩阵乘以10*160048的矩阵

end

Y1 = 20*log10(tempbeamout1/max(tempbeamout1));

Y1_max=max(Y1);

%目标2===========================================================

%目标2位置：r2=5010，z1=15;

%目标2为球形，采用球形的传递函数的幅度反射因子

%2号目标假设为直径a=8的球体

%求回波

%解决三个参数问题：幅度反射因子Am; 时延taom,相位跳变faim

r2=5010;z2=15;%目标二的位置（圆心所在的位置）

R2=4;         %目标2的半径为4

distance2=sqrt((r2-r0)^2+(z2-z0)^2);

Am2=0.5*sqrt(R2^2)/sqrt(1+R2/r2)^2;   %球的传递函数幅度反射因子

taom2=(distance2-R2)/c;

faim2=pi;

pb2_fft=exp(1i*k*r2)/r2*Am2*exp(1i*w*taom2)*exp(1i*faim2).*p_fft;  %目标2回波

pb2=ifft(pb2_fft); real_time = (distance2-R2)/c;                                                            %信号的实际到达时刻 3.3333s

signal_start_time = fix(real_time*fs)+1;                                     %信号第一个采样点的时刻

phase = (signal_start_time*Ts-real_time)/Ts*2*pi;                            %得到信号的第一个采样点的相位

receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =pb2; %模拟采样

pb2 = receive_signal ;

%目标2------------------------------------------------

sanshe2=conv(p,pb2);

sanshe2_fft=fft(sanshe2);

%目标2------------------------------------------------

%阵列接收 rr=0,zz=21:30;

rr=0;

for kkk=1:10    

zz=20+kkk;    

distance2_temp=sqrt((zz-z2)^2+(rr-r2)^2);    %阵元i与目标1的距离    

Am1=0.5*sqrt(R2^2)/sqrt(1+R2/r2)^2;   %球的传递函数幅度反射因子    

taom2=(distance2_temp-R2)/c;    

faim1=pi;    

pb2__fft=exp(1i*k*r2)/r2*Am2*exp(1i*w*taom2)*exp(1i*faim2).*p_fft;   %目标1回波    

pb2__=ifft(pb2__fft);    %1x80049    

%%%%%%%%%%    

real_time =(distance1_temp-R2)/c;                                                    

signal_start_time = fix(real_time*fs)+1;                                        

phase = (signal_start_time*Ts-real_time)/Ts*2*pi;                               

receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =pb2__;     %signal_start_time:33321,sample_num1=50    

pb2_(kkk,:) = receive_signal;  %阵列接收到的信号

end

%得到的阵列接收到的信号与散射信号相关。

for kkk=1:10    

total2(kkk,:)=conv(pb2_(kkk,:),sanshe2);      %10x160048    

total2_fft=fft(total2);

end

 %目标2===================================================================

c = 1500;            %声速             Unit：m/s

D=1;                %阵元的距离

M = 10;             % 阵元数

m=[0:M-1];          %阵元序号

d=0.2;           %系统实孔径0.2

T=5/f0;              %信号脉宽         Unit：s %0.0031

k=d/(c*T);

doa2=atan(5/5010);                                  %doa为线阵法线方向与声源的夹角

theta = linspace(-90,90,360);                        % 扫描方位角   1X360矩阵 %要构造一个160048的矩阵放drive的数值

drive_=zeros(M,length(total2));

for kkk = 1:length(theta)          %1:360    

tow=sin(theta(kkk)*pi/180)*D/c;          

for kk=1:M                           %M个阵元        

Dtheta = exp(-1i*2*pi*f0*tow*(kk-1)); %补偿    %每一个角度，与每一个阵元形成一个信号补偿        

drive(kk,1) = Dtheta;    %10x1矩阵   %在该扫描角度下，阵列接收的信号补偿            

for temp=1:length(total2)                

drive_(kk,temp)=drive(kk,1);              

end                    

end                  

tempbeamout2(kkk) = sum(abs(sum(total2.*drive_))); %输出   %10x160048矩阵乘以10*160048的矩阵

end

Y2 = 20*log10(tempbeamout2/max(tempbeamout2));

Y2_max=max(Y2);

%目标3==============================================

%目标3位置：r3=5050，z3=15;

%目标3为球形，采用球形的传递函数的幅度反射因子

%3号目标假设为直径a=4的球体 %求回波

%解决三个参数问题：幅度反射因子Am; 时延taom,相位跳变faim

r3=5050;z3=15;         %目标3的位置（圆心所在的位置）

R3=2;         %目标3的半径为2

distance3=sqrt((r3-r0)^2+(z3-z0)^2);

Am3=0.5*sqrt(R3^2)/sqrt(1+R3/r3)^2;   %球的传递函数幅度反射因子

taom3=(distance3-R3)/c;

faim3=pi;

pb3_fft=exp(1i*k*r3)/r3*Am3*exp(1i*w*taom3)*exp(1i*faim3).*p_fft;  %目标2回波

pb3=ifft(pb3_fft); real_time = (distance3-R3)/c;                                                            %信号的实际到达时刻 3.3333s

signal_start_time = fix(real_time*fs)+1;                                     %信号第一个采样点的时刻

phase = (signal_start_time*Ts-real_time)/Ts*2*pi;                            %得到信号的第一个采样点的相位

receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =pb3; %模拟采样

pb3= receive_signal ;

%目标3------------------------------------------------

sanshe3=conv(p,pb3); sanshe3_fft=fft(sanshe3);

%目标3------------------------------------------------------------------------

%阵列接收 rr=0,zz=21:30;

rr=0;

for kkk=1:10    

zz=20+kkk;    

distance3_temp=sqrt((zz-z3)^2+(rr-r3)^2);    %阵元i与目标1的距离    

Am3=0.5*sqrt(R3^2)/sqrt(1+R3/r3)^2;   %球的传递函数幅度反射因子    

taom3=(distance3_temp-R3)/c;    

faim3=pi;    

pb3__fft=exp(1i*k*r3)/r3*Am3*exp(1i*w*taom3)*exp(1i*faim3).*p_fft;   %目标1回波    

pb3__=ifft(pb3__fft);    %1x80049    

%%%%%%%%%%    

real_time =(distance3_temp-R3)/c;                                                    

signal_start_time = fix(real_time*fs)+1;                                        

phase = (signal_start_time*Ts-real_time)/Ts*2*pi;                               

receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =pb3__;     %signal_start_time:33321,sample_num1=50    

pb3_(kkk,:) = receive_signal;  %阵列接收到的信号

end

%得到的阵列接收到的信号与散射信号相关。

for kkk=1:10   

  total3(kkk,:)=conv(pb3_(kkk,:),sanshe3);      %10x160048    

total3_fft=fft(total3);

end

 %目标3===================================================================

c = 1500;            %声速             Unit：m/s

D=1;                %阵元的距离

M = 10;             % 阵元数

m=[0:M-1];          %阵元序号

d=0.2;           %系统实孔径0.2

T=5/f0;              %信号脉宽         Unit：s %0.0031

k=d/(c*T);

doa3=atan(5/5050);                                  %doa为线阵法线方向与声源的夹角

theta = linspace(-90,90,360);                        % 扫描方位角   1X360矩阵 %要构造一个160048的矩阵放drive的数值

drive_=zeros(M,length(total1));

for kkk = 1:length(theta)          %1:360    

tow=sin(theta(kkk)*pi/180)*D/c;          

for kk=1:M                           %M个阵元       

  Dtheta = exp(-1i*2*pi*f0*tow*(kk-1)); %补偿    %每一个角度，与每一个阵元形成一个信号补偿       

  drive(kk,1) = Dtheta;    %10x1矩阵   %在该扫描角度下，阵列接收的信号补偿            

for temp=1:length(total1)                

drive_(kk,temp)=drive(kk,1);             

  end                    

end                 

  tempbeamout3(kkk) = sum(abs(sum(total3.*drive_))); %输出   %10x160048矩阵乘以10*160048的矩阵

end

Y3 = 20*log10(tempbeamout3/max(tempbeamout3));

Y3_max=max(Y3);

% %===========================================================================

%画图：

figure(1);

subplot(321);plot(1:length(pb1_fft),abs(pb1_fft)); title('目标1--回波信号--频域');xlabel('omega/pi');ylabel('幅值');

subplot(322);plot(1:length(pb1),pb1); title('目标1--回波信号--时域');xlabel('时间');ylabel('幅值');

subplot(323);plot(1:length(pb2_fft),abs(pb2_fft)); title('目标2---回波信号-频域');xlabel('omega/pi');ylabel('幅值');

subplot(324);plot(1:length(pb2),pb2); title('目标2---回波信号--时域');xlabel('时间');ylabel('幅值');

subplot(325);plot(1:length(pb3_fft),abs(pb3_fft)); title('目标3---回波信号-频域');xlabel('omega/pi');ylabel('幅值');

subplot(326);plot(1:length(pb3),pb3); title('目标3---回波信号--时域');xlabel('时间');ylabel('幅值');

% %===========================================================================

figure(2);

subplot(321);plot(t,sanshe1);xlabel('时间');ylabel('幅值');title('目标1--散射回波--时域');

subplot(322);plot(1:length(sanshe1_fft),sanshe1_fft);xlabel('omega/pi');ylabel('|e^j^omega|');title('目标1--散射回波--频域');

subplot(323);plot(t,sanshe2);xlabel('时间');ylabel('幅值');title('目标2--散射回波--时域');

subplot(324);plot(1:length(sanshe2_fft),sanshe2_fft);xlabel('omega/pi');ylabel('|e^j^omega|');title('目标2--散射回波--频域');

subplot(325);plot(t,sanshe3);xlabel('时间');ylabel('幅值');title('目标3--散射回波--时域');

subplot(326);plot(1:length(sanshe3_fft),sanshe3_fft);xlabel('omega/pi');ylabel('|e^j^omega|');title('目标3--散射回波--频域');

%%============================================================================

figure(3);

subplot(321);plot(1:length(pb1(1,:)),pb1(1,:)); xlabel('omega/pi');ylabel('|e^j^omega|');title('阵元1接收到的信号--时域');

subplot(322);plot(1:length(pb1_fft(1,:)),abs(pb1_fft(1,:))); xlabel('omega/pi');ylabel('|e^j^omega|');title('阵元1接收到的信号--频域');

subplot(323);plot(1:length(pb2(1,:)),pb2(1,:)); xlabel('omega/pi');ylabel('|e^j^omega|');title('阵元2接收到的信号--时域');

subplot(324);plot(1:length(pb2_fft(1,:)),abs(pb2_fft(1,:))); xlabel('omega/pi');ylabel('|e^j^omega|');title('阵元2接收到的信号--频域');

subplot(325);plot(1:length(pb3(1,:)),pb3(1,:)); xlabel('omega/pi');ylabel('|e^j^omega|');title('阵元3接收到的信号--时域');

subplot(326);plot(1:length(pb3_fft(1,:)),abs(pb3_fft(1,:))); xlabel('omega/pi');ylabel('|e^j^omega|');title('阵元3接收到的信号--频域');

%%============================================================================

figure(4);

subplot(321);plot(1:length(total1(1,:)),total1(1,:)); xlabel('时间');ylabel('幅值');title('阵元1接收到的信号与散射信号相关--total1--时域');

subplot(322);plot(1:length(total1_fft(1,:)),abs(total1_fft(1,:))); xlabel('omega/pi');ylabel('|e^j^omega|');title('阵元1接收到的信号与散射信号相关--total1_fft--频域');

subplot(323);plot(1:length(total2(1,:)),total2(1,:)); xlabel('时间');ylabel('幅值');title('阵元2接收到的信号与散射信号相关--total2--时域');

subplot(324);plot(1:length(total2_fft(1,:)),abs(total2_fft(1,:))); xlabel('omega/pi');ylabel('|e^j^omega|');title('阵元2接收到的信号与散射信号相关--total2_fft--频域');

subplot(325);plot(1:length(total3(1,:)),total3(1,:)); xlabel('时间');ylabel('幅值');title('阵元3接收到的信号与散射信号相关--total3--时域');

subplot(326);plot(1:length(total3_fft(1,:)),abs(total3_fft(1,:))); xlabel('omega/pi');ylabel('|e^j^omega|');title('阵元3接收到的信号与散射信号相关--total3_fft--频域');

%%===========================================================================

figure(5);

subplot(311);plot(theta,Y1);title('目标1的指向性');

subplot(312);plot(theta,Y2);title('目标2的指向性');

subplot(313);plot(theta,Y3);title('目标3的指向性');