这次试着处理一下多亮点。
首先把直达，海面反射一次，海底反射一次的程序弄成一个function的函数，进行调用，不然代码太多了，看着很头疼。。
%=================================================================================================================
function receive_signal = Com_Sonar_data_generation(c,fs,f0,Tao,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
% Tao 信号脉宽 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 考虑最大的反射次数
%***************************************输出参数说明******************************************
% receive_signal 以信号发射为起始采样点的水听器接收信号
%==================================================
Ts = 1/fs; %采样时间间隔
sample_num = fix(Sample_time*fs); %采样总点数 FIX：让x向0靠近取整 %
nTs = (0:sample_num-1)/fs; %离散的采样时刻
sample_num1 = fix(Tao*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; %信号的实际到达时刻
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)) = S0*sin(2*pi*f0*nTs1+phase); %模拟采样 1*50050
receive_signal = receive_signal + Noise_var*randn(size(receive_signal)); %加入噪声
if Reflex_num>0 %考虑反射时
i = 1 ; %考虑1次反射的情况
%---------------------------------------------------第一次从水面反射时---------------------------------------------------
% 只考虑一次反射,海底和海面各一个反射点
d1=sqrt(D/(H1+H2)^2+1)*H1;
d2=sqrt(D/(H1+H2)^2+1)*H2; %从海面反射一次后到达的总路程：d1+d2
di=d1+d2;
Si = 1/di*Re_coef_surf; %反射波的声压幅值
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); %模拟采样
%---------------------------------------------------第一次从水底反射时---------------------------------------------------
d3=sqrt(D/(2*H-H1-H2)^2+1)*(H-H1);
d4=sqrt(D/(2*H-H1-H2)^2+1)*(H-H2);
dib = d3+d4; %反射波总路程
Sib = 1/dib*Re_coef_bottom; %反射波的声压幅值
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
%======================================================================================================================
%=========================================================================================================================
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; %考虑最大的反射次数
%=====================================
for temp=1:200 %频率循环
f0=1400+temp; %信号频率 Unit：Hz
fs=16000; %采样率 Unit：Hz
Tao=5/f0; %信号脉宽 Unit：s
%阵要每一次循环完初始化一下，不然运行不成功
signal_12=[]; signal_12_fft=[]; signal_22=[]; signal_22_fft=[]; signal_32=[]; signal_32_fft=[];
%============================亮点1===========================================
%第一过程:发射换能器到亮点1 % r=0;z=10;%发射信号的位置。
H1 = 10; %发射点水深 Unit：m
H2 = 15; %接收点水深 Unit：m
D = 5000; %接收与发射水平距离 Unit：m
signal_11=Com_Sonar_data_generation(c,fs,f0,Tao,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);
signal_11_fft=fft(signal_11);
%======================================================================
%第二过程：亮点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_12(kkk,:)=Com_Sonar_data_generation(c,fs,f0,Tao,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num); % temp：16，50053 signal_12_fft(kkk,:)=fft(signal_12(kkk,:));
end
zhenlie1_fft=zeros(1,length(signal_12_fft(1,:)));
%把该频率下，所有阵元接收到的信号加在一起
for kkk=1:10
K=signal_12_fft(kkk,:);
zhenlie1_fft=K+zhenlie1_fft;
end
% zhenlie1_fft=zhenlie1_fft.;
%=============================================================================
%从发射器到亮点1的频谱*亮点到阵列的信号频谱
signal_1_fft=signal_11_fft.*zhenlie1_fft;
%============================================================================
%============================亮点2===========================================
%第一过程:发射换能器到亮点2 % r=0;z=10;%发射信号的位置。 % r2=5010;z2=15;%2号亮点
H1 = 10; %发射点水深 Unit：m
H2 = 15; %接收点水深 Unit：m
D = 5010; %接收与发射水平距离 Unit：m
signal_21=Com_Sonar_data_generation(c,fs,f0,Tao,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);
signal_21_fft=fft(signal_21);
%======================================================================
%第二过程：亮点2到接收阵列 % r1=5010;z1=15;%2号亮点 % rr=0;zz=20:30;%接收信号
H1=15; %发射点水深 Unit：m
D=5010; %接收与发射水平距离 Unit：m
for kkk=1:10
H2=20+1*kkk; %接收点水深 Unit：m
signal_22(kkk,:)=Com_Sonar_data_generation(c,fs,f0,Tao,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);
signal_22_fft(kkk,:)=fft(signal_22(kkk,:));
end
zhenlie2_fft=zeros(1,length(signal_22_fft(1,:)));
for kkk=1:10
K=signal_22_fft(kkk,:);
zhenlie2_fft=K+zhenlie2_fft;
end
%==================================================
% 该频率下，发射器到亮点2的频谱*亮点2到阵元的频谱
signal_2_fft=signal_21_fft.*zhenlie2_fft;
%===========================================================
%============================亮点3===========================================
%第一过程:发射换能器到亮点3 % r=0;z=10;%发射信号的位置。 % r3=5050;z3=15;%3号亮点
H1 = 10; %发射点水深 Unit：m
H2 = 15; %接收点水深 Unit：m
D = 5050; %接收与发射水平距离 Unit：m
signal_31=Com_Sonar_data_generation(c,fs,f0,Tao,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);
signal_31_fft=fft(signal_31);
%======================================================================
%第二过程：亮点3到接收阵列 % r1=5020;z1=15;%3号亮点 % rr=0;zz=20:30;%接收信号
H1=15; %发射点水深 Unit：m
D=5050; %接收与发射水平距离 Unit：m
for kkk=1:10
H2=20+1*kkk; %接收点水深 Unit：m
signal_32(kkk,:)=Com_Sonar_data_generation(c,fs,f0,Tao,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);
signal_32_fft(kkk,:)=fft(signal_32(kkk,:));
end
zhenlie3_fft=zeros(1,length(signal_32_fft(1,:)));
for kkk=1:10
K=signal_32_fft(kkk,:);
zhenlie3_fft=K+zhenlie3_fft;
end
%==================================================
% 该频率下，发射器到亮点3的频谱*亮点3到阵元的频谱
signal_3_fft=signal_31_fft.*zhenlie3_fft;
%===========================================================
%把该频率下返回的声压全部加在一起
signal_fft(temp,:)=[signal_1_fft zeros(1,55000-length(signal_1_fft))]+[signal_2_fft zeros(1,55000-length(signal_2_fft))]+[signal_3_fft zeros(1,55000-length(signal_3_fft))];
end
%====================================================================================================
total_fft=signal_fft(1,:); for temp=1:200
total_fft =signal_fft(temp,:)+total_fft;
end
figure;
plot(1:length(total_fft),total_fft); total=ifft(total_fft);
figure;
plot(1:length(total),total);
