水声信道仿真_matlab

function [signal]=sdoppler(delay,fading,t,g1a,g2a) sv=5; %*** the speed of transmitor(m) rv=6; %*** the speed of receiver(m) a_c=1520; %*** the speed of acoustic(m) fs=2*length(t); f1=9000; f2=11000; Ts=0.001; %*** the width of the code element N=length(delay); DN=3; L=length(t); f_p=40;f_s=80;R_p=3;R_s=30; %*** parameters of buttord Ws=2*f_s/fs;Wp=2*f_p/fs; %*** frequency normalization [n,Wn]=buttord(Wp,Ws,R_p,R_s); %*** using buttord [b,a]=butter(n,Wn); %*** coefficients of buttord Fs=50000; tt=0:1/Fs:0.13-1/Fs; signal=zeros(1,length(tt)); delay=round(10000*delay)/10000; %delay=delay+3*Ts %delay(1)=delay(1)-0.00001; %N acoustic path for i=1:N R1=raylrnd(40); %*** produce the frequency of doppler R2=raylrnd(40); f1_shift=f1+R1; f2_shift=f2+R2; t1=t+delay(i); %the main signal d_f1=fading(i)*cos(2*pi*f1_shift*t1-2*pi*f1*delay(i)); %*** doppler and propagation delay d_f2=fading(i)*cos(2*pi*f2_shift*t1-2*pi*f2*delay(i)); %*** doppler and propagation delay d_f11=fading(i)*sin(2*pi*f1_shift.*t1-2*pi*f1*delay(i)); d_f22=fading(i)*sin(2*pi*f2_shift.*t1-2*pi*f2*delay(i)); d_signal=g1a.*d_f1+g2a.*d_f2; k=1:3; M=exp(-3*k/DN); nt=zeros(1,L); nnt=zeros(1,L); %*** the multipath signal ***% for k=1:3 Md_f1=fading(i)*cos(2*pi*f1_shift*t1-2*pi*f1*delay(i)-k*Ts); %*** doppler and propagation delay Md_f2=fading(i)*cos(2*pi*f2_shift*t1-2*pi*f2*delay(i)-k*Ts); %*** doppler and propagation delay Md_f11=fading(i)*sin(2*pi*f1_shift.*t1-2*pi*f1*delay(i)-k*Ts); Md_f22=fading(i)*sin(2*pi*f2_shift.*t1-2*pi*f2*delay(i)-k*Ts); Md_signal1=g1a.*Md_f1+g2a.*Md_f2; %*** cophase component of multipath signal Md_signal2=g1a.*Md_f11+g2a.*Md_f22; %*** orthogonal component of multipath signal n1=M(k)*randn(1,L); n1=filter(b,a,n1); ec=sqrt(n1*n1'/L); n1=M(k)/ec*n1; %*** keep the variance unchanged M(k)*M(k) n2=M(k)*randn(1,L); n2=filter(b,a,n2); ec=sqrt(n2*n2'/L); n2=M(k)/ec*n2; if i==1 noise=n2; noiset=Md_signal1.*n1+Md_signal2.*n2; end nt=nt+Md_signal1.*n1+Md_signal2.*n2; end d_signal=d_signal+nt; %*** 由于个本征信号有延迟，所以在信号叠加时，必须在信号前与信号后 ***% %*** 进行补零，使得个本征信号长度一样 ***% zt=0:1/Fs:delay(i); %*** zeros pad before signal zz=0.001*randn(1,length(zt)); d_signal=[zz,d_signal]; zt=0:1/Fs:(delay(i)+0.05); figure(i+2); %subplot(211); plot(zt,d_signal); xlim([0,delay(i)+0.06]); xlabel('time(s)'); ylabel('amplitute'); %zt=0.05+delay(i):1/Fs:0.13-1/Fs; %**** zeros pad afer signal zz=zeros(1,length(signal)-length(d_signal)); d_signal=[d_signal,zz]; %subplot(212); %N=2^(ceil(log2(length(t)))); %n=0:N-1; %f=fs*n/N; %y=fft(d_signal,N); %yy=abs(y); %yy=yy*2/N; %plot(f,yy); signal=signal+d_signal; end %*** 画出个本征信号叠加后总的波形 ***% figure(N+3); plot(tt,signal); xlabel('time(s)'); ylabel('amplitute') %*** 仿真一条再带高斯噪声 ***% %*** 以及信号经过该高斯噪声调制后的波形 ***% m=0:L-1; figure(N+4); plot(m/L*0.05,noise); xlim([0,0.052]); xlabel('time(s)'); ylabel('amplitute') figure(N+5); plot(m/L*0.05,noiset); xlim([0,0.052]); xlabel('time(s)'); ylabel('amplitute') %chongji=xcorr(signal,d_signal); %figure(N+6); %plot(chongji)