大家好，又见面了，我是你们的朋友全栈君。

前言

直接扩频序列调制是用速率很高的伪噪声码序列与信息码序列模二相加（波形相乘）后得到复合码序列，用复合码序列去控制载波相位，从而获得直接扩频序列信号的。直接扩频通信具有低截获概率、抗干扰能力强以及易于实现码分多址等优点，在抗干扰通信及民用移动通信中都得到了广泛的应用。

文章迭代更新

文章进行了大改
发现了为什么振幅是0.2和0.4时，系统无法正常工作的原因。因为bitMultiple函数把运算结果格式转换成了int8，导致精度大量失真！当时转换格式是为内存空间与运行速度做打算，结果今天发现做了负优化！修改后，发现载波振幅对于系统的误码率曲线几乎无影响。
另外，解调函数也做了修改，会根据输入自动计算判决阈值，解决了之前人为设定阈值
的局限，阈值采用正态分布解算法，详见demodulate函数内注释。
新版修改了以下函数：
1. bitMultiple
2. demodulate
3. 新增arrayGroupSum函数
2022/1/25：
修改了以下函数，解决M序列生成bug
1. MseqGen
2. Oct2Bin

仿真流程图

关键技术细节

1. 仿真中，用户码元使用双极性(1和-1)码。
2. 代码的最后通过还原信号与码元做对比，计算误码率，来评价通信质量的好坏。
3. 为了研究信噪比对误码率的影响，代码中用了多线程，针对不同的信噪比进行解调去扰解扩,计算不同信噪比下的误码率。
4. 为了研究载波振幅对信噪比——误码率曲线的影响，增加了振幅的迭代。

扩频解扰解扩部分

1. 采用64阶的walsh码矩阵，由walsh函数产生。
2. 两个用户使用同一扩频矩阵的不同相位(实际上是不同一行)进行扩频。
3. 扩频时，首先对用户码元按照扩频增益，进行周期延拓，然后点乘扩频码，实现扩频。
4. 解扩时，首先对输入的去扰后的码按扩频时一样的相位对应点乘扩频码，然后进行判决。

加扰去扰部分

1. 加扰采用的是5阶的M序列，反馈系数为67(八进制)。
2. 加扰同样采用点乘，输入码元1:1地点乘加扰码。
3. 因为加扰码具有很高的自相关性，所以去扰时只需执行加扰时一样的操作即可。

加扰原理图

调制解调部分

1. 调制采用BPSK(二进制移相键控)调制。
2. 调制所使用的载波为正弦波，通过对其等间隔采样形成离散样值。
3. 调制时，一个加扰后的码元对应相乘一串正弦波的周期采样值，把相乘结果汇总后就是调制后的结果
4. 解调时，把一个码元所对应的一串正弦波的周期采样值与载波采样值进行点乘，根据点乘结果中，正数、负数还是零多，判决这是什么码元。

调制原理图

解调原理图

高斯信道部分

1. 使用matlab自带的函数awgn(input,snr(dB),inputPower)。
2. 只要传入需要加高斯噪声的信号，信噪比还有输入信号的功率，即可返回添加了高斯白噪声的信号。

实验结果

扩频增益为10时，walsh矩阵为64阶时，不同振幅下信噪比误码率曲线

扩频增益为20时，walsh矩阵为64阶时，不同振幅下信噪比误码率曲线

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

误码率随信噪比的变化曲线

观察以上曲线可以发现，随着信噪比的上升，误码率以近似于抛物线的形式在不断下降，并最终等于0，这足以证明本次直接序列扩频系统仿真符合预期结果。

源代码

可以修改的参数

1. 轮回次数。通过修改轮回次数，获得更好的曲线，或者更快的代码运行速度，代码运行参考速度：i7-4790:实测速率9.8秒/轮回 @1280*2&1280码元20&10扩频增益。
2. 用户1的码元数量与扩频增益、用户2的码元数量与扩频增益。修改时注意两用户的码元必须是walsh码阶数的增数倍，同时要保证：用户1的码元数量 * 用户1的扩频增益=用户2的码元数量 * 用户2的扩频增益。
3. walsh码的阶数。但walsh码阶数必须大于扩频增益。
4. 用户1扩频相位、用户2扩频相位。注意相位必须在1到64之间取值，两个用户的取值不能一样。
5. M序列阶数与返回系数。反馈系数取值有相关要求，详见百度，作为参数传入时要以八进制表示。
6. 半个载波周期的采样点数。
7. 最小最大信噪比，以及其尝试步进。
8. 最小最大载波振幅，以及其尝试步进。

！下载链接！

无法下载请私信或评论区评论，看到会补档

链接：https://pan.baidu.com/s/1pz3XopWOLCzqGhcC3rBrSg?pwd=lwvw
提取码：lwvw

核心代码

主函数

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扩频增益
第十三次迭代。
解决加扰函数只使用了M序列第一行的问题。
2022/1/21(以上)
i7-12700K:实测速率 2.4秒/轮回 @1280*2&1280码元20&10扩频增益
Module函数名改名为myModule，防止与Matlab内建函数重名的bug
修改testModulate函数，增加row码元长度，防止码元与扩频码被交换的bug
2022/4/28(以上)
%}
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序列
Mseq1 = Mseq1(:);
Mseq2 = MseqGen(mOrder,feedBack); %用户2加扰用的m序列
Mseq2 = Mseq2(:);
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 = myModulate(user1scarm,carrier);
user2modulate = myModulate(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');

扩频模块

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

调制模块

myModulate.m

function res = myModulate(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(1e4);
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
for i = 1:bitNum*3
if res(i) ~= 0
break;
end
end
res = flip(res(i: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