前言

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

1. 如果想测试多个用户的情况，或者了解原理图如何绘制，欢迎参考我的另一篇文章《8位16位64位等任意数量用户CDMA直接序列扩频通信系统的Matlab仿真》，在其中可以设置用户数量，代码相比于本文也更加简介。
2. 想了解更多关于PSK的非相干调制解调，可以参考我的另外两篇文章基于延迟相乘鉴相的2QPSK调相通信系统仿真和基于延迟相乘鉴相的4Q-DPSK差分通信系统仿真。

文章迭代更新

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

仿真流程图

关键技术细节

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阶时，不同振幅下信噪比误码率曲线

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

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

可见，扩频增益更大的用户2在同信噪比下误码率比用户1更小。

观察以上曲线可以发现，随着信噪比的上升，误码率以近似于抛物线的形式在不断下降，并最终等于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. 最小最大载波振幅，以及其尝试步进。

！下载链接！

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

GitHub：https://github.com/highskyno1/2users-CDMA-simulation
链接：https://pan.baidu.com/s/1I_wOD-OBlZ-li0Ge4wC9Vw?pwd=yq88
提取码：yq88

核心代码

主函数

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(以上)
    i7-12700K:实测速率 0.91秒/轮回 @1280*2&1280码元20&10扩频增益 12线程
    修复解调模块阈值解算的Bug，改为利用"norminv"解算阈值，
    可见运行速度显著提高，同时改为var计算方差，感谢网友的提醒。
2022/12/17(以上)
%}
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 = 40;    %尝试的最大信噪比
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概率一样(即先验概率相等)，
    则合起来为1或-1的概率为1/4，0的概率为1/2，
    为了处理方便，根据大数定理，近似地视为正态分布，
    从而利用正态分布概率密度计算阈值。
%}
%计算方差
res_var = var(input);
%计算平均数
res_mean = mean(input);
%计算概率为1/4时的阈值，小于该阈值判决为-1的阈值
threshhold_negative = norminv(1/4,res_mean,res_var);
%计算概率为3/4时的阈值，大于该阈值判决为1的阈值
threshhold_postive = norminv(3/4,res_mean,res_var);
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