SCMA仿真程序（中文注释）

function [LLR_out] = SCMA_decoder(y,sigma2, config, h) %% SCMA解码器 % @params y 接收信号（K*Nsyms矩阵） % @params sigma2 接收端噪声功率 % @params config 配置数据 % @params h 瑞利衰落（K*V*Nsyms矩阵）,假设完美信道估计 % @return LLR 比特对数似然比（V*N矩阵） % @ps dr=3，即资源节点的度等于3，一般情况以后再更新 % @ps dv=2, 即用户节点的度等于2 % @ps M=4, 即调制指数等于4 LLR_out = zeros(config.V, config.N); LLR = zeros(config.V*log2(config.M),config.Nsyms); % 用于parfor的输出(输出向量 % 减小parfor输入communication overhead M = config.M; K = config.K; V = config.V; F = config.F; Niter = config.Niter; CB = config.CB; parfor n = 1:config.Nsyms % 计算接收信号的条件概率密度指数部分，P(yn|xn) f = zeros(M, M, M, K); for k = 1:K ind = find(F(k,:)==1);% 与资源块k相连的用户index for m1 = 1:M for m2 = 1:M for m3 = 1:M f(m1,m2,m3,k) = -(1/2/sigma2)*abs(y(k,n)- ... (CB(k,m1,ind(1))*h(k,ind(1),n)+ ... CB(k,m2,ind(2))*h(k,ind(2),n)+ ... CB(k,m3,ind(3))*h(k,ind(3),n)) )^2; end end end end % step1: 外信息初始化，后验概率对数向量 Ivg = log(1/M*ones(K, V, M)); Igv = zeros(K, V, M); % 仅是为了预先分配内存提高程序效率 Igv_tmp = zeros(K, V, M);% 归一化前 % step2: 迭代过程 for iter = 1:Niter % 信息传递v->g,源程序把资源节点表示为g，这里沿用之前的变量符号 for k = 1:K % 每一个资源节点单独计算 ind = find(F(k,:)==1); % 与资源节点k相连的用户节点index % 资源节点k到ind(1)的外信息传递 for m1 = 1:M % 外信息向量的每一个元素单独计算 sIgv = zeros(M, M); for m2 = 1:M for m3 = 1:M sIgv(m2,m3) = f(m1,m2,m3,k)+Ivg(k,ind(2),m2)+Ivg(k,ind(3),m3); end end Igv_tmp(k,ind(1),m1) = log_sum_exp(reshape(sIgv,1,[])); end % 归一化 for m1 = 1:M Igv(k,ind(1),m1) = Igv_tmp(k,ind(1),m1) - ... log(sum(exp(Igv_tmp(k,ind(1),:)))); end % 资源节点k到ind(2)的外信息传递 for m2 = 1:M sIgv = zeros(M, M); for m1 = 1:M for m3 = 1:M sIgv(m1,m3) = f(m1,m2,m3,k)+Ivg(k,ind(1),m1)+Ivg(k,ind(3),m3); end end Igv_tmp(k,ind(2),m2) = log_sum_exp(reshape(sIgv,1,[])); end % 归一化 for m2 = 1:M Igv(k,ind(2),m2) = Igv_tmp(k,ind(2),m2) - ... log(sum(exp(Igv_tmp(k,ind(2),:)))); end % 资源节点k到ind(3)的外信息传递 for m3 = 1:M sIgv = zeros(M, M); for m1 = 1:M for m2 = 1:M sIgv(m1,m2) = f(m1,m2,m3,k)+Ivg(k,ind(1),m1)+Ivg(k,ind(2),m2); end end Igv_tmp(k,ind(3),m3) = log_sum_exp(reshape(sIgv,1,[])); end % 归一化 for m3 = 1:M Igv(k,ind(3),m3) = Igv_tmp(k,ind(3),m3) - ... log(sum(exp(Igv_tmp(k,ind(3),:)))); end end % 信息传递g->v for v = 1:V % 每一个用户节点单独计算 ind = find(F(:,v)==1); % 与用户节点v相连的资源节点index % 用户节点v到ind(1)的外信息传递 for m = 1:M Ivg(ind(1),v,m) = Igv(ind(2),v,m); % 因为dv=2，所以不用归一化 end % 用户节点v到ind(2)的外信息传递 for m = 1:M Ivg(ind(2),v,m) = Igv(ind(1),v,m); % 不用归一化 end end end % step3: 计算比特对数似然比 Q = zeros(M, V); ... % 所有用户组成的，每一个符号取值后验概率对数构成的矩阵 LLR_tmp = zeros(V*log2(M),1);% 用于parfor的输出 for v = 1:V ind = find(F(:,v)==1); for m = 1:M Q(m,v) = Igv(ind(1),v,m) + Igv(ind(2),v,m); end Q(:,v) = Q(:,v) - log(sum(exp(Q(:,v)))); % 归一化 LLR_tmp(log2(M)*(v-1)+1) = log(exp(Q(1,v))+exp(Q(2,v))) - ... log(exp(Q(3,v))+exp(Q(4,v))); % 第一个bit LLR_tmp(log2(M)*(v-1)+2) = log(exp(Q(1,v))+exp(Q(3,v))) - ... log(exp(Q(2,v))+exp(Q(4,v))); % 第二个bit end LLR(:,n) = LLR_tmp; end % parfor for n = 1:config.Nsyms LLR_out(:,log2(config.M)*(n-1)+1) = downsample(LLR(:,n),log2(config.M),0); LLR_out(:,log2(config.M)*(n-1)+2) = downsample(LLR(:,n),log2(config.M),1); end end function y = log_sum_exp(x) %% log(exp(a)+exp(b)+...)的近似公式,简化计算 % @params x 由a,b,...组成的行向量 xm = max(x); x = x-repmat(xm, size(x,1), 1); y = xm + log(sum(exp(x))); end