当前位置：首页 > news >正文

LDPC 码 BP 算法性能研究

news 2025/9/17 9:42:12

1. 要点

来源：NASA 技术报告 + IEEE 802.16e 标准矩阵
算法：经典 Belief-Propagation（BP）迭代解码 + 最小和（MS）变体
性能：BER vs Eb/N0、迭代次数、码长/码率影响、BP vs MS 对比
输出：BER 曲线图、迭代收敛图、误码位置分布

2. 结构

LDPC_BP_Study/
├── main_bp_ber.m          % 一键 BER 曲线
├── ldpc_encode.m          % 系统编码（H→G）
├── bp_decode.m            % BP 迭代解码
├── ms_decode.m            % 最小和变体
├── ber_simulation.m       % BER 蒙特卡洛
├── plot_ber_curve.m       % 性能可视化
└── param/└── ieee80216e_H.mat   % 标准 H 矩阵（可换）

3. 核心

① BP 迭代解码（bp_decode.m）

function [bits_hat, iters] = bp_decode(llr_ch, H, maxIter)
% 经典 BP（和积算法）迭代解码
% llr_ch: 信道初始 LLR
% H: 稀疏校验矩阵
% maxIter: 最大迭代
[m, n] = size(H);

llr_v = llr_ch;  % 变量节点 LLR

iter = 0;
for iter = 1:maxIter

% 1. 变量→校验 消息

for j = 1:n

cols = find(H(:,j));

for i = cols'

rows = find(H(i,:));

rows(rows==j) = [];  % 排除自身

msg_v2c(i,j) = llr_v(j) + sum(msg_c2v(i,rows));

end

end
% 2. 校验→变量 消息（tanh 规则）
for i = 1:mrows = find(H(i,:));for j = rows&#39;cols = find(H(:,j));cols(cols==i) = [];S = sign(msg_v2c(cols,j));A = atanh(tanh(abs(msg_v2c(cols,j))/2));msg_c2v(i,j) = 2 * atanh(prod(S) * tanh(sum(A)/2));end
end% 3. 后验 LLR
llr_post = llr_v + sum(msg_c2v, 1);
bits_hat = llr_post &lt; 0;% 4. 早停判据
if all(mod(H * bits_hat&#39;, 2) == 0)break;
end

end

end

② 最小和变体（ms_decode.m）

function [bits_hat, iters] = ms_decode(llr_ch, H, maxIter)
% 最小和（Min-Sum）简化 BP
[m, n] = size(H);
llr_v = llr_ch;
for iter = 1:maxIter% 变量→校验for j = 1:ncols = find(H(:,j));for i = cols'rows = find(H(i,:)); rows(rows==j) = [];msg_v2c(i,j) = llr_v(j) + sum(msg_c2v(i,rows));endend
% 校验→变量（最小和近似）
for i = 1:mrows = find(H(i,:));for j = rows&#39;cols = find(H(:,j)); cols(cols==i) = [];abs_msg = abs(msg_v2c(cols,j));min_val = min(abs_msg);sign_prod = prod(sign(msg_v2c(cols,j)));msg_c2v(i,j) = sign_prod * min_val;end
endllr_post = llr_v + sum(msg_c2v, 1);
bits_hat = llr_post &lt; 0;if all(mod(H * bits_hat&#39;, 2) == 0)break;
end

end

end

③ BER 蒙特卡洛（ber_simulation.m）

function [ber, iters] = ber_simulation(EbN0_dB, H, maxIter, maxErr)
% BER vs Eb/N0 蒙特卡洛
EbN0 = 10^(EbN0_dB/10);
sigma = sqrt(1/(2*EbN0));
[n, k] = size(H);
ber = 0; err = 0; total = 0; iters = [];
while err < maxErr

msg     = randi([0 1], k, 1);           % 信息比特

codeword = ldpc_encode(msg, H);         % 系统编码

rx       = 1 - 2codeword + sigmarandn(n, 1); % AWGN

llr      = 2 * rx / (sigma^2);          % LLR
[dec, iter] = bp_decode(llr, H, maxIter);
err = err + sum(dec ~= codeword);
total = total + k;
iters = [iters, iter];

end

ber = err / total;

end

4. 运行（main_bp_ber.m）

clear; clc; close all; addpath('.');
%% 1. 加载标准 H 矩阵（IEEE 802.16e 0.75 码率）

load param/ieee80216e_H.mat;  % H: 576×432

maxIter = 20; maxErr = 500;   % 蒙特卡洛停止条件

EbN0_dB = 0:0.5:5;            % Eb/N0 范围
%% 2. BER 仿真（BP vs MS）

ber_bp = zeros(size(EbN0_dB));

ber_ms = zeros(size(EbN0_dB));

for k = 1:length(EbN0_dB)

fprintf('Eb/N0 = %.1f dB\n', EbN0_dB(k));

[ber_bp(k), ~] = ber_simulation(EbN0_dB(k), H, maxIter, maxErr);

[ber_ms(k), ~] = ber_simulation(EbN0_dB(k), H, maxIter, maxErr, 'MS');

end
%% 3. 性能对比

plot_ber_curve(EbN0_dB, ber_bp, ber_ms);