8位16位64位等任意数量用户CDMA直接序列扩频通信系统的Matlab仿真
前言模型缺点 原理码元扩频扩频码 叠加二项分布中心极限定理叠加的实现 调制解调解扩&码元判决 实验结果仿真代码代码说明代码下载链接代码可修改的参数主函数双极性码生成模块walsh矩阵生成模块扩频模块解扩模块码元判决模块叠加模块载波生成模块调制模块解调模块误码率计算模块前言
前一篇写直接序列扩频系统仿真的文章中的模型现在发现了严重错误，并且那个模型是针对我们老师布置的课程设计而设计的，对于各位读者的参考性可能不强。于是决定重做这一个模型，同时加入了对任意用户同时通信的支持，实现只要修改参数，就可以对任意名用户同时进行CDMA通信的仿真。
模型缺点
模型所生成的误码率曲线是符合理论预期的，但模型有一个缺点，那就是无论信噪比有多高，误码率始终无法为0，而是趋于某个大于0的常数，该常数会随着用户数量的增加而增加。这是因为n个用户的扩频码叠加会产生n+1种结果，但程序使用了BPSK调制，该调制仅允许-1,0,1三种码元出现，使得扩频码叠加时只能把上述n+1种结果“压缩”成3种结果，导致了系统无法以0误码率传输数据，这也恰恰说明了在用户码元中加入纠错码的重要性。
原理
码元
用户码元采用双极性码，基于随机数法生成。
扩频
扩频原理如下图所示。
假如码元为1，则发送扩频码，假如码元是-1（如果是单极性码，那么是0），则把扩频码取反后发送，把要发送的扩频码连接在一起即为扩频结果。扩频前后频谱对比如下图所示：
可以看到，扩频后的信号在调制后相比不扩频的信号占用更多的频谱资源，其功率谱密度分散在了增个频谱上，这样有利于实现隐蔽通信，也有利于抵抗高斯干扰。
扩频码
扩频码采用了Walsh矩阵法生成。Walsh矩阵的任意两行间都是绝对正交，模型中取Walsh矩阵的一行作为一个用户的扩频码，每个用户取不同的行，从而使得每个用户的扩频码都是正交的。
叠加
扩频后n个用户的扩频码仍是双极性码1和-1，这些扩频码对应相加后，会有n+1种结果。比如有8个用户，那么相加后的可能结果有-8,-6,…,0,2,…,8。因为接下来的调制使用BPSK调制，仅允许-1,1,0三种码元，所以要对相加后的结果进行某种“压缩”处理，把它们“映射”到-1,1,0上。使用什么办法进行“压缩”呢？
二项分布
分析发现，相加后的结果的概率分布满足二项分布。假设在n个用户的结果中，恰有k个用户是1的概率为P{X=k}，因为用户码只有1或-1且它们是等概率出现的，所以：
P { X = k } = C n k ( 1 / 2 ) n P\{X=k\}=C_{n}^{k}(1/2)^{n} P{X=k}=Cnk(1/2)n
当有k个用户是1时，相加后的结果y=2k-n。不难证明，二项分布的数学期望是np，方差是np(1-p)。
可以基于以上二项分布，计算出把n+1种可能的结果进行区分的阈值。但是这样计算出来的阈值同时也是可能的结果，那么在分区是是否包含该阈值呢？比如在n=8的情况下，假如算出来的阈值是-2，那么小于-2的结果当然地区分为-1，那等于-2呢？这种情况下无论是纳入-1还是纳入0，都会对结果造成重大印象，所以有必要按照更好的方法计算分区阈值。
中心极限定理
根据林德伯格——列维(Lindeberg-Levy)中心极限定理，当随机变量足够多时，均值为μ，方差为σ^2>0的独立同分布的随机变量X1,X2,…,Xn的算术平均
X ‾ = 1 n ∑ i = 1 n X i \overline{X}=\cfrac{1}{n}\sum_{i=1}^n X_i X=n1i=1∑nXi
近似地服从均值为μ，方差为σ^2/n的正态分布。
也就是说，上述二项分布在随机变量足够多时，随机变量的算术平均会趋向于服从均值为np，方差为p(1-p)的正态分布。
叠加的实现
叠加的结果有3种，假设这3种各占1/3的概率，那么就可以根据正态分布的概率模型，结合实际数据的算术平均值与标准差，反算出使得概率为1/3的上限值，从而根据该值设置阈值。
实际生成的数据样本是有限的，为了n个用户的n+1种相加结果可以平等地“缩放”到-1,0,1这3种结果上，首先计算上述相加结果的方差与标准差，然后作为参数调用normcdf函数（该函数是对正态分布密度函数的变上限积分），通过解方程normcdf(x,算术平均值,标准差)=1/3解出x作为负阈值；解方程normcdf(x,算术平均值,标准差)=2/3解出x作为正阈值，从而实现阈值随着实际数据分布的变化而变化，公平地把n+1种结果变换到3种结果中。
不过，也正是因为这样变换，使得叠加后发生了“精度丢失”，使得即使信噪比很高，也仍有误码存在。
调制
调制使用BPSK调制，接收调制的有3种码元：-1，0，1。与扩频类似，如果码元是1，则把载波发送出去；如果码元是-1，则把载波取反后发送出去，如果是0，则发送一串载波长度的0码。调制示例图如下图所示：
解调
这是调制的逆过程，采用相干解调法，把接收到的结果与载波循环点对点相乘从而实现鉴相。因为
s i n ( x ) ∗ s i n ( x ) ≥ 0 在 x ⊂ R 恒成立 sin(x)*sin(x) \ge0在x\sub R恒成立 sin(x)∗sin(x)≥0在x⊂R恒成立
所以通过相乘后，分组判决即可判决出调制前是-1还是0还是1。判决采用与叠加类似的方法，通过正态分布方程解算出阈值，实现判决。
解扩&码元判决
解扩采用扩频后的码元与用户相应扩频码点对点相乘的方法，把相乘结果分组相加，根据结果大于0还是小于0，判决出用户的原始码元。
实验结果
实验的用户码元是一千个，64个用户，扩频增益是64，每个用户的误码率曲线如下图所示
平均误码率如下图所示
以上结果与理论是符合的。
仿真代码
代码说明
代码有两条主线：
main.m:这个是实验的代码，会绘制实验结果图；thorydisp.m:原理展示代码，用于原理分析，频谱分析等。
代码下载链接
尽量使用GitHub版本，度娘偶尔会吞文件
GitHub: /highskyno1/CDMA-simulation
度盘：/s/18_SQl7P7Ml9iwtMsWLh4aA?pwd=19q6
提取码：19q6
代码可修改的参数
用户数每个用户的码元数walsh码阶数，也是扩频增益，必须为2的幂次方且必须等于或大于用户数载波频率载波采样率载波采样点数载波的3个参数在修改时，务必注意要能采到一个完整的正弦波波形。
主函数
main.m
%{本程序用于多用户CDMA通信系统仿真，采用Walsh码作为扩频码传输过程使用BPSK调制。%}close all;user_num = 64; %用户数，必须少于或等于walsh码阶数user_code = 1e4; %每个用户的码元数walsh_order = 64; %walsh码阶数,也是扩频增益,必须是2的幂次方carrier_freq = 1e3; %载波频率SampleRate = 25*1e3; %载波采样率SamplePoint = 25; %载波采样点数%产生用户码元userCode = genBipolar(user_num,user_code);%产生Walsh矩阵walshCode = walsh(walsh_order);%扩频res_spreadSpectrum = spreadSpectrum(userCode,walshCode);%扩频后的码元相加res_overlay = overlay(res_spreadSpectrum);%生成载波carrier = carrierGen(carrier_freq,SampleRate,SamplePoint);%调制send = modulate(res_overlay,carrier);%尝试不同的信噪比(dB)snr_max = 20;snr_min = -20;snr_div = 1;snr_array = floor((snr_max-snr_min)/snr_div);%记录每个用户每个信噪比下的用户误码率user_errorRate = zeros(snr_array,user_num);%记录平均误码率mean_errorRate = zeros(1,snr_array);parfor snr_index = 1:snr_array%当前信噪比snr = snr_min + snr_div * snr_index;%加噪receive = awgn(send,snr,powerCnt(carrier));%解调(相干解调)res_demodulate = demodulate(receive,carrier);%解扩res_deSpreadSpectrum = deSpreadSpectrum(res_demodulate,walshCode,user_num);%码元判决res_codeJudge = codeJudge(res_deSpreadSpectrum,walsh_order);%计算误码率error_rate = 1 - compare(res_codeJudge,userCode);%保存误码率user_errorRate(snr_index,:) = error_rate;mean_errorRate(snr_index) = mean(error_rate);fprintf('当前信噪比%f，平均误码率%f\n',snr,mean(error_rate));end%绘制误码率曲线foo_x = linspace(snr_min,snr_max,snr_array);%生成标签legend_plot = cell(1,user_num);figure;for i = 1:user_numplot(foo_x,user_errorRate(:,i));legend_plot{i} = ['用户',num2str(i)];hold on;endtitle('用户误码率曲线');legend(legend_plot);xlabel('信噪比(dB)');ylabel('误码率');%绘制平均误码率曲线figure;%计算平均误码率plot(foo_x,mean_errorRate);title('平均误码率曲线');xlabel('信噪比(dB)');ylabel('误码率');
双极性码生成模块
genBipolar.m
function res = genBipolar(num_user,num_code)% 产生双极性码% num_code:双极性码的规模% num_user:用户数% res:[用户1双极性码;用户二双极性码。。。;]res = rand(num_user,num_code);res = value2Bipolar(res);end
walsh矩阵生成模块
walsh.m
function walsh = walsh(N)% 产生Walsh函数通用函数% N:Walsh函数阶数,当N不是2的幂时,通过向无穷大取整使得所得Walsh阶数为2的幂M = ceil(log2(N));wn = 0;for i = 1:Mw2n = [wn,wn;wn,~wn];wn = w2n;endwn(wn == 0) = -1;walsh = int8(wn);end
扩频模块
spreadSpectrum.m
function res = spreadSpectrum(userCode,spreadCode)% 实现对每个用户进行扩频% userCode:用户码元，每一行为一个用户所有的码元% spreadCode:扩频码，每一行作为一个用户自己的扩频码[row_userCode,col_userCode] = size(userCode);[row_spreadCode,col_spreadCode] = size(spreadCode);if row_spreadCode < row_userCodeerror('扩频码行数不能小于用户数');endres = int8(zeros(row_userCode,col_userCode * col_spreadCode));%对每一个用户进行扩频for i = 1:row_userCodefor j = 1:col_userCoderes(i,(j-1)*col_spreadCode+1:j*col_spreadCode) = userCode(i,j).*spreadCode(i,:);endendend
解扩模块
deSpreadSpectrum.m
function res = deSpreadSpectrum(userSpreadCode,spreadCode,user_num)%本函数用于解扩%userSpreadCode:需要解扩的用户码元%spreadCode:用于扩频的随机码%user_num:用户数量[~,col] = size(userSpreadCode);res = zeros(user_num,col);%为每一个用户分别解扩for i = 1:user_numres(i,:) = bitMultiple(userSpreadCode,spreadCode(i,:));endend
码元判决模块
codeJudge.m
function res = codeJudge(input,spreadSpectrumGain)%codeJudge 码元判决，基于分组相加判决原理%input 解扩后各个用户的矩阵，每一行为一个用户%spreadSpectrumGain 扩频增益[row,col] = size(input);res = int8(ones(row,floor(col/spreadSpectrumGain)));%对于每一个用户for i = 1:row%分组判决res_group = arrayGroupSum(input(i,:),spreadSpectrumGain);for j = 1:length(res_group)if res_group(j) < 0res(i,j) = -1;endendendend
叠加模块
overlay.m
function res = overlay(input)% 码元叠加函数，对输入的每一列进行相加，然后判决相加结果为1、-1或0% input:扩频后的用户码元%对每一列进行相加foo = sum(input(1:end,:));%计算标准差res_std = std(foo);%计算平均数res_mean = mean(foo);%计算概率为1/3时的阈值，也就是判决为-1的阈值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 = int8(zeros(1,length(foo)));for i = 1:length(foo)if foo(i) > threshhold_postiveres(i) = 1;elseif foo(i) < threshhold_negativeres(i) = -1;endendend
载波生成模块
carrierGen.m
function [res] = carrierGen(freq,sampleRate,size)%carrierGen 产出载波% 参数：freq:频率 sampleRate:采样率 size:载波长度n=0:size-1;t=n/sampleRate;%时间序列res=sin(2*pi*freq*t);end
调制模块
modulate.m
function res = modulate(source,carrier)%调制函数%source:被调制信号%carrier:载波信号sizeSource = length(source);sizeCarrier = length(carrier);res = zeros(sizeSource,sizeCarrier);for i = 1:sizeSourceres(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));%计算标准差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_postiveres(i) = 1;elseif res_arrayGroupSum(i) < threshhold_negativeres(i) = -1;endendend
误码率计算模块
compare.m
function accuracy = compare(input1,input2)%比较两个数组%input1:矩阵1，每一行为一个用户所有码元%input2:矩阵2%accuracy:正确率数组[row,~] = size(input1);accuracy = zeros(row,1);for i = 1:rowtemp = (input1(i,:) == input2(i,:));res = length(find(temp == 1));sizeInput = max(length(input1),length(input2));accuracy(i) = res/double(sizeInput);endend
