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的指向性');
