1.1BCH编码(附MATLAB仿真)

1.1BCH编码(附MATLAB仿真)

因工作需要看了很多关于BCH码的资料,但因为有很多的专业名词,所以看了很多天,综合了多方面的资料才大概有些明白BCH编解码的流程,下面也是我的学习笔记,方便以后查看,也可以方便大家一起学习,如有错误之处,还请一起交流讨论,随时修改补充。

写在前面:

只关心具体的流程和实现方法,不太关注为什么要这样做,底层逻辑什么的。

基础知识:

①BCH码是一种纠错码、线性分组码、循环码。

②需要传输信息位数:k

③纠错能力:t

④总码长(信息位+监督位):n

⑤n的长度满足n=2^m – 1时生成的为本原BCH码;n的长度为2^m – 1的因子时为非本原BCH码

(如n=15,n=31,n=63时为本原BCH码;n=21(可被63整除)等时为非本原BCH码)

⑥此外还有加长BCH码和缩短BCH码。

⑦具体的BCH码通常用BCH(n,k)码来表示。

加长BCH码和缩短BCH码:

因为本原BCH码和非本原BCH码要求了n的长度,但很多情况下我们想要的码长并不满足n=2^m – 1或其因子。这时候就需要加长BCH码和缩短BCH码。

(1)缩短BCH码

BCH(50,32)码是扩展域GF(2^6)上BCH(63,45)码的缩短码。BCH(50,32)码和BCH(63,45)码的区别与联系:

①两者纠错能力相同,生成多项式相同。

②缩短码的实现只需要在编译码时在高位上补0,从k = 32凑到k = 45即可。

(2)加长BCH码

在本原BCH码或非本原BCH码的生成多项式中乘因式(x+1),可以得到加长BCH码(n+1,k),加了一个校验位。

编码流程(只关心如何仿真实现,不关心具体流程可以跳过):

在BCH编码时,有三个多项式需要分清楚:本原多项式、极小多项式、生成多项式。编码时首先已知的是本原多项式,其次通过本原多项式算出极小多项式,最后用极小多项式得到最终的生成多项式,已知生成多项式就相当于已经编码成功了!

资料里获取BCH码的生成多项式的步骤如下:

下面用自己的话展开说说:

第一步

确定扩展域GF(P^m)及其本原多项式p(x)。

对于小白来说,不用关心扩展域是什么鬼,只需要知道,在GF(P^m)中只含有P^m个元素:0和a^i,其中i = 0,1,2,…,P^m – 2,(P=2时为二进制,我只关心了二进制)。并且通过本原多项式p(x)可以得到一个映射关系,举个具体的例子来看:

选择本原多项式为p(x)=1+x+x^3的扩展域GF(2^3),用扩展域的元素a来定义p(x)的根,也就是1+a+a^3=0 ——> a^3=1+a(二进制里,加就是减),也就是说,a^3可以表示为更低价a项的加权和。通过这个关系:a^3 = 1+a,可以将扩展域GF(2^3)中的八个元素映射为(本质上就是利用这一关系,把八个元素换算成最高幂次为2的多项式):

0 = 0

a^0 = 1

a^1 = a

a^2 = a^2

a^3 = a + 1

a^4 = a(a+1) = a^2 + a

a^5 = a^2 (a+1) = a^3 + a^2 = a^2 + a + 1

a^6 = a^3 (a+1) = (a+1)^2 = a^2 + a + a + 1 = a^2 + 1

则可知n=P^m – 1,在实例中也就是m = 3,P = 2,n = 7(总码长)。

第二步

因式分解x^n + 1 = f1(x)f2(x)…fm(x),则可获取极小多项式为f1(x),f2(x),…,fm(x)。

在上述实例中也就是对x^7+1因式分解得到:x^7 + 1 = (x+1)(x^3 + x + 1)(x^3 + x^2 + 1),也就是极小多项式有x+1、x^3 + x + 1、x^3 + x^2 + 1三种情况。那如何把极小多项式和a^i相对应呢?就是说我怎么知道哪个是a^i的极小多项式呢?下一步我们来为a^i分配极小多项式,fi(x)表示它为a^i的极小多项式。

第三步:

将fi(x)中代入a^i,其中i = 0,1,2,…,P^m – 2,看a^i是哪个fi(x)的根,fi(x)就是a^i的极小多项式(记得可以利用第一步中的映射关系哦)。还是用例子说话:

举例把a^3代入x+1中,由第一步可知a^3 = 1+a:

f(a^3) = 1+a+1 = a不等于0,不行!

再把a^3代入x^3 + x + 1中,得到:

f(a^3) = (a+1)^3+a+1+1 = (a^2 + a + a + 1)(a + 1) + a = a^3 + a^2 + a + 1 + a = a + 1 + a^2 + 1 = a^2 + a不等于0,不行!

再把a^3代入x^3 + x^2 + 1中,得到:

f(a^3) = (a+1)^3+(a+1)^2+1 = (a^2 + a + a + 1)(a + 1) + a^2 + 1 + a + a + 1 = a^3 + a^2 + a + 1 + a^2 = a + 1 + a + 1 = 0,行!a^3是x^3 + x^2 + 1的根!所以我们给a^3分配的极小多项式是f3(x) = x^3 + x^2 + 1。

同理可得:

f0(x) = x + 1

f1(x) = f2(x) = f4(x) = x^3 + x + 1

f3(x) = f5(x) = f6(x) = x^3 + x^2 + 1

第四步:

然后,我们就可以确定想要纠错几个比特了,假设我们想纠错t bit,则对应的生成多项式就是:g(x) = LCM[f1(x),f2(x),…,f2t(x)]。

最后,我们可以知道:deg g(x) = n – k,即:监督位位数 = g(x)的最高项的幂次。

BCH编码实现(matlab实现):

在上一部分中说了如何从本原多项式到生成多项式的过程,但是其实一般不需要计算获取生成多项式g(x),前人已经把找到的g(x)列成表,因此可以用查表法找到对应的生成多项式。一个特定的n,k和t对应一个g(x),n<=127的二进制本原BCH码生成多项式系数见课本《通信原理》(第七版)(樊昌信、曹丽娜 编著)P353,P354有部分二进制非本原BCH码生成多项式系数,表中给出的生成多项式系数用八进制数字列出,本文最后附录了该课本介绍BCH码的部分,包含表格。

所以想MATLAB仿真实现BCH编码需要:

第一步:确定你需要的n,k和t,对应附录中的生成多项式系数表格,得到生成多项式系数。

第二步:用n,k,t,g(x)系数代入如下代码中仿真。仿真实现了BCH编解码过程。

%-------------------------------------main------------------------------------------%
%% BCH编解码
%% 参数
n = 63;                                                                    % 信息比特+监督比特
t = 2;                                                                     % 纠错能力
k = 51;                                                                    % 要求的信息bit数
k0 = 28;                                                                   % 实际的信息bit数
gx = 12471;                                                                % 生成多项式系数(八进制)
error = 2;                                                                 % 经过信道后的错码bit数
%% 发送序列data与编码endata
disp(['码长 k0 = ',num2str(k0)])
data0 = zeros(1,k0);                                                       % 待编码信息序列(实际长度k0,若k0=k则为本原或非本原BCH码;若k0<k则为缩短BCH码,需要补长到k才可用本原或非本原多项式系数)
if k0 < k
    data = [zeros(1,k-k0) data0];
    data = data0;
endata = BCH_encode1(data,gx,k0);
%% 信道
noi_index = randperm(length(endata));
noi_index = noi_index(1:error);
noi = zeros(1,length(endata));
noi(noi_index) = 1;
noidata = mod(endata+noi,2);
disp(['发送码字为',num2str(endata)])
disp(['接收码字为',num2str(noidata)])
enum=0;
for i=1:length(endata)
    if endata(i)~=noidata(i)
        enum=enum+1;
disp(['通过信道后出错 ',num2str(enum),' 位'])
%% 接收端解码
[dedata,YN] = BCH_decode(noidata,n,k,k0,t);
if YN == 1
    disp('解码后存在误码');
    disp('解码正确');
YN1 = xor(dedata,data0);
if sum(YN1) ~= 0
    disp('实际上解码后存在误码');
    disp('实际上解码正确');
%---------------------------------BCH编码-----------------------------------------------%
function encode = BCH_encode1(code,gx,k0)
bin = oct2bin(gx);
while 0 == bin(1)    
    bin = bin(2:length(bin));
x1=zeros(1,length(bin));                                                    
x1(1)=1;
c1=conv(x1,code);                                                          
[~,r]=deconv(c1,bin);
r=mod(r,2);                                                               
encode = mod(c1+r,2);                                                      
encode = encode(length(code)-k0+1:end);
%----------------------------八进制转二进制------------------------------------------%
function bin = oct2bin(coef)
tmpData = [0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1; 0 0 0];
bin = [];
while coef ~= 0
    tmp = mod(coef,10);
    if  0 == tmp
        tmp = 8;
    bin = [tmpData(tmp,:) bin];
    coef = floor(coef/10);
%--------------------------------BCH解码-------------------------------------------%
function [decode1,YN] = BCH_decode(encode,n,k,k0,t)
% n = 31;                                                                    % 码长(信息码元+监督码元)
% t = 3;                                                                     % 纠错能力
% k = 16;
encode = [zeros(1,k-k0) encode];
m = 0;
while(2^m-1~=n && m<20)                                                      % 计算m
    m = m + 1;
a = gf(2,m);                                                               
s = a + a;                                                                   % 伴随式
for i = 1:2*t
    s(i) = a + a;
    for j = 1:n
        s(i) = s(i) + encode(j)*a^((n-j)*i);
for r = t:-1:1
    A = a + a;
    for i = 1:r
        for j = 1:r
            A(i,j) = s(r+i-j);
    if det(A) ~= 0 
        break;                                                             
d = rank(A);                                                               
B = a + a;
for i = 1:d
    B(i) = s(d+i);
if A == a + a                                                              
    decode = encode;
    E = zeros(1,n);
    sigma = A\B';                                                        
    E = zeros(1,n);
    x = a + a;
    ki = 1;
    for i = 1:n                                                              
        h = a^0;
        for j = 1:d
            h = h + sigma(j)*a^(i*j);
        if h == a + a
            x(k) = a^(n-i);
            E(i) = 1;                                                        
            ki = ki + 1;
    decode = mod(E+encode,2);                                                
decode1 = decode(k-k0+1:k); 
%% 解码后是否仍有误码(检错)
s = a + a;                                                                 
for i = 1:2*t
    s(i) = a + a;
    for j = 1:n
        s(i) = s(i) + decode(j)*a^((n-j)*i);
for r = t:-1:1                                                             
    A = a + a;
    for i = 1:r
        for j = 1:r
            A(i,j) = s(r+i-j);
    if det(A)~=0 
        break;                                                             
d = rank(A);                                                               
if d ~= 0
    YN = 1;
    YN = 0;
end

%------ 2023.3.30代码修改更新start ------%

去年编写代码的时候有一些疏漏:

①BCH译码模块:解码后的检错for i = 1:2*t和for r = t:-1:1两行一开始把t具体话成了5,是我没有看清楚,导致有时译码错误,现在已经直接在上面一版代码中修改过来了;

②上面一版代码还有一个问题(是一位读者发现的耶),在码长较短时,比如n=31或127时可以完美的编译码,但是码长较长,如n=1023时,就会有错误!经过排查,发现算法没有问题,问题出现在了bch编码模块的deconv函数上,deconv用来求余数,求完余数后还会二进制。位数少的时候,deconv求出来的余数有奇有偶,变成二进制有0有1,位数多了,求出来的余数都特别大,比如在n=1023时求出的余数量级在1e63,导致bch编码的校验位全为零(但是我也不知道为什么这种情况下deconv求出来的余数为什么这么大呀?),所以我自己编写了一个二进制除法取余的函数,下面附了二进制除法取余的函数bin_division.m和对bch编码函数的更改,使用此bch编码函数后可实现较长码字n=1023时的编译码:

%------------------------------二进制除法取余---------------------------------------%
function [q,r] = bin_division(A,B)
while 0 == A(1)    
    A = A(2:length(A));
while 0 == B(1)    
    B = B(2:length(B));
NA = length(A); % A长
NB = length(B); % B短
loc1 = 1;
loc2 = loc1 + NB - 1;
q = zeros(1,NA);
r = [];
en = 0;
while loc2 <= NA
    q(loc2) = 1;
    seq = A(loc1:loc2);
    Nr = length(r);
    t = NB - Nr;
    now = [r seq(end-t+1:end)];
    r = xor(now,B);
    if sum(r) == 0
        r = [];
        en = 1;
        while 0 == r(1)    
            r = r(2:length(r));
    Nr = length(r);
    loc2q = loc2;
    loc1 = loc2 - Nr + 1;
    loc2 = loc1 + NB - 1;
    if en == 1
        seq = A(loc1:loc2);
        deta = find(seq==1,1);
        loc1 = loc1 + deta - 1;
        loc2 = loc2 + deta - 1;
        en = 0;
r = [r A(loc2q+1:end)];
while 0 == q(1)    
    q = q(2:length(q));
%---------------------------------BCH编码-----------------------------------------------%
function encode = BCH_encode1(code,gx,k0)
bin = oct2bin(gx);
while 0 == bin(1)    
    bin = bin(2:length(bin));
x1=zeros(1,length(bin));                                                    
x1(1)=1;
c1=conv(x1,code);                                                          
[~,r] = bin_division(c1,bin);