前一篇写直接序列扩频系统仿真的文章中的模型现在发现了严重错误，并且那个模型是针对我们老师布置的课程设计而设计的，对于各位读者的参考性可能不强。于是决定重做这一个模型，同时加入了对任意用户同时通信的支持，实现只要修改参数，就可以对任意名用户同时进行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:原理展示代码，用于原理分析，频谱分析等。链接：https://pan.baidu.com/s/1viNzbBHz3qKz9FJ6obITIQ 提取码：od8o
载波的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_num plot(foo_x,user_errorRate(:,i)); legend_plot{i} = ['用户',num2str(i)]; hold on; end title('用户误码率曲线'); 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); endwalsh.m
function walsh = walsh(N) % 产生Walsh函数通用函数 % N:Walsh函数阶数,当N不是2的幂时,通过向无穷大取整使得所得Walsh阶数为2的幂 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); endspreadSpectrum.m
function res = spreadSpectrum(userCode,spreadCode) % 实现对每个用户进行扩频 % userCode:用户码元，每一行为一个用户所有的码元 % spreadCode:扩频码，每一行作为一个用户自己的扩频码 [row_userCode,col_userCode] = size(userCode); [row_spreadCode,col_spreadCode] = size(spreadCode); if row_spreadCode < row_userCode error('扩频码行数不能小于用户数'); end res = int8(zeros(row_userCode,col_userCode * col_spreadCode)); %对每一个用户进行扩频 for i = 1:row_userCode for j = 1:col_userCode res(i,(j-1)*col_spreadCode+1:j*col_spreadCode) = userCode(i,j).*spreadCode(i,:); end end enddeSpreadSpectrum.m
function res = deSpreadSpectrum(userSpreadCode,spreadCode,user_num) %本函数用于解扩 %userSpreadCode:需要解扩的用户码元 %spreadCode:用于扩频的随机码 %user_num:用户数量 [~,col] = size(userSpreadCode); res = zeros(user_num,col); %为每一个用户分别解扩 for i = 1:user_num res(i,:) = bitMultiple(userSpreadCode,spreadCode(i,:)); end endcodeJudge.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) < 0 res(i,j) = -1; end end end endoverlay.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_postive res(i) = 1; elseif foo(i) < threshhold_negative res(i) = -1; end end endcarrierGen.m
function [res] = carrierGen(freq,sampleRate,size) %carrierGen 产出载波 % 参数：freq:频率 sampleRate:采样率 size:载波长度 n=0:size-1; t=n/sampleRate;%时间序列 res=sin(2*pi*freq*t); endmodulate.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(:); enddemodulate.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_postive res(i) = 1; elseif res_arrayGroupSum(i) < threshhold_negative res(i) = -1; end end endcompare.m
function accuracy = compare(input1,input2) %比较两个数组 %input1:矩阵1，每一行为一个用户所有码元 %input2:矩阵2 %accuracy:正确率数组 [row,~] = size(input1); accuracy = zeros(row,1); for i = 1:row temp = (input1(i,:) == input2(i,:)); res = length(find(temp == 1)); sizeInput = max(length(input1),length(input2)); accuracy(i) = res/double(sizeInput); end end