2021年中国研究生数学建模竞赛E题
 
 
信号干扰下的超宽带（UWB）精确定位问题                                                       
 
 
一、背景
 
UWB（Ultra-Wideband）技术也被称之为“超宽带”，又称之为脉冲无线电技术。这是一种无需任何载波，通过发送纳秒级脉冲而完成数据传输的短距离范围内无线通信技术，并且信号传输过程中的功耗仅仅有几十µW。UWB因其独有的特点，使其在军事、物联网等各个领域都有着广阔的应用。其中，基于UWB的定位技术具备实时的室内外精确跟踪能力，定位精度高，可达到厘米级甚至毫米级定位。UWB在室内精确的定位将会对卫星导航起到一个极好的补充作用，可在军事及民用领域有广泛应用，比如：电力、医疗、化工行业、隧道施工、危险区域管控等。UWB更多应用场景请参见[4—6]。
 
UWB的定位技术有多种方法，本文仅考虑基于飞行时间（Time of Flight, TOF）的测距原理，它是UWB定位法中最常见的定位方法之一。TOF测距技术属于双向测距技术，其通过计算信号在两个模块的飞行时间，再乘以光速求出两个模块之间的距离，这个距离肯定有不同程度的误差，但其精度已经比较高。                                                                                                                              
 
在室内定位的应用中，UWB技术可以实现厘米级的定位精度（一般指2维平面定位）,并具有良好的抗多径干扰和衰弱的性能以及具有较强的穿透能力。但由于室内环境复杂多变UWB 通信信号极易受到遮挡，虽然UWB技术具有穿透能力，但仍然会产生误差，在较强干扰时，数据会发生异常波动（通常是时间延时），基本无法完成室内定位，甚至会造成严重事故。因此，信号干扰下的超宽带（UWB）精确定位问题成为亟待解决的问题。
 
二、问题描述
 
为解决信号干扰下的超宽带（UWB）精确定位问题，我们通过实际场景实测，采集到一定数量的数据，即利用UWB的定位技术（TOF），采集到锚点（anchor）与靶点（Tag）之间的距离，希望通过数学建模（或算法）方法 ，无论信号是否干扰，都可以给出目标物（靶点）的精确定位（3维坐标）。
 
三、实验场景和数据采集
 
如图所示，在5000mm*5000mm*3000mm的测试环境中，分别在4个角落A0，A1，A2，A3放置UWB锚点（anchor），锚点向所有方向发送信号。Tag是UWB标签（靶点），即需要定位的目标（只在测试环境范围内）。Tag接收到4个UWB锚点（anchor）的信号（无论信号是否干扰，Tag一般都可以接收到信号），利用TOF技术，分别解算出对应的4个距离数据。
 
实验在实验场景1中采集了Tag在324个不同位置，在信号无干扰和信号干扰下的UWB数据，即每个位置各测试（采集）2次，一次信号无干扰，另一次信号有干扰（锚点与靶点间有遮挡），注意：每次采集数据时，由于Tag在同一位置会停留一会儿时间，而锚点与Tag之间每0.2—0.3秒之间就会发送、接收信号一次，所以在同一位置点，UWB会采集到多组数据（多组数据都代表同一位置的信息），组数的多少视Tag在同一位置的时间而定，停留的时间越长，组数就越多。数据见文件夹“附件1：UWB数据集”。
 
 
 
图1   实测环境示意图
 
 
  实验场景1：
 
靶点（Tag）范围：5000mm*5000mm*3000mm
 
锚点（anchor）位置（单位：mm）：
 
A0（0，0，1300）、   A1（5000，0，1700）、
 
A2（0，5000，1700）、A3（5000，5000，1300）
 
四、数据文件说明
 
（1）UWB数据集
 
 “附件1：UWB数据集”有2个文件夹和1个文件，1个文件（Tag坐标信息.txt）存放324个不同位置的编号及3维坐标信息，2个文件夹中1个存放信号无干扰下（正常）采集的数据（各文件名为x.正常.txt，x表示对应的位置编号），另1个存放信号有干扰下（异常）采集的数据（各文件名为x.异常.txt，x表示对应的位置编号）。
 
（2）数据文件
 
Tag在每个位置都采集了2个数据文件（1个正常，另1个异常），共有648个数据文件，无论正常、异常数据，数据格式都一样，每个数据文件开头第1行为采集开始行，无实际意义，接下来，每4行为一组，表示UWB采集的一组完整数据（一组数据表示一个样品），如：
 
T:144235622:RR:0:0:950:950:118:1910
 
T:144235622:RR:0:1:2630:2630:118:1910
 
T:144235622:RR:0:2:5120:5120:118:1910
 
T:144235622:RR:0:3:5770:5770:118:1910
 
 
这4行数据的含义分别是:
 
Tag标识：时间戳：Range Report的缩写：Tag ID：锚点ID：该锚点的测距值(mm)：测距值的校验值：数据序列号：数据编号（每个数据之间用“：”分隔）。实际上就是提供了4个锚点到靶点（Tag）的距离，即
 
A0到靶点距离为：950mm 
 
A1到靶点距离为：2630mm 
 
A2到靶点距离为：5120mm
 
A3到靶点距离为：5770mm
 
每个数据文件都有多组数据，表示在同一位置连续时间内UWB自动采集到的多组数据。
 
五、完成任务
 
试根据上述数据，完成如下任务：
 
任务1：数据预处理（清洗）
 
无论是信号无干扰下采集数据，或信号干扰下采集数据，Tag在同一坐标点上都采集多组数据（见附件1中648个数据文件），请用某种方法把每个数据文件相应数值抓取出来，并转换成二维表（矩阵）形式（txt、Excel或其他数据格式），每一行代表一组数据（即一个样品），然后对这些数据文件进行预处理（清洗），删除掉一些“无用”（异常、缺失、相同或相似）的数据（样品）。经处理后，“正常数据”所有数据文件和“异常数据”所有数据文件最后各保留多少组（多少个样品）数据，并重点列出以下4个数据文件，经处理后保留的数据（矩阵形式）；
 
“正常数据”文件夹中：      24.正常.txt、  109.正常.txt
 
“异常数据”文件夹中：      1.异常.txt、   100.异常.txt
 
任务2: 定位模型
 
利用任务1处理后的数据，分别对“正常数据”和“异常数据”，设计合适的数学模型（或算法），估计（或预测）出Tag的精确位置，并说明你所建立的定位模型（或算法）的有效性；同时请利用你的定位模型（或算法）分别对附件2中提供的前5组（信号无干扰）数据和后5组（信号有干扰）数据进行精确定位（3维坐标）；
 
注意：（1）定位模型必须体现实验场景信息；
 
（2）请同时给出定位模型的3维（x，y，z）精度、2维（x，y）精度以及1维的各自精度。
 
任务3：不同场景应用
 
我们的训练数据仅采集于同一实验场景（实验场景1），但定位模型应该能够在不同实际场景上使用，我们希望你所建立的定位模型能够应用于不同场景。附件3中10组数据采集于下面实验场景2（前5组数据信号无干扰，后5组数据信号有干扰），请分别用上述建立的定位模型，对这10 组数据进行精确定位（3维坐标）；
 
实验场景2：
 
靶点（Tag）范围：5000mm*3000mm*3000mm
 
锚点（anchor）位置（单位：mm）：
 
A0（0，0，1200）、   A1（5000，0，1600）、
 
A2（0，3000，1600），A3（5000，3000，1200）
 
任务4: 分类模型
 
上述定位模型是在已知信号有、无干扰的条件下建立的，但UWB在采集数据时并不知道信号有无干扰，所以判断信号有无干扰是UWB精确定位问题的重点和难点。利用任务1处理后的数据，建立数学模型(或算法)，以便区分哪些数据是在信号无干扰下采集的数据，哪些数据是在信号干扰下采集的数据？并说明你所建立的分类模型(或算法)的有效性；同时请用你所建立的分类模型（或算法）判断附件4中提供的10组数据（这10组数据同样采集于实验场景1）是来自信号无干扰或信号干扰下采集的？
 
任务5: 运动轨迹定位
 
运动轨迹定位是UWB重要应用之一，利用静态点的定位模型，加上靶点自身运动规律，希望给出动态靶点的运动轨迹。附件5是对动态靶点采集的数据（一段时间内连续采集的多组数据），请注意，在采集这些数据时，会随机出现信号干扰，请对这个运动轨迹进行精确定位，最终画出这条运动轨迹图（数据采集来自实验场景1）。
 
参考文献
 
[1] Kok M, Hol J D, Schon T B. Indoor Positioning Using Ultrawideband and Inertial Measurements[J]. IEEE Transactions on Vehicular Technology, 2015, 64(4): 1293—1303.
 
[2] Arias-De-Reyna E. A Cooperative Localization Algorithm for UWB Indoor Sensor Networks[J]. Wireless Personal Communications, 2013, 72(1):85—99.
 
[3]缪希仁,范建威等，基站异常情况下基于改进极限学习机的超宽带室内定位方法，传感技术学报，2020，33（10）：1—10.
 
[4]到底什么是 UWB超宽带技术 有哪些应用场景？-电子发烧友网
 
[5]基于UWB的应用场景有哪些？UWB解决方案应用案例！-电子发烧友网
 
[6]全面解析UWB及其应用场景
 
 
 
 
 

（一）python爬虫验证码识别（去除干扰线）
 
1.开发环境与工具
 
python27：sklearn、pytesser、opencv等
pycharm
windows7
 
 
2. 数据集
 

 用request库爬虫抓取某一网站验证码1200张，并做好标注
 
3.验证码识别大概步骤
 
转化成灰度图
去背景噪声
图片分割
 
(1)转化成灰度图
 
im = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
 
(2)去除背景噪声
 
验证码去除干扰线的思想可参考链接：验证码去除干扰线
 本文所使用的去除背景噪声的方法：
      认真观察我们的实验数据，发现根据线降噪方法来去除噪声是不可行的，因为我们的图片干扰线很粗，和数字差不多粗。那再认真观察一下，1个数字的颜色都是一个‘“色”，那么是否可以跟据颜色来分呢？
      博主认真想了下，先转化成灰度图，再通过图像分割把图片分割一下，去除掉边框和部分噪声，这样就分成了4张图，然后统计每张图的灰度直方图（自己设置bins），找到第二大所对应的像素范围，即某一像素范围内像素数第二多所对应的像素范围（像素最多的应该是白色，空白处），取像素范围中位数mode，然后保留（mode±biases）的像素。这样就可以将大部分噪声去除掉啦。
 （这段描述有点复杂，需要一定的图像基础，不懂的，可以看代码del_noise()方法）
    
 
4. 实验方法
 
（1）使用Google的pytesser识别图片方法，准确率在50%左右
        pytesser环境配置参考 https://www.cnblogs.com/lcosima/p/7138091.html
 （2）使用机器学习方法KNN，准确率在99.6%
 
pip install sklearn
 
5.pytesser方法
 
这里im_cut是指：分割好验证码后，传入的子图片
 
def del_noise(im_cut):
    ''' variable：bins：灰度直方图bin的数目
                  num_gray:像素间隔
        method：1.找到灰度直方图中像素第二多所对应的像素，即second_max,因为图像空白处比较多所以第一多的应该是空白，第二多的才是我们想要的内容。
                2.计算mode
                3.除了在mode+-一定范围内的，全部变为空白。
    '''
    bins = 16
    num_gray = math.ceil(256 / bins)
    hist = cv2.calcHist([im_cut], [0], None, [bins], [0, 256])
    lists = []
    for i in range(len(hist)):
        # print hist[i][0]
        lists.append(hist[i][0])
    second_max = sorted(lists)[-2]
    bins_second_max = lists.index(second_max)

    mode = (bins_second_max + 0.5) * num_gray

    for i in range(len(im_cut)):
        for j in range(len(im_cut[0])):
            if im_cut[i][j] < mode - 15 or im_cut[i][j] > mode + 15:
                # print im_cut[i][j]
                im_cut[i][j] = 255
    return im_cut

# 替换文本
def replace_text(text):
    text = text.strip()
    text = text.upper()
    rep = {'O': '0',
           'I': '1',
           'L': '1',
           'Z': '7',
           'A': '4',
           '&': '4',
           'S': '8',
           'Q': '0',
           'T': '7',
           'Y': '7',
           '}': '7',
           'J': '7',
           'F': '7',
           'E': '6',
           ']': '0',
           '?': '7',
           'B': '8',
           '@': '6',
           'G': '0',
           'H': '3',
           '$': '3',
           'C': '0',
           '(': '0',
           '[': '5',
           'X': '7',
           '`': '',
           '\\': '',
           ' ': '',
           '\n': '',
           '-': '',
           '+': '',
           '*': '',
           '.': '',
           ';': ''
           }

    #判断是否有数字，有数字直接返回第一个数字，不需要字符替换
    print text
    if len(text) >= 1:
        pattern = re.compile(u'\d{1}')
        result = pattern.findall(text)
        if len(result) >= 1:
            text = result[0]
        else:
            # 字符替换,替换之后抽取数字返回
            for r in rep:
                text = text.replace(r, rep[r])
            pattern = re.compile(u'\d{1}')
            result = pattern.findall(text)
            if len(result) >= 1:
                text = result[0]

    return text

 
主方法
 
#im_cut = [im_cut_1, im_cut_2, im_cut_3, im_cut_4]
for i in range(4):
	im_temp = del_noise(im_cut[i])
	im_result = Image.fromarray(im_temp.astype('uint8'))
	#使用pytesser识别
	text = image_to_string(im_result)
	#做文本替换处理
	text_rep = replace_text(text)
	#获得预测结果
	pre_text.append(text_rep)
	pre_text = ''.join(pre_text)
 
结果
 
 
6.KNN分类（sklearn）
 
（1）先转成灰度图，去背景噪声，分割1200张已标注好的图片，得到4800张子图片；
 
 （2）用knn训练分类器，训练集：测试集=0.8，训练结果精度达到99%以上；
 （3）使用训练好的模型，进行实际验证码预测，效果不错。
 
# -*-coding:utf-8-*-
import numpy as np
from sklearn import neighbors
import os
from sklearn.preprocessing import LabelBinarizer
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report
from sklearn.externals import joblib

import cv2

if __name__ == '__main__':
    # 读入数据
    data = []
    labels = []
    img_dir = './img_train_cut'
    img_name = os.listdir(img_dir)
    # number = ['0','1', '2','3','4','5','6','7','8','9']
    for i in range(len(img_name)):
        path = os.path.join(img_dir, img_name[i])
        # cv2读进来的图片是RGB3维的，转成灰度图，将图片转化成1维
        image = cv2.imread(path)
        im = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
        image = im.reshape(-1)
        data.append(image)
        y_temp = img_name[i][-5]
        labels.append(y_temp)

    # 标签规范化
    y = LabelBinarizer().fit_transform(labels)

    x = np.array(data)
    y = np.array(y)

    # 拆分训练数据与测试数据
    x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)

    # 训练KNN分类器
    clf = neighbors.KNeighborsClassifier()
    clf.fit(x_train, y_train)

    # 保存分类器模型
    joblib.dump(clf, './knn.pkl')

    # # 测试结果打印
    pre_y_train = clf.predict(x_train)
    pre_y_test = clf.predict(x_test)
    class_name = ['class0', 'class1', 'class2', 'class3', 'class4', 'class5', 'class6', 'class7', 'class8', 'class9']
    print classification_report(y_train, pre_y_train, target_names=class_name)
    print classification_report(y_test, pre_y_test, target_names=class_name)

    # clf = joblib.load('knn.pkl')
    # pre_y_test = clf.predict(x)
    # print pre_y_test
    # print classification_report(y, pre_y_test, target_names=class_name)
 
结果截图：
 
 
ps:在运行过程中，程序还可能会帮你检测出人工标错的哟！
 
 
 
代码已上传至https://download.csdn.net/download/weixin_40267472/11097627。
 
如果喜欢这篇博客，可以小小打赏下小编，表示下支持哟！！！(毕竟好不容易整理出来的) 谢谢哟~~
 
链接：https://pan.baidu.com/s/12ajeAEe_hlH0j-wklvWukg
 提取码：ya85
 

 
基于MATLAB的直接序列扩频通信系统仿真
 
前言
文章迭代更新
  
仿真流程图
关键技术细节
扩频解扰解扩部分
加扰去扰部分
加扰原理图
   
调制解调部分
调制原理图
解调原理图
   
高斯信道部分
  
实验结果
扩频增益为10时，walsh矩阵为64阶时，不同振幅下信噪比误码率曲线
扩频增益为20时，walsh矩阵为64阶时，不同振幅下信噪比误码率曲线
误码率随信噪比的变化曲线
  
源代码
可以修改的参数
！下载链接！
核心代码
主函数
walsh码生成
M序列生成
扩频模块
解扩模块
加扰模块
去扰模块
调制模块
解调模块
   
模块测试代码
扩频解扩模块测试
加扰去扰模块测试
调制解调模块测试
   
非核心代码
八进制转二进制模块
点乘模块
码元比较模块
双极性码生成模块
数字位数获取模块
功率计算模块
周期延拓模块
三极码生成模块
0~1转双极性模块
分组相加模块
  
 

 
前言
 
直接扩频序列调制是用速率很高的伪噪声码序列与信息码序列模二相加（波形相乘）后得到复合码序列，用复合码序列去控制载波相位，从而获得直接扩频序列信号的。直接扩频通信具有低截获概率、抗干扰能力强以及易于实现码分多址等优点，在抗干扰通信及民用移动通信中都得到了广泛的应用。
 
文章迭代更新
 
文章进行了大改
 发现了为什么振幅是0.2和0.4时，系统无法正常工作的原因。因为bitMultiple函数把运算结果格式转换成了int8，导致精度大量失真！当时转换格式是为内存空间与运行速度做打算，结果今天发现做了负优化！修改后，发现载波振幅对于系统的误码率曲线几乎无影响。
 另外，解调函数也做了修改，会根据输入自动计算判决阈值，解决了之前人为设定阈值
 的局限，阈值采用正态分布解算法，详见demodulate函数内注释。
 新版修改了以下函数：
 1. bitMultiple
 2. demodulate
 3. 新增arrayGroupSum函数
 
仿真流程图
 

 
 
关键技术细节
 
仿真中，用户码元使用双极性(1和-1)码。
代码的最后通过还原信号与码元做对比，计算误码率，来评价通信质量的好坏。
为了研究信噪比对误码率的影响，代码中用了多线程，针对不同的信噪比进行解调去扰解扩,计算不同信噪比下的误码率。
为了研究载波振幅对信噪比——误码率曲线的影响，增加了振幅的迭代。
 
扩频解扰解扩部分
 
采用64阶的walsh码矩阵，由walsh函数产生。
两个用户使用同一扩频矩阵的不同相位(实际上是不同一行)进行扩频。
扩频时，首先对用户码元按照扩频增益，进行周期延拓，然后点乘扩频码，实现扩频。
解扩时，首先对输入的去扰后的码按扩频时一样的相位对应点乘扩频码，然后进行判决。
 
加扰去扰部分
 
加扰采用的是5阶的M序列，反馈系数为67(八进制)。
加扰同样采用点乘，输入码元1:1地点乘加扰码。
因为加扰码具有很高的自相关性，所以去扰时只需执行加扰时一样的操作即可。
 
加扰原理图
 
 
调制解调部分
 
调制采用BPSK(二进制移相键控)调制。
调制所使用的载波为正弦波，通过对其等间隔采样形成离散样值。
调制时，一个加扰后的码元对应相乘一串正弦波的周期采样值，把相乘结果汇总后就是调制后的结果
解调时，把一个码元所对应的一串正弦波的周期采样值与载波采样值进行点乘，根据点乘结果中，正数、负数还是零多，判决这是什么码元。
 
调制原理图
 
 
解调原理图
 
 
高斯信道部分
 
使用matlab自带的函数awgn(input,snr(dB),inputPower)。
只要传入需要加高斯噪声的信号，信噪比还有输入信号的功率，即可返回添加了高斯白噪声的信号。
 
实验结果
 
扩频增益为10时，walsh矩阵为64阶时，不同振幅下信噪比误码率曲线
 
 
扩频增益为20时，walsh矩阵为64阶时，不同振幅下信噪比误码率曲线
 

 以上两次实验结果显示，载波振幅对于系统性能无明显影响，这是文章改版前的一个很大的错误！在此笔者向各位读者道歉。
 
误码率随信噪比的变化曲线
 

 观察以上曲线可以发现，随着信噪比的上升，误码率以近似于抛物线的形式在不断下降，并最终等于0，这足以证明本次直接序列扩频系统仿真符合预期结果。
 
源代码
 
可以修改的参数
 
轮回次数。通过修改轮回次数，获得更好的曲线，或者更快的代码运行速度，代码运行参考速度：i7-4790:实测速率9.8秒/轮回 @1280*2&1280码元20&10扩频增益。
用户1的码元数量与扩频增益、用户2的码元数量与扩频增益。修改时注意两用户的码元必须是walsh码阶数的增数倍，同时要保证：用户1的码元数量 * 用户1的扩频增益=用户2的码元数量 * 用户2的扩频增益。
walsh码的阶数。但walsh码阶数必须大于扩频增益。
用户1扩频相位、用户2扩频相位。注意相位必须在1到64之间取值，两个用户的取值不能一样。
M序列阶数与返回系数。反馈系数取值有相关要求，详见百度，作为参数传入时要以八进制表示。
半个载波周期的采样点数。
最小最大信噪比，以及其尝试步进。
最小最大载波振幅，以及其尝试步进。
 
！下载链接！
 
 
 
链接：https://pan.baidu.com/s/14KNoWlD7UyC9WAEza63W5A
 提取码：qg1y
 
 
核心代码
 
主函数
 
main.m
 
%{
    本函数是整个仿真的主函数，用于研究在移动通信时，信噪比对误码率的影响
首先生成随机双极性码，然后经过扩频，加扰，BPSK调制，加高斯白噪声，混合
然后模拟接收端的解调，去扰，解扩，判决。
    通过比较接收端判决输出与原来的码元，计算出误码率，首先通信系统的仿真
以及误码率-信噪比的变化曲线的绘制.
2019/11/22(以上)
    在当前版本中，添加了对信噪比&误码率曲线随BPSK调制载波振幅变化而变化
的相关研究功能。因为实际发现，之前的代码中曲线会在-20dB到0dB处出现平缓现象，
实际排查发现是载波的振幅导致的。振幅研究结果发现：振幅小于0.4时，系统无法正常
工作！其误码率曲线在后期随着信噪比的增加呈现反常膨胀！
^-^想要得到抛物线，只有在载波振幅为1的时候。^-^,不然随着振幅的增加，曲线曲线趋于平缓的信噪比范围会越大；
同时发现大于0.6的曲线会近似相交于某一点，交点前同等信噪比下振幅高的误码率低，但在交点后
***同等的信噪比下振幅小的误码率反而低！！！***
    在本版本中，walsh矩阵的所有码元都被用于扩频，没有使用前一版的矩阵截取方法，
规避了鬼魅版的非严格正交问题，但前面的截取版本恰恰说明非严格正交对本系统影响曲线
的影响并不是很大，但却可以使得系统能在更低的误码率下工作。
    至此，代码已经经历10次版本迭代。
2019/12/1(以上)
    第二次答辩，仍出现部分平滑问题，并且老师说系统性能太好了，与实际系统不符，并且
曲线不应该受振幅的影响，基于此再对代码做修改，计算了输入信号的平均功率，作为第3
个参数传给awgn，否则awgn会把输入信号的功率视为0dBW！现在曲线基本重合。但0.2与0.4
的问题仍未解决！
    至此，程序第11次版本迭代
2019/12/2(以上)
    第十二次迭代。
    发现了为什么振幅是0.2和0.4时，系统无法正常工作的原因，是因为
bitMultiple函数把运算结果格式转换成了int8，导致精度大量失真！当时转换格式是为
内存空间与运行速度做打算，结果今天发现做了负优化！修改后，发现载波振幅对于系统
的误码率曲线几乎无影响。
    另外，解调函数也做了修改，会根据输入自动计算判决阈值，解决了之前人为设定阈值
的局限，阈值采用正态分布解算法，详见demodulate函数内注释。
2020/7/1(以上)
    i7-4790:实测速率9.8秒/轮回 @1280*2&1280码元20&10扩频增益
%}
clear variable;
close all;

mulTimes = 5;  %轮回次数

walshOrder = 64;    %walsh码的阶数，必须大于扩频增益

N1 = 1280*2;   %用户1码元数量
N2 = 1280;   %用户2码元数量

user1SPgain = 10;  %用户1的扩频增益
user2SPgain = 20;  %用户2的扩频增益

%记录两个用户的walsh相位,注意相位必须在1到64之间取值
%两个用户的取值不能一样
user1Phase = 2; %用户1的相位
user2Phase = 16; %用户2的相位

%加扰使用的m序列的参数
mOrder = 5; %级数5级
feedBack = 67;%反馈系数67

%调整时，半个周期的采样点数
samplePiont = 4;

%生成需要使用的walsh码
walshCode = walsh(walshOrder);

%针对低性能电脑做优化，采取以时间换取空间的思路
maxSnr = 20;    %尝试的最大信噪比
minSnr = -20;   %尝试的最小信噪比
div = 1;      %信噪比的尝试步进
maxTime = (maxSnr-minSnr)/div;  %尝试次数
timesUser1Acc = zeros(mulTimes,maxTime);
timesUser2Acc = zeros(mulTimes,maxTime);

%生成双极性码片
user1 = genBipolar(N1);
user2 = genBipolar(N2);   

%扩频
spread1 = spreadSpectrum(user1,walshCode,user1SPgain,user1Phase);
spread2 = spreadSpectrum(user2,walshCode,user2SPgain,user2Phase);

%加扰
Mseq1 = MseqGen(mOrder,feedBack); %用户1加扰用的m序列
Mseq2 = MseqGen(mOrder,feedBack); %用户2加扰用的m序列
user1scarm = scarmbling(spread1,Mseq1);
user2scarm = scarmbling(spread2,Mseq2);

maxAmp = 1.2;    %尝试的最大载波振幅
minAmp = 0.2;   %尝试的最小载波振幅
divAmp = 0.2;      %载波振幅的尝试步进
maxTimesAmp = floor((maxAmp-minAmp)/divAmp);  %振幅尝试次数

ampRecords1 = zeros(maxTimesAmp,maxTime);
ampRecords2 = zeros(maxTimesAmp,maxTime);

for amp = 1:maxTimesAmp %测试不同的载波振幅下的曲线
    
    fprintf('目前正在第%d个振幅\n',amp);

    %调制BPSK
    %生成载波,两用户使用同一个载波
    carrier = (minAmp + divAmp*(amp-1))*sin(0:(pi/samplePiont):(2*pi-2*pi/samplePiont));
    user1modulate = modulate(user1scarm,carrier);
    user2modulate = modulate(user2scarm,carrier);
    
    %计算载波功率
    power = powerCnt(carrier);
    
    for times = 1:mulTimes

        fprintf('目前正在第%d个轮回\n',times);

        user1Acc = zeros(1,maxTime);
        user2Acc = zeros(1,maxTime);

        parfor index = 1:maxTime 
            snr = minSnr + (index-1)*div; %加在发送信号上的高斯噪声的信噪比(dBW)
                        
            %通过高斯信道，添加高斯噪声
            user1send = awgn(user1modulate,snr,power);
            user2send = awgn(user2modulate,snr,power);

            %接收方
            receive = user1send + user2send; %收到混合起来的信号

            %解调
            demodulateRes = demodulate(receive,carrier);

            %去扰
            user1Descarm = deScarmbling(demodulateRes,Mseq1);
            user2Descarm = deScarmbling(demodulateRes,Mseq2);

            %解扩
            user1deDS = deSpreadSpectrum(user1Descarm,walshCode,user1SPgain,user1Phase);
            user2deDS = deSpreadSpectrum(user2Descarm,walshCode,user2SPgain,user2Phase);

            %计算误码率
            [~,user1Accuracy] = compare(user1,user1deDS);
            [~,user2Accuracy] = compare(user2,user2deDS);
            user1Acc(index) = 1-user1Accuracy;
            user2Acc(index) = 1-user2Accuracy;
        end
        timesUser1Acc(times,:) = user1Acc;
        timesUser2Acc(times,:) = user2Acc;
    end
    %总结统计多次实验的结果
    for i = 1:maxTime
        user1Acc(i) = mean(timesUser1Acc(:,i));
        user2Acc(i) = mean(timesUser2Acc(:,i));
    end
    ampRecords1(amp,:) = user1Acc;
    ampRecords2(amp,:) = user2Acc;
end

%误码率随信噪比的变化曲线
figure(1);
X1 = (minSnr:div:maxSnr-div);
semilogy(X1,ampRecords1(5,:),'b');
xlabel('信噪比(dB)');
ylabel('误码率');
title('误码率随信噪比的变化曲线');
hold on;
semilogy(X1,ampRecords2(5,:),'g');
legend('用户1(扩频增益10)','用户2(扩频增益20)');

%用户1振幅不同的误码率随信噪比的变化曲线
figure(2);
for i = 1:maxTimesAmp
    semilogy(X1,ampRecords1(i,:));
    hold on;
end
xlabel('信噪比(dB)');
ylabel('误码率');
title('用户1振幅不同的误码率随信噪比的变化曲线');
legend('0.2','0.4','0.6','0.8','1.0');

%用户2振幅不同的误码率随信噪比的变化曲线
figure(3);
for i = 1:maxTimesAmp
    semilogy(X1,ampRecords2(i,:));
    hold on;
end
xlabel('信噪比(dB)');
ylabel('误码率');
title('用户2振幅不同的误码率随信噪比的变化曲线');
legend('0.2','0.4','0.6','0.8','1.0');
 
walsh码生成
 
walsh.m
 
% 产生 Walsh函数通用函数
% 参数N表示Walsh函数阶数,当N不是2的幂时,通过向无穷大取整使得所得Walsh阶数为2的幂
function [walsh]=walsh(N)
    M = ceil(log2(N));
    wn = 0;
    for i = 1:M
        w2n = [wn,wn;wn,~wn];
        wn = w2n;
    end
    wn(wn == 0) = -1;
    walsh = int8(wn);
end
 
M序列生成
 
MseqGen.m
 
%M序列发生器
%order:阶数
%setBakc:反馈系数
function res = MseqGen(order,feedBack)
    res = zeros(2^order-1,order);
    feedBack = Oct2Bin(feedBack);
    feedBack = feedBack(2:6);
    temp = rand(1,order);
    temp(temp < 0.5) = 0;
    temp(temp >= 0.5) = 1;
    res(1,:) = temp;
    for i = 2:2^order-1
        newBit = sum(res(i-1,feedBack == 1));
        newBit = mod(newBit,2);
        for j = 1:order-1
            res(i,j+1) = res(i-1,j);
        end
        res(i,1) = newBit;
    end
    res(res == 0) = -1;
end
 
扩频模块
 
spreadSpectrum.m
 
%扩频函数
%userCode:需要扩频的用户码元
%PNseq:用于扩频的随机码
%gain:扩频增益
%phase:用户扩频码相位
function res = spreadSpectrum(userCode,PNseq,gain,phase)
    %首先对码元进行周期延拓
    [~,userCode2] = selfCopy(userCode,gain);
    %计算扩频码的行数
    [lineSize,~] = size(PNseq);
    %对扩频码进行重排行序,使初相位位于第一行
    PN = PNseq(phase:lineSize,:);
    PN = [PN;PNseq(1:phase-1,:)];
    res = bitMultiple(userCode2,PN(:)');
end
 
解扩模块
 
deSpreadSpectrum.m
 
%本函数用于解扩
%userCode:需要解扩的用户码元
%PNseq:用于扩频的随机码
%gain:扩频增益
%phase:用户扩频码相位
function res = deSpreadSpectrum(userCode,PNseq,gain,phase)
    [lineSize,~] = size(PNseq);
    PN = PNseq(phase:lineSize,:);
    PN = [PN;PNseq(1:phase-1,:)];
    temps = bitMultiple(userCode,PN(:)');
    %这里是matlab的一个小bug,重排序以列作为索引,我的要求是以行
    %作为索引,所以要行列反写再取转置矩阵
    temps = reshape(temps,gain,length(temps)/gain);
    temps = temps';
    %解扩第二步,码元判决
    res = ones(1,length(temps(:))/gain);
    for i = 1:length(temps(:))/gain
        if sum(temps(i,:)) < 0
            res(i) = -1;
        end
    end
end
 
加扰模块
 
scarmbling.m
 
%加扰函数
%source:被加扰的信号
%PNCode:用于加扰的扰码
%按1:1的数据加扰
function res = scarmbling(source,PNCode)
    res = bitMultiple(source,PNCode);
end
 
去扰模块
 
deScarmbling.m
 
%本函数用于对信号做去扰处理,实际操作与加扰完全一致
%input:需要去扰的信号
%PNseq:用于去扰的随机序列
function res = deScarmbling(input,PNseq)
    res = bitMultiple(input,PNseq);
end
 
调制模块
 
modulate.m
 
function res = modulate(source,carrier)
%调制函数
%source:被调制信号
%carrier:载波信号
    sizeSource = length(source);
    sizeCarrier = length(carrier);
    res = zeros(sizeSource,sizeCarrier);
    for i = 1:sizeSource
        res(i,:) = double(source(i))*carrier;
    end
    %把res按行拼接成一行
    res = res';
    res = res(:);
end
 
解调模块
 
demodulate.m
 
function res = demodulate(input,carrier)
%本函数实现对接收信号的解调
%原理:极性比较
%input:被解调的信号
%carrier:载波信号

%按位循环相乘
res_bitMultiple = bitMultiple(input,carrier);
%分组相加
res_arrayGroupSum = arrayGroupSum(res_bitMultiple,length(carrier));
%{
    以上处理后，结果正负侧近似于泊松分布，为了处理方便，根据大数定理，
    近似地视为正态分布，从而计算阈值。按照正态分布，-1,1,0平分概率，每一个分得1/3的概率。
%}
%计算标准差
res_std = std(res_arrayGroupSum);
%计算平均数
res_mean = mean(res_arrayGroupSum);
%解算阈值
syms temp;
%以下计算时会有一个"Unable to solve symbolically. Returning a numeric solution using
%vpasolve"的警告,让系统不显示
warning off;
threshhold_negative = double(solve(normcdf(temp,res_mean,res_std)==1/3));
%计算概率为2/3时的阈值，也就是判决为1的阈值
threshhold_postive = double(solve(normcdf(temp,res_mean,res_std)==2/3));
res = zeros(1,length(res_arrayGroupSum));
for i = 1:length(res_arrayGroupSum)
    if res_arrayGroupSum(i) > threshhold_postive
        res(i) = 1;
    elseif res_arrayGroupSum(i) < threshhold_negative
        res(i) = -1;
    end
end
end
 
模块测试代码
 
本部分的代码用于单独测试三大模块是否正常工作
 
扩频解扩模块测试
 
testSpreadSpectrum.m
 
%本代码用于测试扩频，解扩模块是否正常工作
%注意扩频与解扩针对的是双极性码
function testSpreadSpectrum()
    raw = genBipolar(64);
    walshCode = walsh(64);
    afterDSSS = spreadSpectrum(raw,walshCode,4,2);
    res= deSpreadSpectrum(afterDSSS,walshCode,4,2);
    [trueNum,accuracy] = compare(raw,res);
    fprintf("正确码元数量:%d\n正确率:%f\n",trueNum,accuracy); 
end
 
加扰去扰模块测试
 
testScarmbling.m
 
%本代码用于测试加扰，去扰模块能否正常工作
function testScarmbling()
    %M序列发生器可以自定义阶数和反馈系数
    Mseq = MseqGen(5,67); %产生加扰用的m序列
    raw = tripleGen(1e6);
    afterScarmb = scarmbling(raw,Mseq);
    res = deScarmbling(afterScarmb,Mseq);
    [trueNum,accuracy] = compare(raw,res);
    fprintf("正确码元数量:%d\n正确率:%f\n",trueNum,accuracy);    
end
 
调制解调模块测试
 
testModulate.m
 
%本代码用于测试BPSK调制与解调模块能否正常工作
function testModulate()
    carrier = 10*sin(0:pi/32:2*pi-pi/32);
    raw = tripleGen(1e6);
    afterModu = modulate(raw,carrier);
    res = demodulate(afterModu,carrier);
    [trueNum,accuracy] = compare(raw,res);
    fprintf("正确码元数量:%d\n正确率:%f\n",trueNum,accuracy);   
end
 
非核心代码
 
八进制转二进制模块
 
Oct2Bin.m
 
%八进制转二进制,返回一个二进制数组
function res = Oct2Bin(value)
    [bitNum,bit] = getPalces(value);
    res = int8(zeros(1,bitNum*3));
    for i = 1:bitNum*3
        bitMain = bit(floor((i-1)/3)+1);
        switch mod(i-1,3)
            case 0
                if bitMain >= 4
                    res(i) = 1;
                end
            case 1
                if bitMain == 2 || bitMain == 3 || bitMain == 6 || bitMain == 7
                   res(i) = 1;
                end
            case 2
                if mod(bitMain,2) ~= 0
                    res(i) = 1;
                end
        end
    end
end
 
点乘模块
 
bitMultiple.m
 
%按bit循环相乘函数,如果信号长度不一致,则较短的信号被循环使用
%signal1:相乘信号1
%signal2:相乘信号2
%res:相乘结果
function res = bitMultiple(signal_1,signal_2)
    sizeSource = length(signal_1);
    sizePNCode = length(signal_2);
    %如果长度不满足条件时,调换顺序递归调用
    if sizeSource < sizePNCode
        res = bitMultiple(signal_2,signal_1);
        return;
    end
    res = zeros(1,sizeSource);
    for i = 1:sizeSource
        res(i) = signal_1(i) * signal_2(mod(i-1,sizePNCode)+1);
    end
    res = int8(res);
end
 
码元比较模块
 
compare.m
 
%比较两个数组
%input1:数组1
%input2:数组2
%res:正确码元数量
%accuracy:正确率
function [res,accuracy] = compare(input1,input2)
    temp = (input1 == input2);
    res = length(find(temp == 1));
    sizeInput = max(length(input1),length(input2));
    accuracy = res/double(sizeInput);
end
 
双极性码生成模块
 
genBipolar.m
 
%产生双极性码
%num:双极性码的规模
%res:满载双极性码的数组
function res = genBipolar(num)
    res = rand(1,num);
    res = value2Bipolar(res);
end
 
数字位数获取模块
 
getPlaces.m
 
%统计一个数字的位数,并返回每位的数字
function [Places,item] = getPalces(value)
    Places = 0;
    temp = abs(value);
    while temp ~= 0
        temp = floor(temp/10);
        Places = Places + 1;
    end
    item = uint8(zeros(1,Places));
    temp = value;
    for i = Places:-1:1
        item(i) = mod(temp,10);
        temp = floor(temp/10);
    end
end
 
功率计算模块
 
powerCnt.m
 
%计算输入信号的平均功率(dBW)
function res = powerCnt(input)
    res = 10*log10(sum(input.*input));
end
 
周期延拓模块
 
selfCopy.m
 
%本函数实现数组的周期延拓
%input:需要延拓的数组
%times-1:延拓的次数
%比如selfCopy([1,2,3],2) = [1,2,3,1,2,3,1,2,3]
function [res,res2] = selfCopy(input,times)
    if times <= 1
        res = input;
    else
        res = zeros(length(input),times);
        for i = 1:times
            res(:,i) = input;
        end
        res2 = res';
        res2 = res2(:)';
        res = res(:)';
    end
end
 
三极码生成模块
 
tripleGen.m
 
%本函数用于生成指定数量的随机三极性码,主要用于做模块测试
%num:数量
%res:满载结果的数组
function res = tripleGen(num)
    raw = rand(1,num);
    for i = 1:num
        if raw(i) >= 0.75
            raw(i) = 1;
        elseif raw(i) <= 0.25
            raw(i) = -1;
        else
            raw(i) = 0;
        end
    end
    res = int8(raw);
end
 
0~1转双极性模块
 
value2Bipolar.m
 
% 本函数用户把随机生成的0~1之间的double数转换成双极性码
function res = value2Bipolar(input)
    res = zeros(1,length(input));
    for index = 1:length(input)
        if(input(index) < 0.5)
            res(index) = -1;
        else
            res(index) = 1;
        end
    end
    res = int8(res);
end
 
分组相加模块
 
arrayGroupSum.m
 
function res = arrayGroupSum(input,group_num)
%arrayGroupSum 数组分组求和,如果数组长度除去group_num无法整除，则截短input数组
%input 目标数组
%group_num 分组数
len = floor(length(input)/group_num);
res = zeros(1,len);
for i = 1:len
    res(i) = sum(input((i-1)*group_num+1:i*group_num));
end
end

https://blog.csdn.net/xiaoyuanwuhui/article/details/88411896
 
用STM32单片机成功驱动一个超声波模块后，接下来便有了疑问如何用stm32单片机驱动多个超声波模块呢？（驱动一个超声波模块可以参考：stm32驱动超声波模块测距）
 超声波模块型号为HC-SR04，模块驱动方式选择利用外部中断的方式驱动（具体细节可参见上面链接）。
 经过一个下午的思考实验发现了两种驱动多个超声波模块的方法：
 
分时测量
 同时测量
 1.分时测量：
 通过类似于实现流水灯的方法，多个超声波模块依次测量并将测量结果依次输出。
 优点：各个超声波模块测量过程中不会出现相互干扰的情况。
 缺点：整体测量周期长，数据采样频率比较低。
 这里有几个问题点需要注意：
 
一个超声波模块需要在测量完成后再开始下一个模块的测量
 超声波模块在测量的过程中有时会跳动明显超出量程的值
 在带电拔掉超声波模块可能出现程序异常
 为解决上面几个问题，需要根据超声波模块的量程（2cm-400cm）进行限辐，这里我根据需要将超声波模块的最远距离设定为255cm。这样最远距离的测量时间为15ms（255cm x 2 / 100 / 340m/s * 1000 = 15ms），这样便可以将定时器计时超过15ms时直接按15ms计算。这样便基本可以解决上面的问题1和问题2，一个模块测量完成的时间不会超过15ms（实际预留时间时建议在此基础上多几毫秒），这样每过一个15ms便可以完成一个模块的测距。问题3可以通过每次测量开始前将超声波数据的初始值初始化为量程的最小值2cm，并将超声波模块的ECHO引脚设为低电平来解决（在上电状态下拔掉超声波模块，ECHO引脚会呈现高电平状态导致一直进外部中断）。这样如果串口上位机距离值一直打印为2cm便很明显知道该超声波模块没有接。
 2.同时测量：
 创建一个系统时钟，每一个超声波模块在测量过程中只记录时间戳，不直接占用定时器，这样便可以实现多个超声波模块同时测量。
 优点：整体测量周期短，数据采样频率高。
 缺点：各个超声波模块在测量过程中可能出现相互干扰情况。
 使用建议：将各个模块尽量大角度朝不同的方向。
 
同步多个超声波传感器
 
超声波传感器是用于物体检测或物位测量的非接触式装置。它们的工作原理是声音通过介质传播，其中换能器以特定频率发出声波。基于反射声脉冲到达换能器所花费的时间，传感器的内部逻辑确定物体是否在标称感测范围内并产生输出信号。
 
传感器发出的声波遵循锥形图案。如果相同换能器频率的多个传感器彼此相邻地操作，则一个传感器的发射脉冲可以被另一个传感器检测为反射脉冲。根据应用程序的不同，此方案可能并不理想。同时发射相同频率的超声波的两个换能器可以产生错误输出。
 
4716259-7408ee48e4895024.jpg
 4716259-4d4207be44182fbe.jpg
 为了防止这种干扰或串扰，简单的布线连接可以同步发射/接收时序模式。例如，UC2000-L2-I-V15超声波传感器上的V15（5针公头）连接器具有一条专用线，专用于同步多个传感器。要手动同步多个UC2000-L2-I-V15传感器，只需将每个传感器的同步输入线互连即可。
 
4716259-628f751de2a5790b.jpg
 由于在连接同步线时自动设定定时，因此一个换能器发出的超声波不会被另一个误接收为接收波。即使多个传感器可以通过同步线互连，每个连接在一起的传感器仍然具有不同的输出。使用同步线的唯一目的是交替每个传感器脉冲的定时。
 
注：当安装超声波传感器，它可能无法坚持的最小间隔距离。Pepperl + Fuchs提供具有同步输入的模型。这可以防止传感器串扰，并可以减少最小间隔距离。
 
使用同步线的缺点是增加了对象移动的响应时间。由于多个超声波传感器交替脉冲，每个传感器以串联模式脉冲。随着更多传感器连接在一起，输出生成将需要更长时间，因为组中的每个传感器必须在下一个传感器跟随之前自行发射和接收。
 
对于某些应用，不能选择布线连接和顺序选通模式。例如，在自动引导车辆（AGV）的碰撞避免中，传感器仅需要检测机器人前方是否存在物体。所有相同频率的超声波传感器可以间隔很远，使得它们不会交叉通话，但是一个传感器的发射脉冲可以被另一个传感器检测为反射脉冲，以可靠地检测AGV前面的当前物体。
 
可以使用接口在某些超声波传感器上编程这些参数; 例如，某些软件允许调整交替脉冲的定时。如果一个阵列中的所有传感器安装得太紧并且同时发生脉冲，则可能发生串扰。响应时间最短，但可以生成错误输出。以顺序模式对所有传感器进行编程大大增加了响应时间，因此通过接口进行编程可以允许阵列上的传感器组根据需要发射/接收。
 
例如，在阵列上订购1-2-3-4-5-6的六个传感器可以分别以1-3-5和2-4-6的组脉冲，以显着缩短响应时间，与脉冲序列相比1-2-3-4-5-6。
 
4716259-7b7aa067a9079478.jpg
 当声锥宽度较窄或传感器依靠彼此的输出来调节定时时，一些超声波传感器应用效果更好。参数调整取决于应用，具体时间是优化响应时间的关键。
 
尽管通过布线连接或接口参数编程使两个或更多个超声波传感器同步，但配置中的每个传感器的输出仍将彼此不同。两种同步方法都有其优点和缺点，但同时使用多个超声波传感器的关键是防止相同换能器频率的传感器之间的干扰。
 
其中有一个小插曲，就是当吧写好的程序烧进去之后，运行时总是出现每次返回一个同样的比正常值小的多的数据，比如说0.034cm，这明显是一个错误的数据，但是刚开始的时候，不知道为什么
 
总是这样，多次复位从新上电总是这一个数据。让我很是苦恼。但是幸运的是，在这样的情况中间，他又会有时出现一两个正常的的数据，让你有点摸不着头脑。
 
上网查了一下才慢慢明白，这种现象叫做“余震”，网上关于余震的解释大致有三种：
 
1、探头的余震。即使是分体式的，发射头工作完后还会继续震一会，这是物理效应，也就是余震。这个余震信号也会向外传播。如果你的设计是发射完毕后立刻切换为接收状态（无盲区），那么这个余震波会通过壳体和周围的空气，直接到达接收头、干扰了检测（注：通常的测距设计里，发射头和接收头的距离很近，在这么短的距离里超声波的检测角度是很大的，可达180度）。
 2、壳体的余震。就像敲钟一样，能量仍来自发射头。发射结束后，壳体的余震会直接传导到接收头，当然这个时间很短，但已形成了干扰。另外，在不同的环境温度下，壳体的硬度和外形会有所变化，其余震有时长、有时短、有时干扰大、有时干扰小，这是设计工业级产品时必须要考虑的问题。
 3、电路串扰。超声波发射时的瞬间电流很大，例如某种工业级连续测距产品瞬间电流会有15A，通常的产品也能达到1A，瞬间这么大的电流会对电源有一定影响，并干扰接收电路。通过改善电源设计可以缓解这种情况，但在低成本设计中很难根除。所以每次发射完毕，接收电路还需要一段时间稳定工作状态。在此期间，其输出的信号很难使用。
 
消除上述现象的方法之一就是在检测的时候多次循环检测，取平均值，也就是加权平均滤波，一个简单的滤波处理。就是下面这一段：
 
            u32 t = 0;
	int i = 0;
	float lengthTemp = 0;
	float sum = 0;
	while(i!=5)
	{
	TRIG_Send = 1;      //发送口高电平输出
	Delay_Us(20);
	TRIG_Send = 0;
	while(ECHO_Reci == 0);      //等待接收口高电平输出
		OpenTimerForHc();        //打开定时器
		i = i + 1;
		while(ECHO_Reci == 1);
		CloseTimerForHc();        //关闭定时器
		t = GetEchoTimer();        //获取时间,分辨率为1US
		lengthTemp = ((float)t/58.0);//cm
		sum = lengthTemp + sum ;
	
}
	lengthTemp = sum/5.0;