基于Matlab的LDPC编解码算法实现及LDPC码性能测试

最新推荐文章于 2024-06-21 19:15:59 发布

天天Matlab代码科研顾问

最新推荐文章于 2024-06-21 19:15:59 发布

阅读量2.3k

点赞数 22

分类专栏：信号处理文章标签： matlab 算法开发语言

本文链接：https://blog.csdn.net/m0_57702748/article/details/137714611

版权

信号处理专栏收录该内容

370 篇文章 57 订阅

订阅专栏

✅作者简介：热爱科研的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代码实现
1. 和积算法（Sum-Product, SP）
2. 最小和算法（Min-Sum, MS）
3. 归一化最小和算法(Normalized Min-Sum, NMS)
4. 偏置最小和算法(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)

2.3.1. 和积算法（SP）

算法流程为：

第一步：初始化

每个比特从信道中接收到的信道信息作为初始置信度，也即变量节点的置信度，其中表示第i个比特从信道接收到的LLR。和积算法初始化步骤的二分图为：

置信度传递

1) 每个变量节点要将自己的置信度传递给具有校验关系的校验节点。第i个变量节点传递给第j个校验节点的置信度计算方法

其中$d_{i}^{v}$表示第i个变量节点的度。和积算法变量节点传递置信度给校验节点步骤的二分图为：

2) 校验节点根据变量节点传递来的置信度，更新每一个变量节点的置信度。第j个校验节点为第i个变量节点更新的置信度为其中$d_{j}^{c}$表示第j个校验节点的度。和积算法校验节点传递置信度给变量节点步骤的二分图为：

判决

1) 置信度更新完毕后做硬判决

2) 如果所有的校验关系得到满足，即$H{\widehat{x}}^{T} = 0$，其中

，则译码结束，取出$\widehat{s} = \left( {\widehat{x}}{N - K + 1},{\widehat{x}}{N - K + 2},\cdots,{\widehat{x}}{N} \right)$即为信息位s的估计；如果不成立，则返回迭代的第二步继续迭代，直至达到迭代次数的上限。

2.3.2. 最小和算法（MS）

因为和积运算的积运算部分需要进行的反双曲正切连乘等非常复杂的运算，在硬件实现中通常只能利用查表等相对繁琐的操作来完成，因此需要进行改进。最小和算法与和积算法大同小异，区别在于积算法不再用

2.3.3. 归一化最小和算法（NMS）

归一化最小和算法建立在最小和算法之上，在计算校验节点为变量节点更新的置信度时引入归一化系数α，以减小Min-Sum译码消息的幅度α的取值在0到1之间，这种算法中通过引入α减小了译码消息的幅度，取值合适时性能逼近和积算法（为了简化要求，本次作业α取值精确到小数点后一位即可）。

2.3.4. 偏置最小和算法（OMS）

偏置最小和算法也是建立在最小和算法之上，在计算校验节点为变量节点更新的置信度时引入偏置量β$\min{k \in { 1,2\cdots d{j}^{c}}\backslash i}(\left| v_{k \rightarrow j} \right|)$ 小于β时，我们将其置为零，也就是消除了其对于下一步信息节点更新中的贡献。β的取值在0到1之间，这种算法使得过小的校验信息不再传递到变量节点中，取值合适时性能逼近和积算法。

3. 仿真流程与结果分析

用MATLAB平台搭建一个简单的仿真系统，验证LDPC编解码的性能。

在本次仿真作业中，首先在Eb/N0为1dB附近时根据BER求出最佳α和β取值；接下来，在2个修正最小和算法中使用最佳的α和β取值，绘制Eb/N0在-1dB到2dB（步进0.5dB）时四种译码算法的BER和FER。

在本次作业提交项目中，编码和解码算法分别位于"Encoder"和"Decoder"中，求解最佳α和β取值的主程序及四种算法性能比较的主程序位于" LDPC_Simulation"文件夹中，下表格说明了各程序功能。

📣 部分代码

clear allclose allclc%% 预定义变量N = 2016;K = 1008;R = K/N;%% 添加工作路径addpath('Encoder')addpath('Decoder')%% H矩阵生成[ H, Hp, Hs ] = HxMatrixGen();%% 仿真Eb_N0_dB = 1.0;beta = 0:0.1:1;BER = zeros(1, length(beta));for beta_i = 1:1:length(beta)        disp(['beta = ' num2str(beta(beta_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_OMS = 0;     ErrorBlocks_OMS = 0;    blocks_OMS = 0;        for i = 1:1:maxBlocks        % 算法2编码（s --> x）        % disp(['    the ' num2str(i) '-th frame of OMS decoding based on beta = ' num2str(beta(beta_i))]);                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_OMS = LDPCDecoder_OMS( H, LLR_y, beta(beta_i), iterMax );        %误比特数、误帧数统计        errorbits_OMS = sum(s ~= v_OMS);        ErrorBits_OMS = ErrorBits_OMS + errorbits_OMS;        blocks_OMS = blocks_OMS + 1;                if errorbits_OMS ~= 0            ErrorBlocks_OMS = ErrorBlocks_OMS + 1;        end        if ErrorBlocks_OMS > maxErrorBlocks            break;        end    end    BER(1, beta_i) = ErrorBits_OMS/(K * blocks_OMS);end% 绘制BERxlswrite('./BERforFindBestBeta.xlsx', BER);semilogy(beta, BER, 'K-^', 'LineWidth', 1.0, 'MarkerSize', 5); % 三角marker 黑线xlabel('\beta'); ylabel('BER');