✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
1. LDPC码背景及概要
LDPC是Low Density Parity Check Code英文缩写，意为低密度奇偶校验码，最早在20世纪60年代由Gallager在他的博士论文中提出，但限于当时的技术条件，缺乏可行的译码算法，此后的35年间基本上被人们忽略。直到1993年Berrou等人发现了Turbo码，在此基础上，1995年前后MacKay和Neal等人对LDPC码重新进行了研究，提出了可行的译码算法，从而进一步发现了LDPC码所具有的良好性能，迅速引起强烈反响和极大关注。经过十几年来的研究和发展，研究人员在各方面都取得了突破性的进展，LDPC码的相关技术也日趋成熟，逐渐有了商业化的应用成果，如今LDPC码已经作为众多新一代通信标准中的信道编码方案：DVB-S2 (Digital Video Broadcasting)、IEEE 802.3an (Ethernet)、IEEE 802.16e (WiMax)、IEEE 802.11n (WiFi)、3GPP 5G标准。
LDPC码是一种稀疏校验矩阵线性分组码，在LDPC编码中，会用到一个叫做H矩阵的校验矩阵(Parity Check Matrix)，该校验矩阵为稀疏矩阵。
LDPC编码分为正则编码和非正则编码。正则编码中，校验矩阵的每行和每列中1的个数是固定的。非正则编码中，校验矩阵的每行和每列中1的个数不固定。
LDPC解码过程中，主要包括了两方面内容：硬解码(Hard Decode)和软解码(Soft Decode). LDPC解码的方法就是收到码字之后，与校验矩阵H相乘，如果是0矩阵，则说明收到的是正确码字。反之，则不正确码字，再根据相乘结果进行进一步纠错解码
2. 仿真要求及内容
本次仿真要求通过MATLAB软件对LDPC码进行编译码仿真，通过计算模拟求解归一化最小和算法α参数的最佳值和偏置最小和算法β参数的最佳值，并仿真出四种不同的译码算法下的误码率和误帧率曲线。实验的要求如下：
- 采用系统码设计，信息序列长1008比特，码长2016比特，码率1/2，即
- 给定的编码矩阵 H_block 是一个18x36的矩阵（Matrix(2016,1008)Block56.mat文件中也是只有18x36个数据），矩阵中每个元素 H(i,j) 是大小为 z*z 的循环移位矩阵（给定z=56），行重为1，它的值表示该矩阵的循环位移偏移量，也是第一行中元素1所处的列的位置，若H(i,j)值为0，表示是全0矩阵。
- 给定 $H_block{18×36}$ 矩阵，根据给定规则求解校验矩阵 $H{1008, 2016}$ 。
- 根据校验矩阵直接进行编码：利用输入信息比特序列$s$和校验矩阵$H$求得校验比特序列 $p$ ， $x = \lbrack p\ s\rbrack$ 即为编码序列。
- 编码序列采用BPSK调制并通过AWGN信道添加噪声。
- 完成以下四种译码算法的MATLAB代码实现
- 和积算法（Sum-Product, SP）
- 最小和算法（Min-Sum, MS）
- 归一化最小和算法(Normalized Min-Sum, NMS)
- 偏置最小和算法(Offset Min-Sum, OMS)
- 对于归一化最小和算法和偏置最小和算法，选定一个Eb/N0（1dB附近），仿真BER得出α和β最佳值。α和β取值范围均为0到1，步进为0.1，α和β的BER曲线各画一张图。
- 仿真得出Eb/N0在-1dB到2dB（步进0.5dB）时四种译码算法的BER和FER，两个修正最小和算法的α和β都取最佳值。应画BER一张图，FER一张图，每一张图都包含四种算法。
2.1. LDPC编码算法
LDPC码通常由校验矩阵H进行定义。尽管线性分组码可以使用传统的生成矩阵进行编码，但是要通过H求解G在实现上较为困难，因此根据校验矩阵直接进行编码。本次作业采用系统码，且校验矩阵H可分为两部分$$H=\left[ H_p|H_s \right] $$，其中$$H_p$$对应校验比特部分，对$$H_s$$应信息比特部分；编码序列x可分为$$x=\left[ p\ s \right] $$，其中$p$为校验比特序列，$s$为信息比特序列。本作业的编码部分有两种编码算法用以实现LDPC编码。
2.1.1. LDPC编码算法1
由于校验矩阵H的性质：
因此可采用基于H矩阵的通用LDPC编码方法：
由上式可知，本算法需要$H{p}$可逆，且上式p为1x(N-K)向量、s为1xK向量、$H{s}^{T}$为Kx(N-K)矩阵、$H_{p}^{T}$为Kx(N-K)矩阵，因此该算法需要$K(N - K) + (N - K)(N - K)$次比特乘法运算和$(K - 1)(N - K) + (N - K - 1)(N - K)$次比特加法运算。
2.1.2. LDPC编码算法2
将基于H矩阵的通用LDPC编码方法进行分解，具体步骤如下：
1) 首先利用信息序列s计算中间结果
2) 利用编码序列x计算校验比特
故可得：
此算法易错处为：上式中运算过程是应由右侧值决定左侧值，即只有当右侧表达式中各变量的值被赋过一次值时才能对左侧值赋值。例如：
- 当i=1时，w1已知，p1=w1
- 当i=2时，w2已知，但p(mb-1)*z并没有值，因此此时还不能立刻求p2，而是应该求p(z+1)=w(z+1) 模二求和 p(1)。即：-当i=(z+1)时，w(z+1)已知，p(1)已知，p(z+1)=w(z+1) 模二求和 p(1)
- ...
3) 组合p与s
将校验比特序列p和信息比特序列s组合即得到编码序列x。
4) 检验由于编码过程利用的是
因此，得到编码序列x后需要对x进行检验，保证其满足上式即表明完成了正确编码工作。
2.1.3. 两种算法比较
两种算法在时间复杂度上的区别如图所示（$\left( H_{p}^{T} \right)^{- 1}$可以离线求解，所以不占用复杂度）， 两种编码算法的时间复杂度对比
算法1适用于$H_{p}^{T}$可逆的通用H矩阵，而算法2适用于本次仿真作业特定的H矩阵。由时间复杂度的比较，本次仿真作业采用时间复杂度低的算法2。
2.2. 调制过信道
本次仿真作业采用实数BPSK调制并AWGN信道下进行误码率和误帧率的性能仿真。
2.2.1. BPSK
实数BPSK调制即将比特0映射为符号1，比特1映射为符号-1，变换式为：
其中d为BPSK调制后的符号（序列），x为编码比特（序列）。
在实数BPSK中，信噪比SNR和Eb/N0的转换关系为：
2.2.2. AWGN
BPSK调制后的信号通过AWGN信道将会对信号添加噪声，matlab中给信号添加AWGN噪声方式为：
其中，d为BPSK调制后的符号（序列），n为噪声信号（序列），y为通过AWGN信道后的信号值（序列）。噪声n在matlab中表示为：
其中$\sigma^{2}$为噪声功率，计算方法为：
2.2.3. 初始置信度
信息序列经过编码、BPSK调制，在AWGN信道下接收信号y的对数似然比（LLR，或初始置信度）为：
其中，
2.3. LDPC解码算法
本次仿真实验要求采用四种译码算法，分别为：
- 和积算法（Sum-Product, SP）
- 最小和算法（Min-Sum, MS）
- 归一化最小和算法(Normalized Min-Sum, NMS)
- 偏置最小和算法(Offset Min-Sum, OMS)
⛄ 部分代码
clear all
close all
clc
%% 预定义变量
N = 2016;
K = 1008;
R = K/N;
%% 添加工作路径
addpath('Encoder')
addpath('Decoder')
%% H矩阵生成
[ H, Hp, Hs ] = HxMatrixGen();
%% 仿真
Eb_N0_dB = 1.5;
alpha = 0:0.1:1;
BER = zeros(1, length(alpha));
for alpha_i = 1:1:length(alpha)
disp(['alpha = ' num2str(alpha(alpha_i)) ' is simulating...']);
% 设定停止条件
if Eb_N0_dB <= 1
maxErrorBlocks = 50;
else
maxErrorBlocks = 3;
end
% 设定译码算法最大迭代次数
iterMax = 30;
%设定每个信噪比下最大仿真帧个数
maxBlocks = 10^6;
% 四种算法的总误码数ErrorBits 和 总误帧数ErrorBlocks 和 所有帧的总循环数blocks 在每个Eb/N0仿真前清零
ErrorBits_NMS = 0;
ErrorBlocks_NMS = 0;
blocks_NMS = 0;
for i = 1:1:maxBlocks
% 算法2编码（s --> x）
s = randi([0, 1], 1, 1008);
x = Encoder2(Hs, Hp, s);
if sum(mod(H*(x'), 2)) > 0
sprintf('> the '+ num2str(i) + ' th encoding is not right');
end
% BPSK调制
d = 1 - 2.*x;
% AWGN
SNR_dB = Eb_N0_dB + 10*log10(R) - 10*log10(1/2);
SNR_linear = 10^(SNR_dB/10);
sigma = sqrt(1/SNR_linear);
y = d + sigma*randn(size(d)); % 加噪声
% 译码端接收
LLR_y = 2*y/(sigma^2);
% NMS译码
v_NMS = LDPCDecoder_NMS( H, LLR_y, alpha(alpha_i), iterMax );
%误比特数、误帧数统计
errorbits_NMS = sum(s ~= v_NMS);
ErrorBits_NMS = ErrorBits_NMS + errorbits_NMS;
blocks_NMS = blocks_NMS + 1;
if errorbits_NMS ~= 0
ErrorBlocks_NMS = ErrorBlocks_NMS + 1;
end
if ErrorBlocks_NMS > maxErrorBlocks
break;
end
end
BER(1, alpha_i) = ErrorBits_NMS/(K * blocks_NMS);
end
% 绘制BER
xlswrite('./BERforFindBestAlpha.xlsx', BER);
semilogy(alpha, BER, 'K-^', 'LineWidth', 1.0, 'MarkerSize', 5); % 三角marker 黑线
xlabel('\alpha'); ylabel('BER')
⛄ 运行结果
⛄ 参考文献
[1] 魏瑶. 准循环LDPC码的编译码技术研究与MATLAB仿真[D]. 河北大学, 2014.
[2] 陈洪雨. 跳频通信系统中LDPC编译码算法及FPGA的实现[D]. 哈尔滨工程大学, 2016.
[3] 张小华, 彭首峰. 基于Matlab的LDPC码研究及实现[J]. 科技传播, 2012(15):2.
