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);