www.gusucode.com > matlab编程LDPC编码的BP译码算法的程序源码 > matlab编程LDPC编码的BP译码算法的程序源码/code/main.m
% % %---------------------------------------------% function [x_hat, success, k] = ldpc_decode(f0,f1,H) % decoding of binary LDPC as in Elec. Letters by MacKay&Neal 13March1997 % For notations see the same reference. % function [x_hat, success, k] = ldpc_decode(y,f0,f1,H) % outputs the estimate x_hat of the ENCODED sequence for % the received vector y with channel likelihoods of '0' and '1's % in f0 and f1 and parity check matrix H. Success==1 signals % successful decoding. Maximum number of iterations(重复迭代) is set to 100. % k returns number of iterations until convergence(收敛). % % Example: % We assume G is systematic G=[A|I] and, obviously, mod(G*H',2)=0 % sigma = 1; % AWGN noise deviation % x = (sign(randn(1,size(G,1)))+1)/2; % random bits % y = mod(x*G,2); % coding % z = 2*y-1; % BPSK modulation % z=z + sigma*randn(1,size(G,2)); % AWGN transmission % % f1=1./(1+exp(-2*z/sigma^2)); % likelihoods % by兰朋朋 ./矩阵元素右除 A./B表示A(i,j)/B(i,j) % f0=1-f1; % [z_hat, success, k] = ldpc_decode(z,f0,f1,H); % x_hat = z_hat(size(G,2)+1-size(G,1):size(G,2)); % x_hat = x_hat'; % Copyright (c) 1999 by Igor Kozintsev igor@ifp.uiuc.edu % $Revision: 1.1 $ $Date: 1999/07/11 $ % fixed high-SNR decoding [m,n] = size(H); if m>n, H=H'; [m,n] = size(H); end if ~issparse(H) % make H sparse if it is not sparse yet,%by兰朋朋 函数的功能是如果该函数不是稀疏矩阵,就把它转化为稀疏矩阵 %by兰朋朋 在MATLAB中,issparse函数用于判断一个稀疏矩阵存储方式是否是sparse storage organization %by兰朋朋 如果稀疏矩阵的存储方式是sparse storage organization,则返回逻辑1;否则返回逻辑0。 [ii,jj,sH] = find(H); %by兰朋朋 返回矩阵H中的非零元素所在的行ii和所在的列jj,和矩阵的非零元素sH H = sparse(ii,jj,sH,m,n); %by兰朋朋 [i,j,s] = find(S); %by兰朋朋 [m,n] = size(S); %by兰朋朋 S = sparse(i,j,s,m,n); end %by兰朋朋 BP译码算法 %initialization(没看懂) [ii,jj] = find(H); % subscript index to nonzero elements of H indx = sub2ind(size(H),ii,jj); % linear index to nonzero elements of H %by兰朋朋 该函数用于把数组中元素下标转换为该元素在数组中对应的索引值 q0 = H * spdiags(f0(:),0,n,n); %by兰朋朋 A=spdiags(B,d,m,n),返回带状稀疏矩阵A。 %by兰朋朋 参数m,n为原带状矩阵的行数与列数。 %by兰朋朋 B为r×p阶矩阵,这里r=min(m,n),p为原带状矩阵所有非零对角线的条数,矩阵B的第i列即为原带状矩阵的第i条非零对角线, d为长度为p的向量。 sq0 = full(q0(indx)); %by兰朋朋 在MATLAB中,该函数用于把一个稀疏矩阵(sparse matrix)转换成一个全矩阵(full matrix)。 sff0 = sq0; q1 = H * spdiags(f1(:),0,n,n); sq1 = full(q1(indx)); sff1 = sq1; %迭代过程 %iterations k=0; success = 0; max_iter = 100; while ((success == 0) & (k < max_iter)), k = k+1; %horizontal step,第一步 sdq = sq0 - sq1; sdq(find(sdq==0)) = 1e-20; % if f0 = f1 = .5 %by兰朋朋 若先验概率一样的话,令sdq为一个很小很小的数 dq = sparse(ii,jj,sdq,m,n); Pdq_v = full(real(exp(sum(spfun('log',dq),2)))); % this is ugly but works :) %将dq在对数域求和,即将dq相乘,构成数组 Pdq_v Pdq = spdiags(Pdq_v(:),0,m,m) * H; sPdq = full(Pdq(indx)); sr0 = (1+sPdq./sdq)./2; sr0(find(abs(sr0) < 1e-20)) = 1e-20; sr1 = (1-sPdq./sdq)./2; sr1(find(abs(sr1) < 1e-20)) = 1e-20; r0 = sparse(ii,jj,sr0,m,n); r1 = sparse(ii,jj,sr1,m,n); %vertical step Pr0_v = full(real(exp(sum(spfun('log',r0),1)))); Pr0 = H * spdiags(Pr0_v(:),0,n,n); sPr0 = full(Pr0(indx)); Q0 = full(sum(sparse(ii,jj,sPr0.*sff0,m,n),1))'; sq0 = sPr0.*sff0./sr0; Pr1_v = full(real(exp(sum(spfun('log',r1),1)))); Pr1 = H * spdiags(Pr1_v(:),0,n,n); sPr1 = full(Pr1(indx)); Q1 = full(sum(sparse(ii,jj,sPr1.*sff1,m,n),1))'; sq1 = sPr1.*sff1./sr1; sqq = sq0+sq1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%自己加的 % if sqq<=0 % sqq=1e-20; % end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% sq0 = (sq0+10*eps)./(sqq+10*eps); sq1 = (sq1+10*eps)./(sqq+10*eps); %tentative decoding QQ = Q0+Q1; Q0 = Q0./QQ; Q1 = Q1./QQ; x_hat = (sign(Q1-Q0)+1)/2; if rem(H*x_hat,2) == 0, success = 1; end end