背景及原理讲解链接: https://blog.csdn.net/qq_44394952/article/details/124489788?spm=1001.2014.3001.5502.
BELLHOPdata.m

clear
clc
filename = 'C:\CMST Software\AcTUP v2.2L\AcTUP\Output\Bounce&Bellhop\Test20211109_00100.arr'; %此处文件名改成自己产生的
Minimum_range=100  %(接收水听器的水平方向上接收范围最小值,m)----R
Maximum_range=1000 %(接收水听器的水平方向上接收范围最大值,m)---RB [ 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,:);figure(1)
stem(x,y)
grid on
xlabel('相对时延/s')
ylabel('幅度')
title('单位冲激响应')%%归一化冲激响应
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(2)
stem(normDelay,normAmp,'^')
grid on
xlabel('相对时延/s')
ylabel('归一化幅度')
title('归一化冲激响应')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 on
end
hold off
grid on
colorbar
xlabel('时延(sec)')
ylabel('Range(m)')
title(filename)

4PSK
PSK4.m

clear all
N=100000;
X=randsrc(1,N,[0,1]);%获得水声信道冲击响应
filename = 'C:\CMST Software\AcTUP v2.2L\AcTUP\Output\Bounce&Bellhop\Test20211109_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
figure;
semilogy(SNR1,Error_rate1);
hold on;
grid on;
legend('4PSK误码率');
xlabel('信噪比SNR(dB)');
ylabel('误码率Pe');
xlim([-25,-10]);
ylim([0.3,0.4]);

PSK4_modulation.m

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);endendend

PSK4_demodulation.m

function [demodulated_X] =PSK4_demodulation(noised_X)
%4PSK解调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;  endendend

16QAM
qam16.m

% 16QAM的调制与解调示例程序。
% 程序输出了上述过程中的图形,并且在控制台中输出了相关的信息。
clc;
clear all;
close all;%==============================
%定义待仿真序列的长度 N
global N
N=10000;
% 首先产生随机二进制序列
source=randsrc(1,N,[1,0]);
%画出序列波形
figure(1);
stairs(source);
xlim([0,40]);
ylim([-0.5,1.5]);
title('随机二进制序列');%==============================
%对产生的二进制序列进行QAM调制
%1.串并转换将原二进制序列转换成两路信号 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路%画出串并变换后序列波形
figure(2);
plot_2way_IQ(source1,source2);
title('串并变换后I路、Q路符号序列');%2.分别对两路信号进行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);%格雷码映射
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路% 画出序列波形
figure(3);
plot_2way_stair(source1,source2);
title('16QAM符号序列'); %画出星座图
figure(4);
plot_astrology(source1,source2);
title('16QAM星座图');%画出频域图
figure(5);plot_2way_dtft(source1,source2);
title('16QAM符号频域图');%===============================
%两路信号source1、source2进行插值,插值的比例为8
sig_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;%画出两路信号的波形图
figure(6);
plot_2way(sig_insert1,sig_insert2,length(sig_insert1),0.5);
title('插值后I、Q两路信号的波形图');%画出频域图
figure(7);
plot_2way_dtft(sig_insert1,sig_insert2);
title('插值16QAM符号频域图');%===============================
%通过低通滤波器
[sig_rcos1,sig_rcos2]=rise_cos(sig_insert1,sig_insert2,0.25,2);%画出两路信号的波形图
figure(8);
plot_2way(sig_rcos1,sig_rcos2,length(sig_rcos1)/4,0.5);
title('通过低通滤波器后两路信号波形图');%画出频域图
figure(9);
plot_2way_dtft(sig_rcos1',sig_rcos2');
title('通过低通滤波器后频域图');%===============================
%将基带信号调制到高频上
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);%对输入信号分别进行插值,相邻的两个点之间加入9个点,且这9个点的值同第0个点的值相同
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_modulate
sig_modulate=yo1.*cos(2*pi*hf*t)-yo2.*sin(2*pi*hf*t);%画出基带信号调制到高频之后的波形图
figure(10);
plot(t,sig_modulate);
xlim([0,800]);
title('基带信号调制到高频之后的波形图');%===============================
%加入高斯白噪声
%设信噪比
snr=10;%符号信噪比
snr1=snr+10*log10(4);%所有信号的平均功率
%var()方差函数
ss=var(sig_rcos1+j*sig_rcos2,1);%加入高斯白噪声
y=awgn([sig_rcos1+j*sig_rcos2],snr1+10*log10(ss/10),'measured');
%合起来的信号分为实部和虚部,实部为I路,虚部为Q路,再整体加高斯
x1=real(y);%对得到的信号,取实部为x1
x2=imag(y);%虚部为x2%   ' 转置矩阵 确保第一个矩阵中的列数与第二个矩阵中的行数匹配
sig_noise1=x1';
sig_noise2=x2';%画出I路和Q路加噪声后的波形图
figure(11);
plot_2way(sig_noise1,sig_noise2,length(sig_noise1)/4,0.5);
title('加入高斯白噪声之后的波形图');%画出频域图
figure(12);
plot_2way_dtft(sig_noise1,sig_noise2);
title('加入高斯白噪声之后频域图');%===============================
%经过匹配滤波器
[sig_match1,sig_match2]=rise_cos(sig_noise1,sig_noise2,0.25,2);%画出经过匹配滤波器之后的波形图
figure(13);
plot_2way(sig_match1,sig_match2,length(sig_match1)/4,0.5);
title('经过匹配滤波器之后的波形图');%画出频域图
figure(14);
plot_2way_dtft(sig_match1',sig_match2');
title('经过匹配滤波器之后DTFT频域图');%===============================
%抽样
sig_pick1=sig_match1(8*3*2+1:8:(length(sig_match1)-8*3*2));
sig_pick2=sig_match2(8*3*2+1:8:(length(sig_match1)-8*3*2));%画出星座图
figure(15);
plot_astrology(sig_pick1,sig_pick2);
title('接收的16QAM符号星座图');%画出序列波形
figure(16);
plot_2way_stair(sig_pick1,sig_pick2);
title('接收的16QAM符号序列');%画出频域图
figure(17);
plot_2way_dtft(sig_pick1',sig_pick2');
title('接收的QAM符号序列频域图');%===============================
%QAM解调
%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;%画出判决后16QAM符号序列
figure(18);
plot_2way_stair(xx1,xx2);
title('判决后16QAM符号序列'); %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路信号按格雷码映射规则还原成0、1信号
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;%画出4/2电平变换后I路、Q路符号序列
figure(19);
plot_2way_IQ(temp1,temp2);
title('4/2电平变换后I路、Q路符号序列');%3.并串转换:将两路0、1信号合成一路
signal=zeros(1,length(temp1)*2);
signal(1:2:length(signal))=temp1; %奇数位为I路
signal(2:2:length(signal))=temp2; %偶数位为Q路%画出接收端恢复的二进制序列
figure(20);
stairs(signal);
xlim([0,40]);
ylim([-0.5,1.5]);
title('接收端恢复的二进制序列');%===============================
%误码率计算
[err_number,Pe] = biterr(source,signal);
fprintf('误码数以及误码率:\n');
fprintf('%d,%f;\n',err_number,Pe);
hold on

rise_cos.m

%x1、x2是两路输入信号,fd是信号信息位的频率,fs是信号的采样频率
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;

plot_2way.m

function y=plot_2way(x1,x2,len,t)
subplot(2,1,1);
plot((1:len)*t,x1(1:len));
axis([0 160 -4 4]);
hold on
plot((1:len)*t,x1(1:len),'.','color','red');
xlabel('实部信号');
hold off
subplot(2,1,2);
plot((1:len)*t,x2(1:len));
axis([0 160 -4 4]);
hold on
plot((1:len)*t,x2(1:len),'.','color','red');
hold off
xlabel('虚部信号');

plot_2way_dtft.m

function y=plot_2way_f_dtft(x1,x2)n=0:1:length(x1)-1;                          %原信号是1行8列的矩阵
w=[-800:1:800]*4*pi/800;            %频域共-800----+800 的长度    X1=x1*exp(-j*(n'*w));               %求dtft变换,采用原始定义的方法,对复指数分量求和而得
subplot(2,1,1);
plot(w/pi,abs(X1));
xlabel('实部信号(单位是π)');X2=x2*exp(-j*(n'*w));
subplot(2,1,2);
plot(w/pi,abs(X2));
xlabel('虚部信号(单位是π)');

plot_2way_IQ.m

function y=plot_2way_IQ(x1,x2) subplot(2,1,1);
hold on;
stairs(x1);
hold off;
xlim([0,20]);
ylim([-0.5,1.5]);
xlabel('I路序列');subplot(2,1,2);
hold on;
stairs(x2);
hold off;
xlim([0,20]);
ylim([-0.5,1.5]);
xlabel('Q路序列');

plot_2way_stair.m

function y=plot_2way_stair(x1,x2)subplot(2,1,1);
hold on;
stairs(x1);
hold off;
xlim([0,20]);
ylim([-4,4]);
xlabel('实部信号');subplot(2,1,2);
hold on;
stairs(x2);
hold off;
xlim([0,20]);
ylim([-4,4]);
xlabel('虚部信号');

plot_astrology.m

function plot_astrology(a,b)
%画出星座图
subplot(1,1,1)
plot(a,b,'*');
axis([-5 5 -5 5]);
line([-5,5],[0,0],'LineWidth',2,'Color','black');
line([0,0],[-5,5],'LineWidth',2,'Color','black');
title('16QAM星座图');

【水声通信】使用Bellohop模型产生水声信道,采用相干检测的方法进行PSK、QAM调制解调【matlab源码】相关推荐

  1. 【水声通信】使用Bellohop模型产生水声信道,采用相干检测的方法进行PSK、QAM调制解调【matlab代码】

    源码 https://blog.csdn.net/qq_44394952/article/details/124490764?spm=1001.2014.3001.5502 1.实验目的 (1)学习并 ...

  2. 【肺实质分割】基于主动轮廓模型和贝叶斯方法识别实现胸膜旁肺实质分割附matlab代码

    ✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信.

  3. 沿全尺寸模型水平轴 MHK 涡轮机(DOE RM1)叶片性能表征计算(Matlab代码实现)

               目录

  4. 灰色预测模型【GM(1,1)模型】 【matlab代码】

    灰色系统介绍 灰色系统是由华中科技大学的邓聚龙教授于80年代初创立,该系统作为新兴的横断学科,在短短的二十年里已得到了长足的发展. 其已经成为社会,经济,科教,科技等很多领域进行预测,决策,评估,规划 ...

  5. matlab中频率选择性衰落信道,浅海水声信道模型分析及频率选择性衰落仿真

    第10卷第12期 Vo l .10,No .12 宜宾学院学报 J ou rnal of Yibin Un i versity 2010年12月 Dec .,2010 收稿日期1修回作者简介贺繇(), ...

  6. matlab对声场仿真,基于声线模型的水声传播MATLAB仿真.pdf

    基于声线模型的水声传播MATLAB仿真.pdf 科技广场 2007.9 基于声线模型的水声传播MATLAB仿真 熊光耀 杨 琴 Ciong GuangyaoYang Qin (江西中医学院计算机系,江 ...

  7. DL之模型调参:深度学习算法模型优化参数之对深度学习模型的超参数采用网格搜索进行模型调优(建议收藏)

    DL之模型调参:深度学习算法模型优化参数之对深度学习模型的超参数采用网格搜索进行模型调优(建议收藏) 目录 神经网络的参数调优 1.神经网络的通病-各种参数随机性 2.评估模型学习能力

  8. 魔术轮胎,dugoff轮胎建模 采用模块化建模方法,搭建非线性魔术轮胎PAC2002,dugoff模型

    魔术轮胎,dugoff轮胎建模 软件使用:Matlab/Simulink 适用场景:采用模块化建模方法,搭建非线性魔术轮胎PAC2002,dugoff模型. 非线性轮胎模型输入: 轮胎侧偏角,轮胎滑移 ...

  9. 该模型在额定以下采用MTPA控制,速度环输出给定电流,然后代入MTPA得到dq电流,电压反馈环输出超前角进行弱磁

    该模型在额定以下采用MTPA控制,速度环输出给定电流,然后代入MTPA得到dq电流,电压反馈环输出超前角进行弱磁. PI控制采用抗积分饱和,SVPWM考虑过调制情况,附带参考文献 ID:5858675 ...

最新文章

  1. pku 2418 Hardwood Species 字典树
  2. java哪个城市的需求量大_4大互联网热门城市Java薪资情况,看完你想去哪个城市发展呢?...
  3. perl 安装html,centos perl 安装HTML-Parser时报错
  4. d3 v5 api select
  5. Spring Security OAuth2 授权失败(401)
  6. 无线传感器网络 | 期末复习知识点1
  7. tab栏切换制作(原生js版本)
  8. Java 使用 POI 对 Excel文件 进行读写操作
  9. Button 自动换行
  10. 秒懂C#通过Emit动态生成代码
  11. QC3.0充电器快充诱骗方法,做个笔记
  12. 基于Labview的串口通信助手,附带免费的exe程序和vi文件(有意见或学习讨论欢迎交流)
  13. 海马玩模拟器无法链接问题处理
  14. 装饰者模式 增加功能;动态代理减少功能 只要完成自己部分功能 (繁杂部分交给他人处理)...
  15. java秃顶_【本人秃顶程序员】在Java中使用函数范式提高代码质量
  16. Medium之1421.净现值查询
  17. Java专题 基础篇--判断(三元表达式,switch等) +个税计算案例
  18. C#DataGridView导入Excel数据,并上传数据
  19. Bloom filters in Python
  20. python人工智能是做什么的,人工智能需要学python吗

热门文章

  1. Matlab------------怎么取一个复数的实部和虚部
  2. 解决百度网盘限速的问题
  3. Redis 事务和事务锁
  4. java实现图片中内容读取实现方案之阿里云OCR
  5. 《角斗士》一个帝国的史诗绝唱
  6. 深度学习CUDA安装失败及解决方案
  7. LaTex 编译中文
  8. 第二届先导杯-在曙光超算平台编译cp2k
  9. “error : unknown filesystem”的解决的方法
  10. oracle如何上个月和下个月的月份