源码
/qq_44394952/article/details/124490764?spm=1001..3001.5502
1.实验目的
（1）学习并熟悉Bellohop模型的原理
（2）使用Bellohop模型产生水声信道，如通信距离、水深、浪高等，进行仿真
（3）实现水声信道调制技术仿真，采用相干检测的方法进行PSK、QAM调制解调，分析水声通信的特点
2.Bellohop模型
Bellohop模型是通过高斯波束跟踪方法，计算水平非均匀环境中的声场。高斯波束跟踪方法对于高频水平变化问题特别有吸引力，这是简正波、波数积分和抛物线模型不可替代的。高斯束射线跟踪法的基本思想是将高斯强度分布与每条声线联系起来，该声线为高斯声束的中心声线,这些声线能较平滑的过渡到声影区，也能较平滑的穿过焦散线，所提供的结果与全波动模型的结果更为一致。在频率范围为600Hz-30kHz时与实验数据及理论模型符合得很好。
Bellohop模型产生信道的过程如下：
（1）BELLHOP在通过相关命令来读取可以描述发送端和接收端的环境和几何特性文件。mallab画图程序(plotssp，plotbty，plotbrc,etc.)用来显示每一个输入文件。输入文件有：
*.env:描述水体环境，必选；
*.ati:描述海面形状，可选；
*.bty:描述海底形状，可选；
*.trc:描述顶部反射系数，可选；
*.brc:描述底部反射系数,可选。
（2）BELLHOP在整体环境文件中通过相关选项来产生不同的输出文件：
*.ray:描述声线和本征声线;
*.shd:描述声线的传播损失;
*.arr:描述声线到达的时间-幅度序列等参数
系统的冲激响应仿真结果图：
3.系统组成及原理
（1）信源产生待传输的序列
信源发射端使用randsrc()产生一定长度的串行的二进制序列
（2）调制
1）16QAM调制
2/4电平转换由二进制变为四进制数字基带信号并进行格雷码映射，取值为-3，-1,1,3。格雷码的特点是邻近的两个信号，其标识仅有一位不同，这样当误符号率一定时，误比特率最小。
插值：在QAM符号之间插入一些0点，和后面的低通滤波联合起来的作用是模拟信号高频调制之后的频谱。分别在I路和Q路中，相邻的两个码字之间添加7个0。
低通：用于调整QAM符号的波形，使两路QAM符号分别经过平方根升余弦滤波器，形成等待调制到高频的信号。
两路QAM多进制信号分别和两路正交载波相乘。
2）4PSK调制
4PSK利用4种不同的相位表征数字信息，将二进制数字序列中每两个比特分为一组，00、01、10、11分别对应相位1/4π、3/4π、5/4π、7/4π。
（3）通过信道
将信号与Bellohop模型水声信道产生的单位冲激响应信号卷积，并通过awgn()函数加入高斯白噪声
（4）解调
1）16QAM解调
接入端接收到的QAM信号分别使用两路相干载波进行解调：
抽样：用于还原为调制时插值之前的QAM符号序列数目。经过插值，4QAM序列的点数是原本序列的4倍，16QAM序列的点数是原本序列的8倍，需要去除这些冗余点。
判决：还原QAM符号序列为二进制码元序列。通过根据最大后验概率准则，得到最小检测距离。
按照调制中的格雷码映射规则，经过“L-2电平变换”转换成还原为原来的I路和Q路码元信息。
对两路信号进行并/串变换，将两路码元数据合并后形成完整的比特流数据，
I路信号作为还原后信号中序号为奇数的码元，Q路则作为还原后信号中序号为偶数的码元。
2）4PSK解调
接收端根据判断出的相位，确定发送端发送的数字序列。若收到的相位靠近1/4π，则认为发送端发送的是00；若收到的相位靠近3/4π，则认为发送端发送的是01；若收到的相位靠近5/4π，则认为发送端发送的是10；若收到的相位靠近7/4π，则认为发送端发送的是11。
（5） 计算误码率
使用biterr()函数比较解调后输出的比特流与发送端产生的比特流，计算误码率，绘制仿真误符号率与理论误符号率曲线，进行对比分析。
4.模块源代码及各模块产生波形
（1）Bellohop模型水声信道产生单位冲激响应信号
代码：
filename = 'C:\CMST Software\AcTUP v2.2L\AcTUP\Output\Bounce&Bellhop\Test1109_00100.arr'; Minimum_range=100 %（接收水听器的水平方向上接收范围最小值，m）Maximum_range=1000 %（接收水听器的水平方向上接收范围最大值，m）[ amp1, delay, SrcAngle, RcvrAngle, NumTopBnc, NumBotBnc, narrmat, Pos ]... = read_arrivals_asc(filename);%%单位冲激响应[m,n]=size(amp1);amp = abs(amp1); %取模 x = delay(m,:); %获取第50个接收机的时延和幅值y = amp(m,:);%%归一化冲激响应Amp_Delay = [x;y];Amp_Delay(:,all(Amp_Delay==0,1))=[]; %去掉0值Amp_Delay=sortrows(Amp_Delay',1); %按照时延从小到大排序normDelay = Amp_Delay(:,1)-Amp_Delay(1,1);%归一化时延normAmp = Amp_Delay(:,2)/Amp_Delay(1,2);%归一化幅度%距离时延figure(3)mum=1:m;ReRange = Minimum_range+(Maximum_range-Minimum_range)/m*mum;for i=1:min(narrmat)plot(delay(:,i),ReRange,'o')hold onendhold offgrid oncolorbarxlabel('时延(sec)')ylabel('Range(m)')title(filename)
波形：
（2）16QAM调制解调
1）调制
<1>产生随机二进制序列
%定义待仿真序列的长度 Nglobal NN=10000;% 首先产生随机二进制序列source=randsrc(1,N,[1,0]);
<2>串并转换
I路：第1,3,5,7...元素 Q路：第2,4,6,8...元素N=length(source);a=1:2:N;source1=source(a);%I路source2=source(a+1);%Q路
<3>2/4电平变换
a=1:2:N/2;temp1=source1(a);temp2=source1(a+1);y11=temp1*2+temp2;%I路 用十进制表示二进制temp1=source2(a);temp2=source2(a+1);y22=temp1*2+temp2;%Q路 用十进制表示二进制a=1:N/4;source1=(y11*2-1-4)*1.*cos(2*pi*a);source2=(y22*2-1-4)*1.*cos(2*pi*a);
<4>格雷码映射
%格雷码映射source1(y11==0)=-3; source1(y11==1)=-1;source1(y11==3)=1;source1(y11==2)=3;%I路source2(y22==0)=-3;source2(y22==1)=-1;source2(y22==3)=1;source2(y22==2)=3;%Q路
<5>插值
%两路信号source1、source2进行插值，插值的比例为8sig_insert1=zeros(1,8*length(source1));a=1:8:length(sig_insert1);sig_insert1(a)=source1;sig_insert2=zeros(1,8*length(source2));a=1:8:length(sig_insert2);sig_insert2(a)=source2;
<6>低通
[sig_rcos1,sig_rcos2]=rise_cos(sig_insert1,sig_insert2,0.25,2);function [y1,y2]=rise_cos(x1,x2,fd,fs)%生成平方根升余弦滤波器[yf, tf]=rcosine(fd,fs, 'fir/sqrt');%使用平方根升余弦滤波器对两路信号进行滤波[yo1, to1]=rcosflt(x1, fd,fs,'filter/Fs',yf);[yo2, to2]=rcosflt(x2, fd,fs,'filter/Fs',yf);y1=yo1;y2=yo2;
<7>将基带信号调制到高频上
f=0.25; %输入信号的频率hf=2.5; %载波的频率%用于存储插值后的输入信号yo1=zeros(1,length(sig_rcos1)*hf/f*10);yo2=zeros(1,length(sig_rcos2)*hf/f*10);n=1:length(yo1);yo1(n)=sig_rcos1(floor((n-1)/(hf/f*10))+1);yo2(n)=sig_rcos2(floor((n-1)/(hf/f*10))+1);t=(1:length(yo1))/hf*f/10;sig_modulatesig_modulate=yo1.*cos(2*pi*hf*t)-yo2.*sin(2*pi*hf*t);
2）解调
<1>判决
%对x1路信号进行判决xx1(sig_pick1>=2)=3;xx1((sig_pick1<2)&(sig_pick1>=0))=1;xx1((sig_pick1>=-2)&(sig_pick1<0))=-1;xx1(sig_pick1<-2)=-3;%对x2路信号进行判决xx2(sig_pick2>=2)=3;xx2((sig_pick2<2)&(sig_pick2>=0))=1;xx2((sig_pick2>=-2)&(sig_pick2<0))=-1;xx2(sig_pick2<-2)=-3;
<2>逆格雷码映射
%将x1路信号按格雷码映射规则还原成0、1信号temp1=zeros(1,length(xx1)*2); %创建全0矩阵temp1(find(xx1==-1)*2)=1; %矩阵第二列置一temp1(find(xx1==1)*2-1)=1; %矩阵第一列置一temp1(find(xx1==1)*2)=1;temp1(find(xx1==3)*2-1)=1;%将x2路temp2=zeros(1,length(xx2)*2);temp2(find(xx2==-1)*2)=1;temp2(find(xx2==1)*2-1)=1;temp2(find(xx2==1)*2)=1;temp2(find(xx2==3)*2-1)=1;
<3>并串转换
signal=zeros(1,length(temp1)*2);signal(1:2:length(signal))=temp1; %奇数位为I路signal(2:2:length(signal))=temp2; %偶数位为Q路
（3）4PSK调制解调
1）调制
代码：
function [modulated_X] = PSK4_modulation(X)modulated_X=zeros(1,length(X)/2);jj=0;for ii=1:2:length(X)jj=jj+1;if X(ii)==0&&X(ii+1)==0modulated_X(jj)=exp(1i*pi/4);elseif X(ii)==0&&X(ii+1)==1modulated_X(jj)=exp(1i*3*pi/4);elseif X(ii)==1&&X(ii+1)==0modulated_X(jj)=exp(1i*5*pi/4);else modulated_X(jj)=exp(1i*7*pi/4);endend end
2）解调
代码：
function [demodulated_X] =PSK4_demodulation(noised_X)demodulated_X=zeros(1,length(noised_X)*2);for k=1:length(noised_X)if(real(noised_X(k))>=0)&&(imag(noised_X(k))>=0)demodulated_X(k*2-1)=0;demodulated_X(k*2)=0;elseif(real(noised_X(k))<0)&&(imag(noised_X(k))>=0)demodulated_X(k*2-1)=0;demodulated_X(k*2)=1; elseif(real(noised_X(k))<0)&&(imag(noised_X(k))<0)demodulated_X(k*2-1)=1;demodulated_X(k*2)=0; elsedemodulated_X(k*2-1)=1;demodulated_X(k*2)=1; end endend
波形：
3）误码率
代码：
N=100000;X=randsrc(1,N,[0,1]);%获得水声信道冲击响应filename = 'C:\CMST Software\AcTUP v2.2L\AcTUP\Output\Bounce&Bellhop\Test1109_00100.arr'; [ amp1, delay1, SrcAngle, RcvrAngle, NumTopBnc, NumBotBnc, narrmat, Pos ]... = read_arrivals_asc(filename);%%单位冲激响应[m,n]=size(amp1);amp=abs(amp1(m,:));%取模 delay=delay1(m,:);%%归一化冲激响应Amp_Delay=[delay;amp];Amp_Delay(:,all(Amp_Delay==0,1))=[];%去掉0值Amp_Delay=sortrows(Amp_Delay',1);%按照时延从小到大排序normDelay=Amp_Delay(:,1)-Amp_Delay(1,1);%归一化时延normAmp=Amp_Delay(:,2)/Amp_Delay(1,2);%归一化幅度%4PSK调制modulated_X1=PSK4_modulation(X);Error_rate1=zeros(1,4000);SNR1=zeros(1,4000);modulated_X2=conv(modulated_X1,normAmp);modulated_X=modulated_X2(1:length(modulated_X1));%4PSK解调并画误码率曲线for n=1:4000noised_X1=awgn(modulated_X,n/16,'measured');power_of_signal1=sum(modulated_X.^2)/length(modulated_X);power_of_noise1=sum(abs(noised_X1-modulated_X).^2)/length(modulated_X);SNR1(n)=10*log10(power_of_signal1/power_of_noise1);demodulated_X1=PSK4_demodulation(noised_X1);Error_number1=sum(abs(demodulated_X1-X));Error_rate1(n)=Error_number1/N;end
5.仿真结果分析及结论
仿真的条件：
声源和接收水听器的水平距离：1000m
表面粗糙度的均方根值：2.0m
深度范围，即距海平面的垂直距离：[0 20]
横波声速：[1540 1543]
该层海水密度：[1021 1021]
横波吸收系数：[69.2912 69.4261]
把海底建模为声学弹性半空间，声速不变，密度均匀
忽略底部水体表面的粗糙度
沉积物层深度：20m
沉积物层的声速：2000m/s
泥沙密度：1810kg/m3
横波吸收系数：0.5
声波频率：100Hz
声源深度：10m
接收水听器的深度：10m
接收水听器的水平方向上接收范围最小值：100m
接收水听器的水平方向上接收范围最大值：1000m
接收水听器的水平方向上个数：50
波束类型：高斯波束
声线个数：10
第一个声线的到达角度：-15
最后一个声线的到达角度：15
结论:
由水声信道的单位冲激响应和归一化冲激响应，可以看出水声信道的距离时延以及衰减情况。在距离改变时会出现多径效应，会使接收端接收到的信号产生畸变。水声通信在不使用其他技术直接传输时，误码率非常高，在仅添加调制、解调技术时误码率也比较高，想要降低误码率还需要再加入均衡、分集、扩频等技术。而且在相同的信噪比条件下，使用4PSK比使用16QAM误码率低，4PSK抗噪声性能更好。
6.心得体会
通过本次水声通信课程实验，我巩固了课堂上老师讲过的知识，尤其是水声通信中的调制解调部分，同时还提升了matlab编程能力。在实践应用过程中，了解到了BELLHOP模型，BELLHOP模型是利用高斯波束跟踪方法代替几何波束跟踪方法是该模型与传统射线模型的主要区别。它克服传统射线模型中声影区强度为0和焦散线截面为0处声强度为无穷大的缺陷。缺点是确定性模型，不能反映系统的时变特性，不能处理水平变化的声速问题。本次水声通信课程实验使我受益良多。
【水声通信】使用Bellohop模型产生水声信道 采用相干检测的方法进行PSK QAM调制解调【matlab代码】
如果觉得《【水声通信】使用Bellohop模型产生水声信道 采用相干检测的方法进行PSK QAM调制解调》对你有帮助，请点赞、收藏，并留下你的观点哦！
