【BCH码2】BCH码简化的BM迭代译码原理详解及MATLAB实现(不使用MATLAB库函数)

2024-04-05 12:44:24

关注公号【逆向通信猿】更精彩！！！

校验矩阵 H \boldsymbol{H} H

零空间的定义：

H \boldsymbol H H的零空间或核记为
N ( H ) ≜ { x ∈ R n : H x = 0 } \mathcal{N}(\boldsymbol{H}) \triangleq\left\{\boldsymbol{x} \in \mathbb{R}^{n}: \boldsymbol{H} \boldsymbol{x}=\mathbf{0}\right\} N(H)≜{x∈Rn:Hx=0}

所以，BCH码空间是 H \boldsymbol H H的零空间。

译码原理

具体的原理推导不再赘述，请参考文献[1]，在此仅归纳总结出其核心思想。

译码思想

1 方程组的建立

根据伴随式 S i S_i Si与错误图样 e ( X ) e(X) e(X)的关系
S i = e ( α i ) S_i=\boldsymbol e(\alpha ^i) Si=e(αi)

得到错误位置 j 1 , j 2 , … , j v j_1,j_2,\ldots,j_v j1,j2,…,jv和 2 t 2t 2t个伴随式的关系方程组：

解方程组得到错误位置即可进行纠错译码，然而方程为非线性的，直接求解不可行，因此采用间接求解方法。

2 牛顿恒等式

定义

观察上面的式子，可以看到，利用完全平方展开式，这个构造能够将伴随式分量与错误位置多项式的系数联系起来，得到牛顿恒等式

通过牛顿恒等式解方程，得到错误位置多项式 σ ( X ) \sigma(X) σ(X)及其系数，进而得到 σ ( X ) \sigma(X) σ(X)的根，而根的倒数正是错误位置。由此得到译码过程如下。

3 译码过程

4 牛顿恒等式的迭代求解（BM算法核心）

译码过程（2）中牛顿恒等式的求解，可通过Berlekamp Massey迭代算法（BM算法）进行求解。

BM算法通过 2 t 2t 2t次迭代求解出错误位置多项式 σ ( X ) \sigma(X) σ(X)，有点类似归纳法。

第 k k k次得到一个次数最小的 σ ( k ) ( X ) \sigma^{(k)}(X) σ(k)(X)，满足前 k k k个方程。第 k + 1 k+1 k+1次，基于 σ ( k ) ( X ) \sigma^{(k)}(X) σ(k)(X)找到下一个次数最小的错误位置多项式 σ ( k + 1 ) ( X ) \sigma^{(k+1)}(X) σ(k+1)(X)，使其满足前 k + 1 k+1 k+1个方程。

首先，检查 σ ( k ) ( X ) \sigma^{(k)}(X) σ(k)(X)的系数是否满足第 k + 1 k+1 k+1个方程，如果满足，则
σ ( k + 1 ) ( X ) = σ ( k ) ( X ) \sigma^{(k+1)}(X)=\sigma^{(k)}(X) σ(k+1)(X)=σ(k)(X)

否则，对 σ ( k ) ( X ) \sigma^{(k)}(X) σ(k)(X)添加校正项得到 σ ( k + 1 ) ( X ) \sigma^{(k+1)}(X) σ(k+1)(X)。具体如下：

计算

译码算法描述

q q q元BCH码BM迭代算法的简化

二元BCH码BM迭代算法的简化

用一个实例还原BM迭代译码过程

以(15,5,3)BCH码为例，码长为15，纠错能力 t = 3 t=3 t=3。域生成多项式为 p ( X ) = X 4 + X + 1 p(X)=X^4+X+1 p(X)=X4+X+1（用次数表示为[1 0 0 1 1]，用八进制数表示为23，下文中统一用八进制表示多项式），码生成多项式用八进制数表示为2467。

生成的域 G F ( 2 4 ) GF(2^4) GF(24)如下表

首先进行初始化，如下

i i i的初始值为-1， i − l i = − 1 i-l_i=-1 i−li=−1，开始迭代：
① k = 1 k=1 k=1时，由于 d 0 = 1 ≠ 0 d_0=1\ne0 d0=1=0，所以计算校正项为 d 0 d − 1 − 1 X 0 − ( − 1 ) σ ( − 1 ) ( X ) = X d_0d_{-1}^{-1}X^{0-(-1)}\sigma^{(-1)}(X)=X d0d−1−1X0−(−1)σ(−1)(X)=X，所以 σ ( 1 ) ( X ) = 1 + X \sigma^{(1)}(X)=1+X σ(1)(X)=1+X，由此得到 d 1 = S 1 + 1 + σ 1 ( 1 ) S 1 = S 2 + 1 ⋅ S 1 = 0 d_1=S_{1+1}+\sigma_1^{(1)}S_1=S_2+1\cdot S_1=0 d1=S1+1+σ1(1)S1=S2+1⋅S1=0，所以 l 1 = 1 l_1=1 l1=1， 1 − l 1 = 0 1-l_1=0 1−l1=0， i = 0 i=0 i=0， i − l i = 0 i-l_i=0 i−li=0。
表格更新为

② k = 2 k=2 k=2时，由于 d 1 = 0 d_1=0 d1=0，所以 σ ( 2 ) ( X ) = σ ( 1 ) ( X ) = 1 + X \sigma^{(2)}(X)=\sigma^{(1)}(X)=1+X σ(2)(X)=σ(1)(X)=1+X，由此得到 d 2 = S 2 + 1 + σ 1 ( 2 ) S 2 = S 3 + 1 ⋅ S 2 = α 10 + 1 = α 5 d_2=S_{2+1}+\sigma_1^{(2)}S_2=S_3+1\cdot S_2=\alpha^{10}+1=\alpha^5 d2=S2+1+σ1(2)S2=S3+1⋅S2=α10+1=α5，所以 l 2 = 1 l_2=1 l2=1， 2 − l 2 = 1 2-l_2=1 2−l2=1， i = 0 i=0 i=0， i − l i = 0 i-l_i=0 i−li=0。
表格更新为

③ k = 3 k=3 k=3时，由于 d 2 = α 5 ≠ 0 d_2=\alpha^5\ne0 d2=α5=0，所以计算校正项为 d 2 d 0 − 1 X 2 − 0 σ ( 0 ) ( X ) = α 5 X 2 d_2d_{0}^{-1}X^{2-0}\sigma^{(0)}(X)=\alpha^5X^2 d2d0−1X2−0σ(0)(X)=α5X2，所以 σ ( 3 ) ( X ) = σ ( 2 ) ( X ) + α 5 X 2 σ ( 0 ) ( X ) = 1 + X + α 5 X 2 \sigma^{(3)}(X)=\sigma^{(2)}(X)+\alpha^5X^2\sigma^{(0)}(X)=1+X+\alpha^5X^2 σ(3)(X)=σ(2)(X)+α5X2σ(0)(X)=1+X+α5X2，由此得到 d 3 = S 3 + 1 + σ 1 ( 3 ) S 3 + σ 2 ( 3 ) S 2 = S 4 + 1 ⋅ S 3 + α 5 S 2 = 1 + α 10 + α 5 = 0 d_3=S_{3+1}+\sigma_1^{(3)}S_3+\sigma_2^{(3)}S_2=S_4+1\cdot S_3+\alpha^5S_2=1+\alpha^{10}+\alpha^5=0 d3=S3+1+σ1(3)S3+σ2(3)S2=S4+1⋅S3+α5S2=1+α10+α5=0，所以 l 3 = 2 l_3=2 l3=2， 3 − l 3 = 1 3-l_3=1 3−l3=1， i = 2 i=2 i=2， i − l i = 1 i-l_i=1 i−li=1。
表格更新为

④以此类推，最后得到 2 t 2t 2t次的迭代结果如下：

简化的BM算法的结果为：

仿真性能与分析

仿真结果

主程序代码

文件名：bch_codec.m

clear; close all; clc

%% (15,11,1) BCH码
% px=[1 0 0 1 1];  % 23本原多项式
% p = length(px)-1;
% [a,a_dec,log_a] = sub_gf_gen( px ); % 有限域
% gx=[1 0 0 1 1]; % 23
% t=1;

%% (15,7,2) BCH码
% gx=[1 1 1 0 1 0 0 0 1]; % 721
% t=2;

%% (31,26,1) BCH码
% px=[1 0 0 1 0 1];  % 45
% p = length(px)-1;
% [a,a_dec,log_a] = sub_gf_gen( px ); % 有限域
% gx=[1 0 0 1 0 1]; % 45
% t=1;
%% (31,21,2) BCH码
% gx=[1 1 1 0 1 1 0 1 0 0 1]; % 3551
% t=2;

%% (63,57,1) BCH码
px=[1 0 0 0 0 1 1];  % 103
p = length(px)-1;
[a,a_dec,log_a] = sub_gf_gen( px ); % 有限域
% gx=[1 0 0 0 0 1 1]; % 103
% t=1;

%% (63,51,2) BCH码
gx=[1 0 1 0 1 0 0 1 1 1 0 0 1]; % 12471
t=2;

m = length(px)-1; 
n=2^m-1; % 码长
pg = length(gx) -1; % 生成多项式次数
k=n-pg; % 信息长

codnum = 100000;
snr = 5:0.5:10;
snrlen = length(snr);
num_bef = zeros(1,snrlen);
num_aft = zeros(1,snrlen);
tic
parfor i = 1:snrlen
    SNRi = snr(i);
    for j = 1:codnum
%         rng('shuffle');
        u = randi([0,1],1,k);
        v = sub_encod_bch(u,gx,0);
        vn = awgn(1-2*v,SNRi,'measured');
        r = double(vn<0);
        num_bef(i) = num_bef(i)+biterr(r,v);
        [rdec,pos_err,flag]=sub_decod_bch(r,t,a,log_a,p);
        num_aft(i) = num_aft(i)+biterr(rdec,v);        
    end
    disp(['信噪比 ',num2str(SNRi),' 完成']);
end
toc
ber_bef = num_bef/n/codnum;
ber_aft = num_aft/n/codnum;

画图代码

文件名：plot_bch_dec.m

clear; clc; close all;

snr=5:0.5:10;
len = length(snr);
lw = 2;

%% 译码前误码率
load('ber_bef.mat'); 
figure;semilogy(snr,ber_bef,'Marker','x','color',[0.49 0.18 0.56],'LineWidth',lw);hold on;

%% (15,11,1)BCH码
load('ber_aft_bch_15_11_1.mat'); 
semilogy(snr,ber_aft,'Marker','s','color',[0 161/255 59/255],'LineWidth',2);hold on;

%% (15,7,2)BCH码
load('ber_aft_bch_15_7_2.mat'); 
semilogy(snr,ber_aft,'o-b','LineWidth',lw);hold on;

%% (31,26,1)BCH码
load('ber_aft_bch_31_26_1.mat'); 
semilogy(snr,ber_aft,'Marker','d','color',[0.85 0.33 0.10],'LineWidth',lw);hold on;

%% (31,21,2)BCH码
load('ber_aft_bch_31_21_2.mat'); 
semilogy(snr,ber_aft,'^-r','LineWidth',lw);hold on;

%% (63,57,1)BCH码
load('ber_aft_bch_63_57_1.mat'); 
semilogy(snr(4:end),ber_aft(4:end),'Marker','v','color',[0 0 0],'LineWidth',lw);hold on;

%% (63,51,2)BCH码
load('ber_aft_bch_63_51_2.mat'); 
semilogy(snr,ber_aft,'Marker','p','color',[0 0.4470 0.7410],'LineWidth',lw);hold on;

xlabel('SNR (dB)');ylabel('译码BER'); grid on;
legend('译码前','(15,11,1)','(15,7,2)','(31,26,1)','(31,21,2)','(63,57,1)','(63,51,2)');
set(gca,'FontSize',13,'Fontname', '微软雅黑');

部分子程序代码

有限域乘法子程序

文件名：sub_gf_mul.m

function [ z ] = sub_gf_mul( x, y, n )
% 伽罗华域元素乘法，参数含义：
%% input:
% x y: 表示两个元素的幂次，取值范围-inf或[0,n-1]
% n：2^m-1

%% output:
% z：表示乘积的幂次

if x>=0 && y>=0
    z = mod(x+y,n);
else
    z=-inf;
end

有限域加法子程序

文件名：sub_gf_add.m

function [ z ] = sub_gf_add( x, y, a, log_a )
% 伽罗华域元素加法，参数含义：
%% input:
% x y: 表示两个元素的幂次，取值范围-inf或[0,n-1]
% n：2^m-1

%% output:
% z：表示两个元素和的幂次，和为0时规定幂次为-inf，方便计算

if x<0 
    z=y;
elseif y<0
    z=x;
elseif x==y
    z=-inf;
else
    sum=mod(a(x+1,:)+a(y+1,:),2);
    sum_deci = bi2de(sum, 'left-msb');
    z = log_a(sum_deci);
end

所有代码下载

文件名	功能
bch_codec.m	主程序
sub_gf_gen.m	有限域生成子程序
sub_gf_add.m	有限域加法子程序
sub_gf_mul.m	有限域乘法子程序
sub_poly_div.m	多项式除法子程序
sub_encod_bch.m	BCH编码子程序
sub_decod_bch.m	BCH码简化的BM迭代译码子程序
plot_bch_dec.m	画图子程序
ber_aft_bch_15_7_2.mat	(15,7,2)码译码误码率数据文件
ber_aft_bch_15_11_1.mat	(15,11,1)码译码误码率数据文件
ber_aft_bch_31_26_1.mat	(31,26,1)码译码误码率数据文件
ber_aft_bch_31_21_2.mat	(31,21,2)码译码误码率数据文件
ber_aft_bch_63_57_1.mat	(63,57,1)码译码误码率数据文件
ber_aft_bch_63_51_2.mat	(63,51,2)码译码误码率数据文件

a)土豪请随意

直接点击以下链接，到博主资源中进行下载所有代码（带注释）和画图数据

所有代码下载：https://download.csdn.net/download/wlwdecs_dn/45111931

还可订阅《信道编码》专栏，查阅各子程序的详解

本原多项式生成有限域的原理及MATLAB实现

二元多项式除法电路原理及MATLAB实现

【BCH码1】系统BCH码编码原理及MATLAB实现（不使用MATLAB库函数）

b)支持博主的好心人请看这

关注公号【逆向通信猿】，口令：BCH，简单支持一波即可获取

c)伸手白嫖党

无

参考文献

[1]William E Ryan, Shu Lin. 信道编码：经典与现代[M]. 电子工业出版社, 2017.

码农公寓