第1页 / 共8页
第2页 / 共8页
第3页 / 共8页
第4页 / 共8页
第5页 / 共8页
第6页 / 共8页
第7页 / 共8页
第8页 / 共8页
74汉明码硬判决最大似然和积算法SPA仿真程序-hammingcodedecoding.doc
现代数字通信和编码理论
一、（7，4）汉明码误比特性能仿真曲线
分别采用硬判决、最大似然译码（MLD）、以及和积算法（SPA）三种译码
方法对（7，4）汉明码进行译码，得到该汉明码的译码性能曲线如图（1）所示。
图（1）硬判决，最大似然译码及和积译码的误码性能
为了节省仿真时间，对随机产生 8*105 个二进制信息进行编译码，仿真结果
表明，在加性高斯信道下，得到在误码率为 10-4 时 （7，4）汉明码的最大似然
译码较硬判决译码多出近 3dB 的编码增益，采用和积算法的迭代译码当迭代次
数为 100 时，误码性能非常接近最大似然译码，即迭代译码方式与最佳的译码方
式的性能相当。
二、译码原理概述
对任意正整数 m≥3，存在具有如下参数的汉明码：
码长： n=2m-1
信息符号数：k=2m-m-1
校验符号数：n-k=m
纠错能力：t=1（dmin=3）
本次实验中 n=7，k=4；即（7，4）汉明码。
假设在一个 AWGN 信道中采用（7，4）汉明码 C 作为差错控制编码，信道
的噪声均值为 0，双边功率谱密度为 0 / 2N 。采用 BPSK 调制，每个码元的发射
能量为单位能量。一个码字
v
(
,
v v
0
1
,...,
)n
v
1
在发射之前被映射到一个二进制序
列
x
(
,
x x
0
1
,...,
)n
x
1
， 其 中
x
l
2
v
l
， 当 lv =1 时 lx =1 ， 当 lv =0 时 lx =-1 ，
1
0
。令
n
1
l
r
(
,
r r
0
1
,...,
)n
r
1
为接收端匹配滤波器输出的软判决接收序列。对
于 0
r
， l
n
1
l
x
l
，其中 ln 是均值为 0，方差为 0 / 2N 的高斯随机变量。
n
l
对汉明码进行译码时分为硬判决和软判决两种方案。对于硬判决，接收解调
器的匹配滤波器将接收到的码字按照预先规定的判决门限直接将接收信号量化
为两个电平：0 和 1；然后计算每一组判决后的码字的校正子，得到对应的错误
图样，对接收码字进行纠错，如果超过其纠错能力则无法纠正其错误。对于（7，
4）汉明码可以纠正一组码字中的 1 个错误。
对于软判决译码可以得到比硬判决译码更好的误码性能，但是运算的复杂度
往往较高。软判决译码可以采用很多译码方案，本题的（7，4）汉明码采用最大
似然译码(ML)和和积算法(SPA)进行比较。最大似然译码是指：
(
P R C
\
c）
(r /
=
i
P
i
n
j=1
)ij
（1）
若能在 2k 个码字 C 中选择某一个 Ci 使式（1）取得最大，则这种译码规则称为
最大似然译码（MLD）。
P R C）为似然函数。对于 DMC 信道，最大似然译码
(
\
i
是使译码错误概率最小的一种最佳译码准则。为了简化译码运算复杂度，最大似
然译码可以等效为最小平方欧氏距离做为译码量度，即接收序列 r 被译为使平方
欧氏距离 2 ( , )
Ed r c 最小的码字 v 。而平方欧氏距离最小又可以等效为使下式最大：
( , )
m r v m r c
( , )
n
1
i
0
rc
i
i
（2）
将此求和项称为 r 与码字序列 c 之间的相关。因此软判决译码可以使用相关作为
译码度量，即将 r 译为使 ( , )
m r v 最大的码字 v 。
和积算法（SPA）是一种基于置信度传播的迭代译码，它也是一种逐符号的、
软输入软输出译码算法，基于奇偶校验和提高每个符号的可靠度。每个译码符号
的可靠度的量度可以采用其边缘后验概率、对数似然比或者对应的接受符号值。
每次译码迭代结束时得到的码符号的可靠度量度的计算结果被用于下一次迭代
的输入。译码迭代过程持续进行，直到满足某个特定条件（或准则），然后基于
码符号的可靠度量度的计算结果，作出硬判决。
附：源程序
%硬判决，采用校正子进行纠错，超过其纠错范围则译码失败；
clear;
Eb_N0dB=0:1:10;
Eb_N0=10.^(0.1*Eb_N0dB);
Es_N0=(4/7).*Eb_N0;%将 Eb_N0 转化为 Es_N0
m=3;n=2^m-1;k=n-m;
M=200000;
zero_matrix=[0 0 0];
%产生 U,表示所有可能的输入比特序列，对于（7，4）hamming 码是一个 16*4 的矩阵
for i=1:2^k;
for j=k:-1:1
if rem(i-1,2^(-j+k+1))>=2^(-j+k)
u(i,j)=1;
u(i,j)=0;
else
end
end
end
%产生生成矩阵 parmat 和校验矩阵 genmat
[parmat,genmat]=hammgen(m);
%产生编码后的码字空间 C
C=rem(u*genmat,2);
%msg1 是一个 M*k 个随机产生的信息位
msg1=randint(M*k,1,[0,1]);
%将 msg1 变成矩阵的形式，每行含有 k 个信息位
msg2=vec2mat(msg1,k);
%对信息序列进行编码
code1=encode(msg1,n,k,'hamming/binary');
code2=encode(msg2,n,k,'hamming/binary');%每行含有 n=7 个元素；
%对信息序列进行 BPSK 调制
bp_code1=2*code1-1;
bp_code2=2*code2-1;
%加入高斯白噪声，并采用硬判决译码
for i=1:1:length(Eb_N0dB);
noise_aver=0;noise_sigma=sqrt(1/(2* Es_N0(i)));
noise1=normrnd(noise_aver,noise_sigma,size(bp_code1));
noise2=vec2mat(noise1,n);
noise_code1=bp_code1+noise1;
noise_code2=bp_code2+noise2;
for j=1:1:length(noise_code1);
if noise_code1(j)>0,de_code1(j)=1;
else
end
de_code1(j)=0;
end
de_code2=vec2mat(de_code1,n);
%对 de_code2 采用校正子译码,分别对 M 行接收序列进行纠错；
for j=1:1:M
error=[0 0 0 0 0 0 0];
if(de_code2(j,:)*(parmat')==zero_matrix),de_code3(j,:)=de_code2(j,:);
else
S=rem((de_code2(j,:)*(parmat')),2);%caculate the 校正子
S_DE=bi2de(S);
parmat_de=bi2de(parmat');
index1=find(parmat_de==S_DE);
if(index1)%计算错误图样
error(index1)=1;
end
de_code3(j,:)=de_code2(j,:)+error;%，
end
end
for a=1:1:M;
for b=4:1:n;
in_decode(a,b-3)=de_code3(a,b);
end
end
[number(i),ratio(i)]=biterr(in_decode,msg2);
end
plot(Eb_N0dB,ratio,'k-');
hold on;
%软判决译码，采用最大似然译码准则译码；
for i=1:1:length(Eb_N0dB);
noise_aver=0;noise_sigma=sqrt(1/(2* Es_N0 (i)));
noise1=normrnd(noise_aver,noise_sigma,size(bp_code1));
noise2=vec2mat(noise1,n);
noise_code1=bp_code1+noise1;
noise_code2=bp_code2+noise2;
%对接收到的软序列，进行最大似然译码（等效于最小距离译码）；
for j=1:1:M
w=noise_code2(j,:)';
product=C*w;
index2=find(product==max(product));
de_code3(j,:)=C(index2(1),:);
end
for a=1:1:M;
for b=4:1:n;
in_decode(a,b-3)=de_code3(a,b);
end
end
[number(i),ratio(i)]=biterr(in_decode,msg2);
end
plot(Eb_N0dB,ratio,'k-o');
hold on;
%和积算法（SPA），迭代 100 次
for i=1:1:length(Eb_N0dB);
noise_aver=0;noise_sigma=sqrt(1/(2* Es_N00(i)));sigma2=noise_sigma*noise_sigma;
noise1=normrnd(noise_aver,noise_sigma,size(bp_code1));
noise2=vec2mat(noise1,n);
noise_code1=bp_code1+noise1;
noise_code2=bp_code2+noise2;
%转换成接收序列的先验概率信息，p_zero（）为每一位取 0 的概率；
for j=1:1:M;
for l=1:1:n;
p_zero(j,l)=1/(1+exp(2*noise_code2(j,l)/sigma2));
p_one(j,l)=1/(1+exp(-2*noise_code2(j,l)/sigma2));
end
end
%以 7 个信息位为一个单位，分别对 M*7 个码组依次进行迭代译码；
for j=1:1:M;%共 M 行 7 位的信息组
Imax=1;%最大迭代次数初始化；
for a=1:1:m;
for b=1:1:n;%第一次迭代，变量节点到校验节点传递的信息
if parmat(a,b)==0,q_zero(a,b)=0;q_one(a,b)=0;%初始化；
else q_zero(a,b)=p_zero(j,b);q_one(a,b)=p_one(j,b);
end
end
end
while Imax<=100
Imax=Imax+1;
for a=1:1:m; %校验节点传递到变量节点的信息；
w=parmat(a,:);%计算从第一个校验节点开始的信息传递；
index1=find(w);%找出在对应的第 a 行中检验矩阵为 1 的元素的索引；
for b=1:1:n;%计算校验节点传递到所有变量节点的信息；
if parmat(a,b)==0,S_zero(a,b)=0;S_one(a,b)=0;
%去掉要传向信息位的索引
else index2=index1;x=find(index2==b);index2(x)=[];
S_zero(a,b)=0;S_one(a,b)=0;%初始化
for x1=1:2^length(index2);%产生剩余节点可能对应取值的
for x2=length(index2):-1:1
if
rem(x1-1,2^(-x2+length(index2)+1))>=
2^(-x2+length(index2))
一个矩阵
v(x1,x2)=1;
else
v(x1,x2)=0;
end
end
end
%对剩余变量节点所有可能的取值进行判断；
for x1=1:2^length(index2);
check=0;P_zero=1;P_one=1;%清零和初始化；
for x2=1:1:length(index2);
check=rem(v(x1,x2)+check,2);
end
if check==0
%check=0 计算一种可能取值的概率
for x2=1:1:length(index2);
if
v(x1,x2)==0,
P_zero=P_zero*q_zero(a,index2(x2));
else P_zero=P_zero*q_one(a,index2(x2));
end
end
S_zero(a,b)=S_zero(a,b)+P_zero;
Else
%check=1 计算另一种可能取值的概率；
for x2=1:1:length(index2);
if
v(x1,x2)==0,
P_one=P_one*q_zero(a,index2(x2));
else P_one=P_one*q_one(a,index2(x2));
end
end
S_one(a,b)=S_one(a,b)+P_one;
end
end%对应 for x1=1:2^length(index2);
end%对应 else index2=index1;x=find(index2==b);index2(x)=[];
end%对应 for b=1:1:n;
end%对应 a=1:1:m;
%更新 Q 值，计算更新后变量节点到校验节点传递的信息；
for a=1:1:m;
for b=1:1:n;
if parmat(a,b)==0,q_zero(a,b)=0;q_one(a,b)=0;
else w=parmat(:,b);
index2=find(w);y=find(index2==a);index2(y)=[];
if index2,q_zero(a,b)=1;q_one(a,b)=1;%初始化
for l=1:1:length(index2);
q_zero(a,b)=q_zero(a,b)*S_zero(index2(l),b);
q_one(a,b)=q_one(a,b)*S_one(index2(l),b);
end
q_zero(a,b)=q_zero(a,b)*p_zero(j,b);
q_one(a,b)=q_one(a,b)*p_one(j,b);
else
end
q_zero(a,b)=p_zero(j,b);
q_one(a,b)=p_one(j,b);
end
end
end
%将 q_zero(a,b);q_one(a,b)标准化；
for a=1:1:m;
for b=1:1:n;
if q_zero(a,b)==0,sw_q_zero(a,b)=0;sw_q_one(a,b)=0;
else sw_q_zero(a,b)=q_zero(a,b)/(q_zero(a,b)+q_one(a,b));
sw_q_one(a,b)=q_one(a,b)/(q_zero(a,b)+q_one(a,b));
end
end
end
%为下一次迭代准备
q_zero=sw_q_zero;
q_one=sw_q_one;
%计算本次迭代的伪后验概率；
for b=1:1:n;
Pse_p_zero(b)=1;Pse_p_one(b)=1;%后验概率初始化；
w=parmat(:,b);index4=find(w);
for l=1:1:length(index4);
Pse_p_zero(b)=Pse_p_zero(b)*S_zero(index4(l),b);
Pse_p_one(b)=Pse_p_one(b)*S_one(index4(l),b);
end
Pse_p_zero(b)=Pse_p_zero(b)*p_zero(j,b);
Pse_p_one(b)=Pse_p_one(b)*p_one(j,b);
end
%伪后验概率归一化；
for b=1:1:n;
Nor_p_zero(b)=Pse_p_zero(b)/( Pse_p_zero(b)+Pse_p_one(b));
Nor_p_one(b)=Pse_p_one(b)/( Pse_p_zero(b)+Pse_p_one(b));
end
%基于本次迭代的伪后验概率，做译码判决；
for b=1:1:n;
if Nor_p_one(b)>0.5,de_code(j,b)=1;
else de_code(j,b)=0;
end
end
%计算校正子，若正确则停止迭代
z=de_code(j,:);sy=rem(z*(parmat'),2);
if sy==zero_matrix,break;
end
end
end
%去掉校验位，得到信息位并计算误码率
for a=1:1:M;
for b=4:1:n;
in_decode(a,b-3)=de_code(a,b);
end
end
[number(i),ratio(i)]=biterr(in_decode,msg2);
end
plot(Eb_N0dB,ratio,'k-v');
xlabel('EbN0(dB)');
ylabel('误比特率');
legend('硬判决','最大似然','SPA-迭代 100 次');
grid on;
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
