水声多径信道下LMS仿真(Python代码)

2023-10-30 19:10
文章标签 python 代码 仿真 水声 信道 多径 lms

本文主要是介绍水声多径信道下LMS仿真(Python代码)，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

1. 模块的导入

做仿真，numpy应该都知道
arlpy是水声通信工具箱，不仅可以产生信道的冲激响应，还有一些通信相关的函数
bokeh是用来画图的，配合jupyter notebook使用很爽！！（Matplotlib是常用的绘图工具包，个人感觉bokeh更舒服…）

# 数组运算工具包
import numpy as np
# 水声通信工具包
import arlpy.uwa as au
import arlpy.uwapm as aup
import arlpy.comms as ac
# 绘图工具包
from bokeh.plotting import show, figure
from bokeh.layouts import gridplot
from bokeh.io import output_notebook
output_notebook()
import warnings
warnings.filterwarnings('ignore')

2. 准备工作

2.1 水声信道的冲激响应

如何设置参数，得到冲激响应、衰减、时延等等在arlpy的官方文档中，都有介绍，可自行查阅
https://arlpy.readthedocs.io/en/latest/

# 设置环境参数
env = aup.create_env2d()
env['bottom_soundspeed'] = 1800
env['bottom_absorption'] = 0
env['depth'] = 100
env['max_angle'] = 30
env['min_angle'] = -30
env['frequency'] = 100
env['nbeams'] = 10
env['tx_depth'] = 20
env['rx_depth'] = 20
env['rx_range'] = 1000
# 假设声速在浅海环境下是递增的
sound_speed = np.zeros(shape=[6, 2])
for i in range(6):sound_speed[i, 0] = 20 * isound_speed[i, 1] = au.soundspeed(depth=20 * i)
env['soundspeed'] = sound_speedarrivals_sd = aup.compute_arrivals(env) # 声线到达的时间
delay_sd = min(arrivals_sd.time_of_arrival) # 端到端的传输时延
multidelay_sd = max(arrivals_sd.time_of_arrival) - \min(arrivals_sd.time_of_arrival) # 多径时延
fs = 100  # 信道的采样频率
cir_sd = aup.arrivals_to_impulse_response(arrivals_sd, fs=fs)
cir_sd = np.real(cir_sd)
cir_sd[(cir_sd < 0)] *= -1 # 信道的冲激响应

2.2 归一化及均衡代码

考虑到信道衰减比较大，在对信号进行处理的时候，需要考虑归一化的问题

def nor(x):'''param x: 带归一化信号return: 已归一化信号'''alpha = np.sum(np.square(x))/x.sizereturn x/np.sqrt(alpha)def equlization(sig_input, sig_expect, mu, M):'''param sig_input: 输入待均衡信号param sig_expect: 期望信号param mu: 均衡器步长param M: 均衡器长度return: (均衡后信号, 均方误差)'''n = sig_expect.sizeW = np.zeros(M)error = np.zeros(n)sig_output = np.zeros(n)if sig_input.size < n + M:sig_input = np.pad(sig_input, (0, n+M-sig_input.size))for i in range(n):input1 = sig_input[i:i+M]output = np.dot(W.T, input1)sig_output[i] = outputerror[i] = sig_expect[i] - outputW = W + mu * error[i] * input1J = np.square(error)return sig_output, J

3. 仿真

仿真的具体参数在调试时，有稍微改动，可自行调试

3.1 发射端

发射端使用了根升余弦滤波器进行滤波

N = 10000  # 数据点数
span = 4  # 升余弦滤波器范围
sps = 6  # 每符号采样率
rrc = ac.rcosfir(beta=0.5, sps=sps, span=span)
delay_rrc = int((rrc.size-1)/2)  # 经过滤波器之后，需要考虑滤波器的延迟sig_ori = ac.random_data(size=N, m=2) * 2 - 1
sig_exp = ac.upconvert(x=sig_ori, sps=sps, fs=0, fc=0)
sig_exp = nor(np.real(sig_exp))  # 使用ac.upconvert产生虚部为0的部分，为方便发图，取其实部
sig_tra = np.convolve(sig_exp, rrc)  # 发射滤波器的输出信号TOOLS = 'pan,box_zoom,xwheel_zoom,ywheel_zoom,crosshair,save,reset'  # 使用bokeh绘图时使用的交互工具
p1 = figure(tools=TOOLS, plot_width=800, plot_height=500,toolbar_location='above', toolbar_sticky=False)
p1.line(np.linspace(0, sig_exp.size-1, sig_exp.size),np.real(sig_exp), legend='sig_exp', line_color='red')
p1.line(np.linspace(0, sig_tra.size-1, sig_tra.size),np.real(sig_tra), legend='sig_tra-未考虑滤波器延迟',line_color='blue', line_dash='dashed')
p1.line(np.linspace(0, sig_tra.size-1-delay_rrc, sig_tra.size-delay_rrc),np.real(sig_tra[delay_rrc:]), legend='sig_tra-考虑滤波器延迟',line_color='green')
show(p1)

这个图是可以动的！

这个图是可以动的，横着拉，竖着拉，局部放大…（下面的图也是…）

3.2 多径信道

这里设置了三个不同的采样频率分别为100Hz，500Hz，1000Hz，得到的信道抽头个数分别为12,56,110，由仿真图可观察到采样频率越大，失真越严重

# fs = 100，信道长度为12
sig_isi1 = np.convolve(cir_sd, sig_tra)
sig_isi1 = nor(sig_isi1)p2 = figure(tools=TOOLS, plot_width=800, plot_height=300, title='信道长度为12',toolbar_location='left', toolbar_sticky=False)
p2.line(np.linspace(0, sig_exp.size-1, sig_exp.size),np.real(sig_exp), legend='sig_exp', line_color='red')
p2.line(np.linspace(0, sig_isi1.size-1-delay_rrc, sig_isi1.size-delay_rrc),np.real(sig_isi1[delay_rrc:]), legend='sig_isi-考虑滤波器延迟',line_color='green')fs = 500  # 信道的采样频率
cir_sd2 = aup.arrivals_to_impulse_response(arrivals_sd, fs=fs)
cir_sd2 = np.real(cir_sd2)
cir_sd2[(cir_sd2 < 0)] *= -1  # 信道的冲激响应
sig_isi2 = np.convolve(cir_sd2, sig_tra)
sig_isi2 = nor(sig_isi2)
p3 = figure(tools=TOOLS, plot_width=800, plot_height=300, title='信道长度为56',toolbar_location='left', toolbar_sticky=False)
p3.line(np.linspace(0, sig_exp.size-1, sig_exp.size),np.real(sig_exp), legend='sig_exp', line_color='red')
p3.line(np.linspace(0, sig_isi2.size-1-delay_rrc, sig_isi2.size-delay_rrc),np.real(sig_isi2[delay_rrc:]), legend='sig_isi-考虑滤波器延迟',line_color='green')fs = 1000  # 信道的采样频率
cir_sd3 = aup.arrivals_to_impulse_response(arrivals_sd, fs=fs)
cir_sd3 = np.real(cir_sd3)
cir_sd3[(cir_sd3 < 0)] *= -1  # 信道的冲激响应
sig_isi3 = np.convolve(cir_sd3, sig_tra)
sig_isi3 = nor(sig_isi3)
p4 = figure(tools=TOOLS, plot_width=800, plot_height=300, title='信道长度为110',toolbar_location='left', toolbar_sticky=False)
p4.line(np.linspace(0, sig_exp.size-1, sig_exp.size),np.real(sig_exp), legend='sig_exp', line_color='red')
p4.line(np.linspace(0, sig_isi3.size-1-delay_rrc, sig_isi3.size-delay_rrc),np.real(sig_isi3[delay_rrc:]), legend='sig_isi-考虑滤波器延迟',line_color='green')
show(p2)
show(p3)
show(p4)

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

为了方便观察，使用抽头长度为12的信道

3.3 接收端

接收方式：（做以下两种方式的测试）
（1）先滤波再均衡
（2）先均衡再滤波

M = int(cir_sd3.size)
mu = 0.0001
sig_rec1 = np.convolve(rrc, sig_isi3)
[sig_equ1, temp] = equlization(sig_rec1[delay_rrc*2:], sig_exp, mu, M)p5 = figure(tools=TOOLS, plot_width=800, plot_height=300,toolbar_location='above', toolbar_sticky=False)
p5.line(np.linspace(0, sig_exp.size-1, sig_exp.size),np.real(sig_exp), legend='sig_exp', line_color='red')
p5.line(np.linspace(0, sig_equ1.size-1, sig_equ1.size),np.real(sig_equ1), legend='sig-先滤波再均衡',line_color='green')[sig_equ2, temp] = equlization(sig_isi3[delay_rrc:], sig_exp, mu, M)
sig_rec2 = np.convolve(sig_equ2, rrc)[delay_rrc:]
p5.line(np.linspace(0, sig_rec2.size-1, sig_rec2.size),np.real(sig_rec2), legend='sig-先均衡再滤波',line_color='blue')
show(p5)

在这里插入图片描述

3.4 判决并得到误码率

3.4 判决并得到误码率

1、只与原序列进行比较判决
2、在下采样之前就进行判决

def decision1(sig_ori, sig_res, sps):err = 0for i in range(sig_ori.size):if sig_res[i*sps+int(sps/2)] < 0:temp = -1else:temp = 1if temp != sig_ori[i]:err += 1return errdef decision2(sig_exp,sig_res):err = 0for i in range(sig_exp.size):if sig_res[i] < 0:temp = -1else:temp = 1if temp != sig_exp[i]:err  += 1return errprint(decision1(sig_ori,sig_equ1,sps)/N,decision1(sig_ori,sig_rec2,sps)/N)
print(decision2(sig_exp,sig_equ1)/N/sps,decision2(sig_ori,sig_rec2)/N/sps)

0.0298 0.0699
1.0 0.08368333333333333

在不同的采样频率下，得到的信道抽头个数就不同，具体均衡器的步长需要进一步设定。
经测试，调整步长的之后，误码率的结果还看的过去（比上面好看！）…

一花

第一次发博，关于Python、仿真或者其他地方有不好的地方，请批评我！！！

这篇关于水声多径信道下LMS仿真(Python代码)的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/309984。

相关文章

python: 多模块(.py)中全局变量的导入

python: 多模块(.py)中全局变量的导入

文章目录 global关键字可变类型和不可变类型数据的内存地址单模块(单个py文件)的全局变量示例总结 多模块(多个py文件)的全局变量from x import x导入全局变量示例 import x导入全局变量示例 总结 global关键字 global 的作用范围是模块(.py)级别： 当你在一个模块（文件）中使用 global 声明变量时，这个变量只在该模块的全局命名空
阅读更多...
活用c4d官方开发文档查询代码

活用c4d官方开发文档查询代码

当你问AI助手比如豆包，如何用python禁止掉xpresso标签时候，它会提示到 这时候要用到两个东西。https://developers.maxon.net/论坛搜索和开发文档 比如这里我就在官方找到正确的id描述 然后我就把参数标签换过来
阅读更多...
【Python编程】Linux创建虚拟环境并配置与notebook相连接

【Python编程】Linux创建虚拟环境并配置与notebook相连接

1.创建 使用 venv 创建虚拟环境。例如，在当前目录下创建一个名为 myenv 的虚拟环境： python3 -m venv myenv 2.激活 激活虚拟环境使其成为当前终端会话的活动环境。运行： source myenv/bin/activate 3.与notebook连接 在虚拟环境中，使用 pip 安装 Jupyter 和 ipykernel： pip instal
阅读更多...
poj 1258 Agri-Net（最小生成树模板代码）

poj 1258 Agri-Net（最小生成树模板代码）

感觉用这题来当模板更适合。 题意就是给你邻接矩阵求最小生成树啦。~ prim代码：效率很高。172k...0ms。 #include<stdio.h>#include<algorithm>using namespace std;const int MaxN = 101;const int INF = 0x3f3f3f3f;int g[MaxN][MaxN];int n
阅读更多...
【机器学习】高斯过程的基本概念和应用领域以及在python中的实例

【机器学习】高斯过程的基本概念和应用领域以及在python中的实例

引言 高斯过程（Gaussian Process，简称GP）是一种概率模型，用于描述一组随机变量的联合概率分布，其中任何一个有限维度的子集都具有高斯分布 文章目录 引言一、高斯过程1.1 基本定义1.1.1 随机过程1.1.2 高斯分布 1.2 高斯过程的特性1.2.1 联合高斯性1.2.2 均值函数1.2.3 协方差函数（或核函数） 1.3 核函数1.4 高斯过程回归（Gauss
阅读更多...
【学习笔记】 陈强-机器学习-Python-Ch15 人工神经网络（1）sklearn

【学习笔记】 陈强-机器学习-Python-Ch15 人工神经网络（1）sklearn

系列文章目录 监督学习：参数方法 【学习笔记】 陈强-机器学习-Python-Ch4 线性回归 【学习笔记】 陈强-机器学习-Python-Ch5 逻辑回归 【课后题练习】 陈强-机器学习-Python-Ch5 逻辑回归（SAheart.csv） 【学习笔记】 陈强-机器学习-Python-Ch6 多项逻辑回归 【学习笔记 及 课后题练习】 陈强-机器学习-Python-Ch7 判别分析 【学
阅读更多...
计算机毕业设计 大学志愿填报系统 Java+SpringBoot+Vue 前后端分离 文档报告 代码讲解 安装调试

计算机毕业设计 大学志愿填报系统 Java+SpringBoot+Vue 前后端分离 文档报告 代码讲解 安装调试

🍊作者：计算机编程-吉哥 🍊简介：专业从事JavaWeb程序开发，微信小程序开发，定制化项目、 源码、代码讲解、文档撰写、ppt制作。做自己喜欢的事，生活就是快乐的。 🍊心愿：点赞 👍 收藏 ⭐评论 📝 🍅 文末获取源码联系 👇🏻 精彩专栏推荐订阅 👇🏻 不然下次找不到哟~Java毕业设计项目~热门选题推荐《1000套》 目录 1.技术选型 2.开发工具 3.功能
阅读更多...
nudepy，一个有趣的 Python 库！

nudepy，一个有趣的 Python 库！

更多资料获取 📚 个人网站：ipengtao.com 大家好，今天为大家分享一个有趣的 Python 库 - nudepy。 Github地址：https://github.com/hhatto/nude.py 在图像处理和计算机视觉应用中，检测图像中的不适当内容（例如裸露图像）是一个重要的任务。nudepy 是一个基于 Python 的库，专门用于检测图像中的不适当内容。该
阅读更多...
基于UE5和ROS2的激光雷达+深度RGBD相机小车的仿真指南(五)：Blender锥桶建模

基于UE5和ROS2的激光雷达+深度RGBD相机小车的仿真指南(五)：Blender锥桶建模

前言 本系列教程旨在使用UE5配置一个具备激光雷达+深度摄像机的仿真小车，并使用通过跨平台的方式进行ROS2和UE5仿真的通讯，达到小车自主导航的目的。本教程默认有ROS2导航及其gazebo仿真相关方面基础，Nav2相关的学习教程可以参考本人的其他博客Nav2代价地图实现和原理–Nav2源码解读之CostMap2D(上)-CSDN博客往期教程： 第一期：基于UE5和ROS2的激光雷达+深度RG
阅读更多...
代码随想录冲冲冲 Day39 动态规划Part7

代码随想录冲冲冲 Day39 动态规划Part7

198. 打家劫舍 dp数组的意义是在第i位的时候偷的最大钱数是多少 如果nums的size为0 总价值当然就是0 如果nums的size为1 总价值是nums[0] 遍历顺序就是从小到大遍历 之后是递推公式 对于dp[i]的最大价值来说有两种可能 1.偷第i个 那么最大价值就是dp[i-2]+nums[i] 2.不偷第i个 那么价值就是dp[i-1] 之后取这两个的最大值就是d
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2