ldpc.zip_LDPC_LDPC 通信系统_ldpc码信道

%************ preparation ****************************** M=100; N=200; iter = 5; onePerCol=3; sum=0; %*********** creat H matrix ***************************** %H=makeLdpc(M, N, 1, 1, onePerCol); if N/M ~= 2 fprintf('Code rate must be 1/2\n'); end onePerRow = (N/M)*onePerCol; fprintf('Creating LDPC matrix...\n'); %case {1} % Distribute 1s uniformly at random within column for i = 1:N onesInCol=zeros(M,N); onesInCol(: , i)= randperm(M)'; end % Create non zero elements (1s) index r = reshape(onesInCol(1:onePerCol, :), N*onePerCol, 1); tmp = repmat([1:N], onePerCol, 1); c = reshape(tmp, N*onePerCol, 1); % Make the number of 1s between rows as uniform as possible % Order row index [r, ix] = sort(r);%升续排列这60个值 % Order column index based on row index for i = 1:N*onePerCol cSort(i, :) = c(ix(i)); end % Create new row index with uniform weight tmp = repmat([1:M], onePerRow, 1); r = reshape(tmp, N*onePerCol, 1); % Create sparse matrix H % Remove any duplicate non zero elements index using logical AND S = and(sparse(r, cSort, 1, M, N), ones(M, N)); H = full(S); %*********** code *************************************** seldata=round(rand(1,M)); checkcode=seldata*H; for i=1:M checkcode(1,i)=mod(checkcode(1,i),2); end code=zeros(1,N); code(1,M+1:N)=checkcode(1,1:M); code(1,1:M)=seldata; %1*N rx=code; %*********** decode ************************************** %vHat=decodeProbDomain(rx,H,iter); [M N] = size(H); % Prior probabilities P1 = ones(size(rx))./(1 + exp(-2*rx)); P0 = 1 - P1; K0 = zeros(M, N); K1 = zeros(M, N); rji0 = zeros(M, N); rji1 = zeros(M, N); qij0 = H*repmat(P0', 1, N);%M应改成N qij1 = H*repmat(P1', 1, N); % ----- Horizontal step ----- for i = 1:M % Find non-zeros in the column c1 = find(H(i, :));%找到非零，即1的位置 for k = 1:length(c1) % Get column products of drji\c1(l) drji = 1; for l = 1:length(c1) if l~= k drji = drji*(qij0(i, c1(l)) - qij1(i, c1(l))); end end % for l rji0(i, c1(k)) = (1 + drji)/2; rji1(i, c1(k)) = (1 - drji)/2; end % for k end % for i % ------ Vertical step ------ for j = 1:N % Find non-zeros in the row r1 = find(H(:, j)); for k = 1:length(r1) % Get row products of prodOfrij\ri(l) prodOfrij0 = 1; prodOfrij1 = 1; for l = 1:length(r1) if l~= k prodOfrij0 = prodOfrij0*rji0(r1(l), j); prodOfrij1 = prodOfrij1*rji1(r1(l), j); end end % for l % Update constants K0(r1(k), j) = P0(j)*prodOfrij0; K1(r1(k), j) = P1(j)*prodOfrij1; % Update qij0 and qij1 qij0(r1(k), j) = K0(r1(k), j)./(K0(r1(k), j) + K1(r1(k), j)); qij1(r1(k), j) = K1(r1(k), j)./(K0(r1(k), j) + K1(r1(k), j)); end % for k % Update constants Ki0 = P0(j)*prod(rji0(r1, j)); Ki1 = P1(j)*prod(rji1(r1, j)); % Get Qj Qi0 = Ki0/(Ki0 + Ki1); Qi1 = Ki1/(Ki0 + Ki1); % Decode Qj if Qi1 > Qi0 vHat(j) = 1; else vHat(j) = 0; end end % for j for i=1:M chazhi=abs(vHat(1,i)-seldata(1,i)); sum=sum+chazhi; end