课程作业记录5：LDPC码在删除/高斯信道下的误帧率计算/Matlab仿真

由于该过程涉及到了多个过程，为增强代码的可读性和可移植性，我们将其分为了主函数和子函数，其各自功能叙述如下：
1） 主函数：整个仿真过程的主架构。
2） 子函数1，MakeLdpc：用于生成一个LDPC矩阵。
3） 子函数2，MakeParityChk：用于生成校验矩阵。
4） 子函数3，BPBEC：删除信道下的BP译码过程。

下面分块给出代码：

1.主函数

% 删除信道下的误帧率计算
clc; 
clear all;
% 给定LDPC码矩阵大小
M = 64;
N = 128;
% 给定每列1的个数
onePerCol = 3;
% 删除概率del
del = 0.2:0.05:0.4;
% 给定迭代次数
iter = 50;
% 给定帧数（64比特一帧）
frame = 500;
% 给定删除元素
% 生成LDPC矩阵
H = makeLdpc(M, N, onePerCol);
for i=1:1:length(del)
   fer(i) = 0;
   error(i) = 0;
   % 随机生成0/1作为源码
   dSource = round(rand(M, frame));
   for j = 1:frame  
      % 编码
      [c, newH] = makeParityChk(dSource(:, j), H);
      u = [c; dSource(:, j)];
      % BPSK调制
      bpskMod = 2*u - 1;
      % 删除信道
      delcheck=randperm(128,128);
      tx=((delcheck-128*del(i))>0)'.*bpskMod+((delcheck-128*del(i))<=0)'*0.1;
      % 解码
      vhat = BPBEC(tx, newH, del(i), iter);
      % 累计帧错误
      error(i)=(sum(u~=vhat')~=0)+error(i);
   end
   fer(i)=error(i)/frame;
end
semilogy(del, fer,'o-');
xlabel('DEL');
ylabel('FER');
title('FER vs DEL with BEC channel');

2.makeLdpc

function H = makeLdpc(M, N, onePerCol)
% 该子函数用于生成一个LDPC矩阵
% 参数说明：
%  M        : 行数
%  N        : 列数
%  onePerCol: 每列1的个数
%  H        : LDPC矩阵                
%每行1的个数
onePerRow = (N/M)*onePerCol;
for i = 1:N
    onesInCol(:, i) = randperm(M)';
end
% 生成非零元素索引
r = reshape(onesInCol(1:onePerCol, :), N*onePerCol, 1);
tmp = repmat([1:N], onePerCol, 1);
c = reshape(tmp, N*onePerCol, 1);   
% 行索引
[r, ix] = sort(r);
% 根据行索引生成列索引
for i = 1:N*onePerCol
    cSort(i, :) = c(ix(i));
end     
% 生成新的行索引
tmp = repmat([1:M], onePerRow, 1);
r = reshape(tmp, N*onePerCol, 1);   
% 生成稀疏矩阵H
% 消除任意两个连续的非零元素
S = and(sparse(r, cSort, 1, M, N), ones(M, N));
H = full(S);      
% 检查矩阵是否存在只有一个1或没有1的行
for i = 1:M
    n = randperm(N);
    % 如果该行没有1，加上两个1
    if length(find(r == i)) == 0
        H(i, n(1)) = 1;
        H(i, n(2)) = 1;
    % 如果该行只有一个1，加上一个1  
    elseif length(find(r == i)) == 1
        H(i, n(1)) = 1;
    end
end % for i
%消除四环
for i = 1:M
    % 寻找行列对
    for j = (i + 1):M         
        w = and(H(i, :), H(j, :));
        c1 = find(w);
        lc = length(c1);
        if lc > 1
            % 如果找到四环，用0跳过连续的1
            if length(find(H(i, :))) < length(find(H(j, :)))
                % 重复该过程 
                for cc = 1:lc - 1
                    H(j, c1(cc)) = 0;
                end
            else
               for cc = 1:lc - 1
                  H(i, c1(cc)) = 0;
               end
            end % if            
         end % if
      end % for j
      end % for i

3.makeParityChk

function [c, newH] = makeParityChk(dSource, H)
% 该子函数为生成校验矩阵H函数
% 参数说明：
%  dSource : 信息源（0/1）
%  H       : LDPC矩阵        
%  c       : 校验比特
%
% 矩阵维度
[M, N] = size(H);
% 生成一个新的矩阵F，便于L、U分解
F = H;
% L、U矩阵
L = zeros(M, N - M);
U = zeros(M, N - M);
% 重新排序M x (N - M)复矩阵
for i = 1:M
    % 对角线非零元素
    [r, c] = find(F(:, i:end));
    colWeight = sum(F(:, i:end), 1) - 1;
    rowWeight = sum(F(i, :), 2) - 1;
    % 非零元素行索引
    rowIndex = find(r == i);
    % 最小积
    [x, ix] = min(colWeight(c(rowIndex))*rowWeight);
    % 按照索引重新排序矩阵F
    chosenCol = c(rowIndex(ix)) + (i - 1);
    % 矩阵H、F列重新排序
    tmp1 = F(:, i);
    tmp2 = H(:, i);
    F(:, i) = F(:, chosenCol);
    H(:, i) = H(:, chosenCol);
    F(:, chosenCol) = tmp1;
    H(:, chosenCol) = tmp2;
    % 按列填充L、U矩阵
    L(i:end, i) = F(i:end, i);
    U(1:i, i) = F(1:i, i);
    if i < M           
        % 在第i列中，寻找下一个非零元素所在位置
        [r2, c2] = find(F((i + 1):end, i));          
        % 找到后，插入一行
        F((i + r2), :) = mod(F((i + r2), :) + repmat(F(i, :), length(r2), 1), 2);
    end % if 
end % for i
% 源码Z 
z = mod(H(:, (N - M) + 1:end)*dSource, 2);
% 校验比特
c = mod(U(Lz), 2); 
% 返回值为校验矩阵 
newH = H;

4.BPBEC

function vHat = BPBEC(rx, H, del, iteration)
% 该子函数为对数域和积算法（BP译码）
% 参数说明：
%  rx        : 收到信号列向量
%  H         : LDPC矩阵
%  del       : 删除概率
%  iteration : 迭代次数
%  vHat      : 解码向量（0/1）
[M N] = size(H);
% 先验对数概率
Lci = (-4*rx./del)';
% 初始化
Lrji = zeros(M, N);
Pibetaij = zeros(M, N);
% 将L（ci）矩阵与H矩阵相关
Lqij = H.*repmat(Lci, M, 1);
% 得到非零元素
[r, c] = find(H);
% 迭代
for n = 1:iteration
   % L（qij）的符号与角度  
   alphaij = sign(Lqij);   
   betaij = abs(Lqij);
   for l = 1:length(r)
      Pibetaij(r(l), c(l)) = log((exp(betaij(r(l), c(l))) + 1)/(exp(betaij(r(l), c(l))) - 1));
   end
   for i = 1:M
      % 列非零元素
      c1 = find(H(i, :));
      % Pi(betaij)的总和       
      for k = 1:length(c1)
         sumOfPibetaij = 0;
         prodOfalphaij = 1;
         % Pi(betaij)c1(k)的总和
         sumOfPibetaij = sum(Pibetaij(i, c1)) - Pibetaij(i, c1(k));
         % 为避免0做除数，我们计算Pi(sum(Pi(betaij)))
         if sumOfPibetaij < 1e-20
            sumOfPibetaij = 1e-10;
         end         
         PiSumOfPibetaij = log((exp(sumOfPibetaij) + 1)/(exp(sumOfPibetaij) - 1));
         % 乘上alphaijc1(k) 
         prodOfalphaij = prod(alphaij(i, c1))*alphaij(i, c1(k));
         % 更新L(rji)
         Lrji(i, c1(k)) = prodOfalphaij*PiSumOfPibetaij;
      end % for k
   end % for i
   for j = 1:N
      % 行非零元素
      r1 = find(H(:, j));
      for k = 1:length(r1)        
         % 用L(rij)r1(k)的和更新L(qij)
         Lqij(r1(k), j) = Lci(j) + sum(Lrji(r1, j)) - Lrji(r1(k), j);
      end % for k
      LQi = Lci(j) + sum(Lrji(r1, j));
      % 判定L(Qi)
      if LQi < 0
         vHat(j) = 1;
      else
         vHat(j) = 0;
      end                
   end % for j
end % for n

对于高斯信道而言，该仿真过程代码仅有少数部分不同，我们将改动的地方列举如下：

1. 主函数信道部分，由删除信道变为了高斯信道。

% 高斯白噪声
N0 = 1/(exp(EbN0(i)*log(10)/10));
tx = bpskMod + sqrt(N0/2)*randn(size(bpskMod));

2. BP译码子函数部分，输入的变量与先验概率有所变化。

function vHat = decodeLogDomain(rx, H, N0, iteration)
% 先验对数概率
Lci = (-4*rx./N0)';

3. 主函数的输入量变为了Eb/N0。

% EbN0 in db
EbN0 = [0 0.5 1 1.5];

最后

以上就是爱笑小松鼠为你收集整理的课程作业记录5：LDPC码在删除/高斯信道下的误帧率计算/Matlab仿真的全部内容，希望文章能够帮你解决课程作业记录5：LDPC码在删除/高斯信道下的误帧率计算/Matlab仿真所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。