基于matlab实现的多普勒频移海底混响点散射模型程序

完整程序：

 %有多普勒频移的海底混响点散射模型
clear all; close all; clc    

H=100;              %海水深度
D=50;               %合置声纳深度
c=1500;             %声速
azm=pi/6;          %水平方位角
u=-27;              %垂直散射系数
v=20;              %声纳运动速度
fs=25000;           %采样频率
f0=4000;            %中心频率
T=0.001;          %仿真步长
rou=0.5; % 0.85;    %散射体密度

startt=0.08;        %混响起始时间
endt=0.5;          %混响起始时间
t=0:1/fs:endt-1/fs;
Nt=length(t);
Rb1=zeros(size(t));

R1=zeros(size(t));

t1=0.05;

ts=0:1/fs:t1-1/fs;
s1=exp(1i*2.*pi.*f0.*ts);                 %CW信号

Ns=length(ts);
% 信号时域
figure(1);
plot(ts,real(s1));
axis([0 0.06 -1 1]);    %axis([xmin xmax ymin ymax])
xlabel('时间/s');
ylabel('幅度');
title('CW信号');

%信号频谱
X1=abs(fft(s1));

figure(2);

f=fs*(0:Ns/2)/Ns;
plot(f,X1(1:Ns/2+1));
xlabel('频率/Hz');
title('CW信号频谱')
grid on;%添加网格

m=20; %m次运算取平均
for i=1:m
    for tr=startt:T:endt
        r1=tr*c/2;
        r2=(tr+T)*c/2;
        As=(r2^2-r1^2)*azm/2;
        Na=poissrnd(As*rou,1,1);  %产生均值为面积的按泊松分布的一个数
        r=r1+rand(1,Na)*(r2-r1);
        phi=rand(1,Na)*2*pi;
        alf=-pi/12+pi/6.*rand(1,Na);
        for n=1:Na
                d=r(n);
                %计算多普勒频移
                cosbt=(sqrt((r(n))^2-D^2))/r(n);
                fd1=2*v/c*f0*cos(alf(n))*cosbt;
                s1p=exp(1i*2.*pi.*(f0+fd1).*ts);%多普勒频移相当于对入射波进行了频移

                Nrt=fix(d/c*2*fs);%取整
                Sr=sqrt(10^(u/10))*(H-D)/d;     %散射损失，采用兰伯特（Lambert）定律
                Fr=1/d;                         %传播损失
                h=Sr*Fr^2*exp(1i*phi(n));                
                if Nrt+Ns>Nt
                    Rb1(Nrt+1:Nt)=Rb1(Nrt+1:Nt)+s1p(1:Nt-Nrt).*h;
                else
                    Rb1(Nrt+1:Nrt+Ns)=Rb1(Nrt+1:Nrt+Ns)+s1p.*h;
                end
        end
    end
    R1=R1+Rb1;
end
%混响时域信号
R1=R1/m;
Rr1=real(R1);

figure(3);
plot(t,Rr1);
xlabel('时间/s');
ylabel('归一化瞬时值');
title('CW信号的混响时域信号');

%混响频谱
F1=abs(fft(Rr1));                
figure(4);
f1=fs*(0:Nt/2)/Nt;
plot(f1,F1(1:Nt/2+1));
xlabel('频率/Hz');
ylabel('归一化瞬时值');
title('CW信号的混响频域信号');

%%%%%%求多普勒%%%%%%%%%%%%%%%%%%%%%%%
tao=0.002;
fai=pi/6;
%%%找底位置%%%%%%%
tc=0:1/fs:endt-1/fs;
cpy=exp(1i*2*pi*f0*tc).*(tc>0 & tc<t1);
CPY=conj(fft(cpy));

R2=fft(R1);
crr=ifft(R2.*CPY);

figure(5)
plot(t,real(crr))

[~,bt]=max(real(crr));

tn=tao*fs;
x1=R1(bt:end-tn);
x2=R1(bt+tn:end);
xout=x1*x2';
fb=abs(angle(xout)/2/pi/tao)
vout=fb*c/2/cos(fai)/f0