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

时间：2019-05-11 13:41:05 阅读：119 评论：0 收藏：0 [点我收藏+]

标签：第一个亮点 star 反射 mamicode size info 角度 tar

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

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

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

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--频域‘);

技术图片