我貌似前面写的两篇博客里面的理论都是错误的。。。
程序说明：
发射信号点为（r=0,z=10）,
亮点1位置为（r1=5000,z1=15）;
接收信号位置为：（rr=0,zz=20:30）;
条件：
声源发射信号，假设海面与最大深度处，信号只经历一次反射
整个过程分为三个过程：
- 声源到亮点1，得到信号 signal_11（时域）
- 在亮点1处散射得到 P1
R1=conv(P1,signal_11), R1为散射信号。
3.散射信号返回到接收信号处，从z=20到z=30,共10个阵元。假设海面和海底还是经历一次反射，10个阵元接收到10个信号，得到signal_13 (这是一个10行的矩阵，
每行代表一个阵元)；
Signal_13的每行数据分别和 散射信号R1进行相关，得到一个10行的矩阵
Zhenlie_1(1:10,:)=conv(signal_13,R1)
从-90度扫描到90度，进行时域波束形成，得到指向性图，但是可能是距离过远，得到的指向性图不够明显，所以进行加权，得到加权指向性图。
%==============================================================================================
首先编写function函数，简化程序的篇幅：
function receive_signal= Com_Sonar_data_generation_2(c,fs,f0,T,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num)
%**********************************通信声呐水听器数据生成函数************************************
%几点假设：
%（1）基于射线声学理论
%（2）几何衰减按球面波传播衰减规律衰减，不考虑吸收衰减
%（3）考虑水面和水底的反射（浅海声信道？）
%（4）考虑在高斯白噪声背景下
%（5）整个空间声速分布均匀
%***************************************输入参数说明******************************************
% c 声速 Unit：m/s
% fs 采样率 Unit：Hz
% f0 信号频率 Unit：Hz
% T 信号脉宽 Unit：s
% Sample_time 采样时长 Unit：s 假设信号发射的时刻为零时刻
% SNR 信噪比 Unit：dB
% H 水深 Unit：m
% H1 发射换能器水深 Unit：m
% H2 接收换能器水深 Unit：m
% D 接收与发射换能器水平距离 Unit：m
% Re_coef_surf 水面反射系数
% Re_coef_bottom 水底反射系数 存在半波损失
% Reflex_num 考虑最大的反射次数
%zhenfu 发射信号的振幅
%phasee 发射信号的相位
%***************************************输出参数说明******************************************
% receive_signal 以信号发射为起始采样点的水听器接收信号
% receive_signal=[];
%==================================================
H1b = H - H1; %发射换能器到水底的距离 Unit：m
H2b = H - H2; %接收换能器到水底的距离 Unit：m
Ts = 1/fs; %采样时间间隔
sample_num = fix(Sample_time*fs); %采样总点数 FIX：让x向0靠近取整
nTs = (0:sample_num-1)/fs; %离散的采样时刻
sample_num1 = fix(T*fs); %信号的采样点数
nTs1 = (0:sample_num1-1)/fs; %信号的离散的采样时刻
receive_signal = zeros(size(nTs)); %用于存储水听器接收的信号
d0 = sqrt((H1-H2)^2+D^2); %两个换能器的直线距离
S0 = 1/d0; %直达波的声压幅值
Noise_var = 10^(-SNR/20)*S0; %白噪声的方差
%==========================================================
%% 计算到达回波信号
real_time = d0/c; %信号的实际到达时刻 3.3354
signal_start_time = fix(real_time*fs)+1; %信号第一个采样点的时刻 53366(应该是说第53366个采样点开始才是到达的时刻的数据)
phase = (signal_start_time*Ts-real_time)/Ts*2*pi; %得到信号的第一个采样点的相位
receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) = S0*sin(2*pi*f0*nTs1+phase); %模拟采样
receive_signal = receive_signal + Noise_var*randn(size(receive_signal)); %加入噪声
if Reflex_num>0 %考虑反射时
for i = 1:Reflex_num %考虑i次反射的情况
if mod(i,2) == 0 %偶数次时
%---------------------------------------------------第一次从水面反射时---------------------------------------------------
M = i*H+H1-H2; %中间系数
k1 = sqrt(1+(D/M)^2); %中间系数
d1 = k1*H1; %发射换能器距离水面第一个反射点的距离
d2 = k1*H; %界面相邻两点间的距离
d3 = k1*(H-H2); %水听器距离边界最后一个反射点最近的距离
di = d1+d2*(i-1)+d3; %反射波总路程
Si = 1/di*Re_coef_surf^(i/2)*Re_coef_bottom^(i/2); %反射波的声压幅值
real_time = di/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)) =...
receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) + Si*sin(2*pi*f0*nTs1 + phase); %模拟采样
%---------------------------------------------------第一次从水底反射时---------------------------------------------------
Mb = i*H+H1b-H2b; %中间系数
k1b = sqrt(1+(D/Mb)^2); %中间系数
d1b = k1b*H1b; %发射换能器距离水底第一个反射点的距离
d2b = k1b*H; %界面相邻两点间的距离
d3b = k1b*(H-H2b); %水听器距离边界最后一个反射点最近的距离
dib = d1b+d2b*(i-1)+d3b; %反射波总路程
Sib = 1/dib*Re_coef_surf^(i/2)*Re_coef_bottom^(i/2); %反射波的声压幅值
real_time = dib/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)) =...
receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) + Sib*sin(2*pi*f0*nTs1 + phase); %模拟采样
else %奇数次时
%---------------------------------------------------第一次从水面反射时---------------------------------------------------
M = (i-1)*H+H1+H2; %中间系数
k1 = sqrt(1+(D/M)^2); %中间系数
d1 = k1*H1; %发射换能器距离水面第一个反射点的距离
d2 = k1*H; %界面相邻两点间的距离
d3 = k1*H2; %水听器距离边界最后一个反射点最近的距离
di = d1+d2*(i-1)+d3; %反射波总路程
Si = 1/di*Re_coef_surf^((i-1)/2+1)*Re_coef_bottom^((i-1)/2); %反射波的声压幅值
real_time = di/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)) =...
receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) + Si*sin(2*pi*f0*nTs1 + phase); %模拟采样
%---------------------------------------------------第一次从水底反射时---------------------------------------------------
Mb = (i-1)*H+H1b+H2b; %中间系数
k1b = sqrt(1+(D/Mb)^2); %中间系数
d1b = k1b*H1b; %发射换能器距离水底第一个反射点的距离
d2b = k1b*H; %界面相邻两点间的距离
d3b = k1b*H2b; %水听器距离边界最后一个反射点最近的距离
dib = d1b+d2b*(i-1)+d3b; %反射波总路程
Sib = 1/dib*Re_coef_surf^((i-1)/2)*Re_coef_bottom^((i-1)/2+1); %反射波的声压幅值
real_time = dib/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))=receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) + Sib*sin(2*pi*f0*nTs1 + phase); %模拟采样
end
end
end
***************************************************快乐的分割线*****************************************************************************************
主程序============================================================================================
clc;close all;clear all;
% r=0;z=10;%发射信号的位置。 % r1=5000;z1=15;%1号亮点 % r2=5010;z2=15;%2号亮点 % r3=5050;z3=15;%3号亮点 % rr=0;zz=20:30;%接收信号
c = 1500; %声速 Unit：m/s
SNR = 60; %信噪比 Unit：dB
H = 100; %水深 Unit：m
%Sample_time = 0.1; %采样时长 Unit：s 假设信号发射的时刻为零时刻
Re_coef_surf = -1; %水面反射系数
Re_coef_bottom = 0.8; %水底反射系数
Reflex_num =1; %考虑最大的反射次数
%=====================================
f0=1600; %信号频率 Unit：Hz
fs=16000; %采样率 Unit：Hz
T=5/f0; %信号脉宽 Unit：s %0.0031
ts=0:1/fs:T-1/fs; %一个脉冲包络的持续时间序列
Ns=length(ts); %50
startt=0.5;endt=3.5;%混响持续时间
Sample_time=endt;
azm=pi/6;
%============================亮点1===========================================
%第一过程:发射换能器到亮点1 % r=0;z=10;%发射信号的位置。
H1 = 10; %发射点水深 Unit：m
H2 = 15; %接收点水深 Unit：m
D = 5000; %接收与发射水平距离 Unit：m
%****************************************************************************
signal_11=Com_Sonar_data_generation_2(c,fs,f0,T,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num); signal_11_fft=fft(signal_11); %56000
X1=abs(fft(signal_11));
X1=X1/max(X1);
figure(1);
plot(signal_11);
xlabel('时间s');ylabel('幅度'); string = ['信号发生器—亮点1，信噪比SNR=',num2str(SNR),'dB,考虑的最大反射次数',num2str(Reflex_num),'时域'];title(string);
figure(2);
plot(1:length(signal_11_fft),signal_11_fft);
xlabel('频率omega');ylabel('幅度'); string = ['信号发生器—亮点1，信噪比SNR=',num2str(SNR),'dB,考虑的最大反射次数',num2str(Reflex_num),'频域'];title(string);
figure(3);
N=length(signal_11); f=1:N/2; plot(f,signal_11_fft(1:N/2)); xlabel('频率/Hz');ylabel('归一化幅值');title('频谱的一半');grid on;
%==========================================================
%过程2：从亮点散射
Np=(endt-startt)*fs; %散射系数的点数48000个点，发射信号持续时间内采样的点
for th=startt:1/fs:endt-1/fs %混响时间 Np48000个点,th共循环48000次
r1=c*(th-T)/2; %T=0.005
r2=c*th/2; %c=1500
N=round((r2^2-r1^2)*azm/2); %azm=0.5236（pi/6）方位角 %随着混响时间的增加，混响点数也在增加
ain=randn(1,N); %幅度
fai=randn(1,N)*2*pi; %相位
for i=1:N
p1(i)=ain(i)*exp(1i*fai(i));
end
P1(round((th-startt)*fs)+1)=sum(p1)/r1^2; %48000
end
% 得到混响信号
R1=conv(signal_11,P1); % length(signal_11)+length(P1)-1; 56000+48000-1=103999
Rr1=real(R1);
figure(4);
t=startt:1/fs:endt+1240*T-6002/fs; %: 共103999个点,这个是凑出来的
plot(t,Rr1/max(Rr1));
axis([startt endt+1240*T-6002/fs -1 1]);xlabel('时间/s');ylabel('归一化幅值');title('亮点1处--混响');
%===================================================================
figure(5);
R1_fft=fft(R1);
plot(1:length(R1_fft),R1_fft);
xlabel('omega/pi');ylabel('归一化幅值|e^j^omega|');title('亮点1处--混响---相关后傅里叶');
%=====================================================================
figure(6);
R1_fft=fft(R1);
plot(1:length(R1_fft)/2,R1_fft(1:length(R1_fft)/2));
xlabel('omega/pi');ylabel('归一化幅值|e^j^omega|');title('亮点1处--混响---相关后傅里叶--频谱的一半');
%=================================================================================
%重新思考第三过程： %选取阵列第一个点为参考点，zz=20。zz=20:30,rr=0; %亮点1散射后到达阵列1 r1=5000;z1=15;%1号亮点 rr=0;zz=20:30;%接收信号
H1 = 15; %发射点水深 Unit：m
D = 5000; %接收与发射水平距离 Unit：m
for kkk=1:10; %接收点水深 Unit：m
H2=20+1*kkk;
signal_13(kkk,:)=Com_Sonar_data_generation_2(c,fs,f0,T,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);
end
for kkk=1:10
zhenlie_1(kkk,:)=conv(signal_13(kkk,:),R1);
end
total_1=zhenlie_1(1,:);
for kkk=2:10
total_1=zhenlie_1(kkk,:)+total_1;
end
figure(7);
plot(1:length(total_1),total_1);
xlabel('时间‘s');ylabel('幅度');title('阵列接收到的亮点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);
doa=atan(5/5000); %doa为线阵法线方向与声源的夹角
w = exp(j*2*pi*k*m'*sin(doa*pi/180)); %每个阵元的时延tow=k*m‘*sin（）不同
temp=zeros(M,length(total_1));
for kkk=1:10 temp(kkk,1)=w(kkk,1); end ww=temp*total_1';
theta = linspace(-90,90,360); % 扫描方位角 a = exp(-j*2*pi*k*m'*sin(theta*pi/180)); %theta是扫描角度% 导向适量
Y = abs(ww'*a); Y = 20*log10(Y/max(Y)); Y_max=max(Y);
figure(8);
plot(theta,Y);
title('未加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');
%=============================================================
%加权
quanzhi=[zeros(1,90) 0.5*ones(1,45) 1*ones(1,90) 0.5*ones(1,45) zeros(1,90)];
tempp=zeros(360,360);
for kkk=1:360
tempp(1,kkk)=quanzhi(1,kkk);
end
YY = abs(Y*tempp);
YY = 20*log10(YY/max(YY));
figure(9);
plot(theta,YY);
title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');
%============================================================================
YYY=2.^(3*Y);
YYY = 20*log10(YYY/max(YYY));
YYY_max=max(YYY);
figure(10);
plot(theta,YYY);
title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');
%==========================================================
然后依次写亮点2和亮点3的程序，然后把得到的total_1,total_2,total_3画在一张图上，得到了
感觉生无可恋。。。。。。。。