Viterbi译码的MATLAB仿真研究

BUPT

卷积码编码及Viterbi译码

班级：07114 学号：070422 姓名：吴希龙 指导老师：彭岳星

邮箱：FusionBupt@163.com

0

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

1. 序言

卷积码最早于1955年由Elias提出，稍后，1957年Wozencraft提出了一种有效地译码方法即序列译码。1963年Massey提出了一种性能稍差但是比较实用的门限译码方法，使得卷积码开始走向实用化。而后1967年Viterbi提出了最大似然译码算法，它对存储级数较小的卷积码很容易实现，被称作Viterbi译码算法，广泛的应用于现代通信中。

2. 卷积码编码及译码原理

2.1 卷积码编码原理

卷积码是一种性能优越的信道编码，它的编码器和解码器都比较易于实现，同时还具有较强的纠错能力，这使得它的使用越来越广泛。卷积码一般表示为(n,k,K)的形式，即将k各信息比特编码为n个比特的码组，K为编码约束长度，说明编码过程中相互约束的码段个数。卷积码编码后的n各码元不经与当前组的k个信息比特有关，还与前K-1个输入组的信息比特有关。编码过程中相互关联的码元有K*n个。R=k/n是编码效率。编码效率和约束长度是衡量卷积码的两个重要参数。典型的卷积码一般选n,k较小，但K值可取较大(>10)，以获得简单而高性能的卷积码。

卷积码的编码描述方式有很多种：冲激响应描述法、生成矩阵描述法、多项式乘积描述法、状态图描述，树图描述，网格图描述等。

2.1.1 卷积码解析表示法

卷积码的解析表示发大致可以分为离散卷积法，生成矩阵法，码多项式法。下面以离散卷积为例进行说明。

卷积码的编码器一般比较简单，为一个具有k个输入端，n个输出端，m级移位寄存器的有限状态有记忆系统。下图所示为(2,1,7)卷积码的编码器。

若输入序列为u=(u0u1u2u3……),

则对应两个码字序列𝐜①=(c0c1c2c3……)和𝐜②=(c0c1c2c3……)

相应的编码方程可写为𝐜①=𝐮∗𝐠①，𝐜②=𝐮∗𝐠②，c=(𝐜①,𝐜②)。 “∗” 符号表示卷积运算，𝐠①，𝐠②表示编码器的两个冲激响应，即编码器的输出可以由输入序列和编码器的两个冲击响应卷积而得到，故称为卷积码。这里的冲激响应指：当输入为[1 0 0 0 0 … … ]序列时，所观察到的两个输出序列值。由于上图K值为7，故冲激响应至

1

①①①①

②②②②

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

多可持续到第7位，可写为

𝐠①=[1 1 1 1 0 0 1]，𝐠②=[1 0 1 1 0 1 1]

然后将两个输出端的码字序列合并为一个码字序列为 𝐜=(c0c0c1c1c2c2……) 若输入信息序列为[1 1 0 1];

①②①②①②

则𝐜①=[1 0 0 1 0 1 0 1 0 1]，𝐜②=[1 1 1 1 1 0 1 1 1 1] c=[1 1 0 1 0 1 1 1 0 1 1 0 0 1 1 1 0 1 1 1]。

下图所示为(2,1,3)卷积码的编码器，也是本次课程设计所研究的卷积码编码器，由于其生成冲激响应分别为[1 1 1]和[1 0 1]，故被称为(7,5)码。

2.1.2 卷积码图形表示法

除了用解析法描述卷积码的编码外，还可以使用比较形象的图形法来表示卷积码。比较常用的有状态图法，网格图法和树图法。下面简介状态图法和网格图法。 状态图法：

由于卷积码编码器在下一时刻的输出取决于编码器的当前状态和下一时刻的输入，而编码器当前状态取决于编码器当前各移位寄存器的存储内容。称编码器当前各移位寄存器存储内容(0或)为编码器在该时刻的状态(此状态代表记忆以前的输入信息)。随着信息序列的不断输入，编码器不断从一个状态转移到另外一个状态，并且输出相应的编码序列。编码器的总可能状态数为2mk个。对(7,5)码的编码器来说，n=2，k=1，K=3，m=2。共有四个可能状态，其状态图如下：

图中四个方块表示状态，状态间的连线与箭头表示转移方向，连线上的数字表示是状态发生转移的到来比特，斜杠后的数字由一个状态到另一个状态转移时的输出码字。如当前状态为11，输入信息为0，则转移到01状态并输出01码字，若输入信息为1，则依然为11

2

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

状态并输出10码字。

网格图法： 网格图可以描述卷积码的状态随时间推移而转移的情况。该图纵坐标表示所有状态，横坐标表示时间。网格图在卷积码的概率译码，特别是Viterbi译码中非常重要，它综合了状态图法直观简单和树图法时序关系清晰的特点。如下图

图中实线表示输入0时所走分支，虚线表示输入1时所走分支，编码时只需从起始状态开始依次选择路线并读出输出即可。假设从a状态开始，输入为[1 0 1 1]，则可由图中读出输出为[11 10 10 01]。

2.2 卷积码译码原理

卷积码的译码方式主要有三种:

1).1963年Massey提出的门限译码，这是一种基于码代数结构的代数译码，类似于分组码中的大数逻辑译码。

2).1963年有Fano改进的序列译码，这是基于码的树状图结构的一种准最佳概率译码。 3).1967年Viterbi提出的Viterbi算法，基于码的网格图基础上的最大似然译码算法，是一种最佳概率译码。

其中，代数译码，利用编码本身的代数结构进行译码，不考虑信道本身的统计特性。该方法的硬件实现简单，但性能较差，其中具有典型意义的是门限译码。另一类是概率译码，这种译码通常建立在最大似然准则的基础上。由于计算是用到了信道的统计特性.因而提高了译码性能，但这种性能的提高是以增加硬件的复杂度为代价的。常用的概率译码方法有维特比译码和序列译码。维特比译码具有最佳性能，但硬件实现复杂；门限译码性能最差，但硬件简单；序列译码在性能和硬件方面介于维特比译码和门限译码之间。

2.2.1 Viterbi译码

卷积码概率译码的基本思路是:以接收码流为基础，逐个计算它与其他所有可能出现的、连续的网格图路径的距离，选出其中可能性最大的一条作为译码估值输出。概率最大在大多数场合可解释为距离最小，这种最小距离译码体现的正是最大似然的准则。卷积码的最大似然译码与分组码的最大似然译码在原理上是一样的，但实现方法上略有不同。主要区别在于:分组码是孤立地求解单个码组的相似度，而卷积码是求码字序列之间的相似度。基于网格图搜索的译码是实现最大似然判决的重要方法和途径。用格图描述时，由于路径的汇聚消除了树状图中的多余度，译码过程中只需考虑整个路径集合中那些使似然函数最大的路径。如果

3

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

在某一点上发现某条路径已不可能获得最大对数似然函数，就放弃这条路径，然后在剩下的“幸存”路径中重新选择路径。这样一直进行到最后第L级(L为发送序列的长度)。由于这种方法较早地丢弃了那些不可能的路径，从而减轻了译码的工作量，Viterbi译码正是基于这种想法。 对于(n, k, K )卷积码，其网格图中共2kL种状态。由网格图的前K-1条连续支路构成的路径互不相交，即最初2k_1条路径各不相同，当接收到第K条支路时，每条路径都有2条支路延伸到第K级上，而第K级上的每两条支路又都汇聚在一个节点上。在Viterbi译码算法中，把汇聚在每个节点上的两条路径的对数似然函数累加值进行比较，然后把具有较大对数似然函数累加值的路径保存下来，而丢弃另一条路径，经挑选后第K级只留下2K条幸存路径。选出的路径同它们的对数似然函数的累加值将一起被存储起来。由于每个节点引出两条支路，因此以后各级中路径的延伸都增大一倍，但比较它们的似然函数累加值后，丢弃一半，结果留存下来的路径总数保持常数。由此可见，上述译码过程中的基本操作是，“加-比-选”，即每级求出对数似然函数的累加值，然后两两比较后作出选择。有时会出现两条路径的对数似然函数累加值相等的情形，在这种情况下可以任意选择其中一条作为“幸存”路径。

卷积码的编码器从全零状态出发，最后又回到全零状态时所输出的码序列，称为结尾卷积码。因此，当序列发送完毕后，要在网格图的终结处加上（K-1）个己知的信息作为结束信息。在结束信息到来时，由于每一状态中只有与已知发送信息相符的那条支路被延伸，因而在每级比较后，幸存路径减少一半。因此，在接收到（K-1）个己知信息后，在整个网格图中就只有唯一的一条幸存路径保留下来，这就是译码所得的路径。也就是说，在己知接收到的序列的情况下，这条译码路径和发送序列是最相似的。

3. MATLAB仿真

在本次课程设计中，我们对整个通信过程进行了仿真，其过程如下： BPSK序列 信道 产生 编码 调制 AWGN 信道传输 ViterbiBPSK信息 输出 译码 解调

下面将分别对每一部分进行仿真

4

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

3.1卷积码编码仿真

在程序设计中，我们没有采用MATLAB自带的编码函数而是采用了自己的编码函数codec，其参数m为输入信息序列，g1，g2为两个输出端口的冲激响应序列。

function cod=codec(m,g1,g2)%g1,g2为两输出端口的冲激响应序列。 m1=conv(m,g1); m2=conv(m,g2); l=length(m1); for i=1:l;

cod([2*i-1])=rem(m1([i]),2); cod([2*i])=rem(m2([i]),2);

end

下为试运行编码结果，g1=[1 1 1],g2=[1 0 1]。

clear all g1=[1 1 1]; g2=[1 0 1]; msg=[1 1 0 1]; cod=codec(msg,g1,g2)

输出为：

cod =110101001011 符合预期结果。

3.2信道传输过程仿真

为了方便起见，我们采用了BPSK调制，为了了解整个通信系统，我们对BPSK的调制过程作了研究，用一个简短的程序对BPSK的全过程进行了观察。

function [bpsk_output]=bpsk_1(bpsk_input); %g=[1 0 1 1 1 0 0 1];基带信号 g=bpsk_input;

f=100; %载波频率 t=0:2*pi/99:2*pi; cp=[];sp=[];

mod=[];mod1=[];bit=[];

for n=1:length(g); if g(n)==0;

die=-ones(1,100); %Modulante se=ones(1,100); % end c=sin(f*t); cp=[cp die]; mod=[mod c]; bit=[bit se]; end

5

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

bpsk=cp.*mod;

subplot(2,1,1);plot(bit,'LineWidth',1.5);grid on; title('Binary Signal');

axis([0 100*length(g) -2.5 2.5]); bpsk_output=bpsk;

subplot(2,1,2);plot(bpsk,'LineWidth',1.5);grid on; title('BPSK modulation');

axis([0 100*length(g) -2.5 2.5]);

观察到如下波形

上图很好的符合了BPSK的调制特性。

信道中存在噪声，我们对高斯白噪声信道做了模拟，采用

noise=randn(1,length(bpsk))； sigma=sqrt(5)*10^(-w/20); recv=bpsk+3*sigma*noise;

来产生噪声并叠加到BPSK调制信号上，w为信噪比。

由于高斯白噪声的固有性质。若X= 0n t ∗φ(t)dt,则X为高斯随机变量，数学期望0，

N02

T

方差为σ2X=

T

0φ2(t)dt。则对于BPSK系统，接收信号为

T

T

T

y(Tb)= 0b[s t +nw t ]∗s1(t)dt= 0bs(t)∗s1(t)dt+ 0bnw∗s1(t)dt=±Eb+Z

其中s(t)= ±s1(t)。则知Z也为高斯随机变量。

由以上分析，为了简便起见，在主程序中可以直接将高斯白噪声叠加到输入[-1,1]序列上，即将调制过程简化为

bpsk=cod*2-1;

noise=randn(1,length(bpsk));

6

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

recv=bpsk+3*sigma*noise;

其中cod为编码序列。

解调时由于涉及软硬判决的信道信息问题，故留在Viterbi译码时一并讨论。

3.3译码过程仿真

维特比译码可分为软判决和硬判决两种，下面我们对其分别进行仿真研究。

3.3.1 硬判决译码

硬判决译码即比较汉明距作为累计度量，故直接将接收序列转为[0,1]序列即可，此时的判决代码为：

for j=1:length(recv); if recv(j)>0; recv(j)=1; else

recv(j)=0; end end

相应的译码程序为：

function dec=Hviterbi(cod) leth=length(cod)/2;

nextstate=[0 2;0 2;1 3;1 3];%储存下一个到达状态信息 outputs=[0 3;3 0;1 2;2 1];%储存状态转移时的输出信息。 suropt=zeros(4,leth-2);%使用2维数组保存幸存路径信息。 suroptem=suropt;

current=[0 0 0 0];%保存到达当前路径的累计汉明距值

current(1)=HMdist([0 0 0 0],cod([1 2 3 4]));%初始化前两步的码重度量信息 current(2)=HMdist([1 1 0 1],cod([1 2 3 4]));% current(3)=HMdist([0 0 1 1],cod([1 2 3 4]));% current(4)=HMdist([1 1 1 0],cod([1 2 3 4]));% suropt(1,[1 2])=[0 0];%前两步的输出信息 suropt(2,[1 2])=[1 0]; suropt(3,[1 2])=[0 1]; suropt(4,[1 2])=[1 1]; currenttem=[0 0 0 0];

ls=[1 0;2 0;3 0;4 0;1 1;2 1;3 1;4 1]; %保存能到达当前状态的前两个状态 for i=3:(leth-2); %开始译码

for x=1:4;

c=HMdist(deci2bin(outputs(ls(x*2-1,1),ls(x*2-1,2)+1),2),cod([i*2-1 i*2]));%有两个状态能转移到当前状态，计算其中一个转移过来时的输出与当前接收码的汉明距。

d=HMdist(deci2bin(outputs(ls(x*2,1),ls(x*2,2)+1),2),cod([i*2-1 i*2])); 计算另外一个状态转移过来时的输出与当前接收码的汉明距。

if (current(ls(x*2-1,1))+c)<(current(ls(x*2,1))+d);%比较找出累计码重较小的一个

currentem(x)=current(ls(x*2-1,1))+c;

7

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

suroptem(x,:)=suropt(ls(x*2-1,1),:); suroptem(x,i)=ls(x*2-1,2); else%=

currentem(x)=current(ls(x*2,1))+d;% suroptem(x,:)=suropt(ls(x*2,1),:); suroptem(x,i)=ls(x*2,2); end end

current=currentem; suropt=suroptem; end

a=[1 9999];

for s=1:4;%找出累计度量最小的路径，即最佳路径 if current(s)dec=suropt(a(1),:);%输出。

其中HMdist为比较两个序列的汉明距，deci2bin为将十进制数转换为二进制。 下图为在0~30dB信噪比下，硬判决维特比译码与无信道编码直接判决的对比。

由图中可见，在信噪比非常低的情况下，卷积码编码译码后的误码率反而比较高，因为误比特太多导致接受信息几乎无效，信噪比稍高后卷积码编码译码的误码率就大大降低。

8

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

程序流程图如下：

程序初始化，将路径信息和当前累计度量值清零，统计码长 对前两步进行汉明距度量并保存路径信息 判断是否解码完成 是 否 计算能到达当前状态的两条路径 将两条路径a,b输出与输入比较汉明距，并累加到路径度量 否 a路径累计度量是否小于b 是 选a为当前点的最优到达路径。 选b为当前点的最优到达路径。 选取度量值最小路径输出

9

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

3.3.2 软判决译码

由于硬判决在判决过程中将接收直接判决为[1 0],而损失了相应的近似信息，导致了接收信息在某种意义上的浪费。软判决将接收电平量化，尽可能保存相应的信息，故可以提高译码性能。 理想情况下软判决根据信道特性而选取量化方式，而对于实际产品如手机来说，信道特性是可变的，未知的，故在仿真中采取了均匀量化，相应的接收和量化程序如下。

msg=(randn(1,96)>0); sigma=sqrt(5)*10^(-w/20); cod=codec(msg,g1,g2); bpsk=cod*2-1;

noise=randn(1,length(bpsk)); recv=bpsk+3*sigma*noise; de=Sviterbi(recv,3)

即不判决而将直接将接收电平序列输入Viterbi软判决译码器。

量化程序如下：

function o=quantization(s,l); o=0; if s>1; o=2^l; elseif s<-1; o=-2^l; else c=abs(s); for i=1:2^l;

if c>1/2^(l+1)+(i-1)/2^l; o=o+1; end end if s<0; o=-o; end end

o=(o+2^l)/2^(l+1);

l为量化比特数，s为输入电平序列。程序将输入电平均匀量化为0~1之间的一级数。一共2l级电平。

软判决代码如下：

function dec=Sviterbi(recv,l); leth=length(recv)/2; for i=1:length(recv);

qr(i)=quantization(recv(i),l); end

nextstate=[0 2;0 2;1 3;1 3]; outputs=[0 3;3 0;1 2;2 1]; suropt=zeros(4,leth-2); suroptem=suropt;

10

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

current=[0 0 0 0];

current(1)=EUdist([0 0 0 0],qr([1 2 3 4])); current(2)=EUdist([1 1 0 1],qr([1 2 3 4])); current(3)=EUdist([0 0 1 1],qr([1 2 3 4])); current(4)=EUdist([1 1 1 0],qr([1 2 3 4])); suropt(1,[1 2])=[0 0]; suropt(2,[1 2])=[1 0]; suropt(3,[1 2])=[0 1]; suropt(4,[1 2])=[1 1]; currenttem=[0 0 0 0];

ls=[1 0 0;2 0 0;3 0 0;4 0 0;1 0 1;2 0 1;3 0 1;4 0 1]; for i=3:(leth-2); for x=1:4;

c=EUdist(deci2bin(outputs(ls(x*2-1,1),ls(x*2-1,3)+1),2),qr([i*2-1 i*2]));

d=EUdist(deci2bin(outputs(ls(x*2,1),ls(x*2,3)+1),2),qr([i*2-1 i*2]));

if (current(ls(x*2-1,1))+c)<(current(ls(x*2,1))+d); currentem(x)=current(ls(x*2-1,1))+c; suroptem(x,:)=suropt(ls(x*2-1,1),:); suroptem(x,i)=ls(x*2-1,3); else

currentem(x)=current(ls(x*2,1))+d; suroptem(x,:)=suropt(ls(x*2,1),:); suroptem(x,i)=ls(x*2,3); end end

current=currentem; suropt=suroptem; end

a=[1 9999]; for s=1:4

if current(s)dec=suropt(a(1),:);

其与硬判决的差异仅在于其比较欧氏距而非汉明距，并以欧氏距度量累计作为选择路径的标准。其理论上比硬判决有着3dB左右的优势。EUdist为欧氏距比较函数。 其与硬判决的差错率比较如下图

11

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

由图可见，软判决的确比硬判决有着一定的优势，在信噪比低的情况下尤其如此。而在现代通信环境下，软判决译码将比硬判决译码有着相当大的优势。

其总过程程序为

function mainprocedureS function mainprocedureS g1=[1 0 1]; g2=[1 1 1];

for w=1:20%选择1~20dB信噪比 q=0;

for i=1:200;%循环200次

msg=(randn(1,96)>0);%每次96比特信息 sigma=sqrt(5)*10^(-w/20); cod=codec(msg,g1,g2); bpsk=cod*2-1;

noise=randn(1,length(bpsk)); recv=bpsk+3*sigma*noise; %叠加噪声 de=Sviterbi(recv,3);%直接将电平送入译码器 q=q+HMdist(msg,de);%计算总误码数 end q=q+0

c(w)=10*log10(q/19200);%将误码率转为dB形式 end

for w=1:20;% a=0;

12

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

for i=1:200;

sigma=sqrt(5)*10^(-w/20); msg=(randn(1,96)>0); cod=codec(msg,g1,g2); bpsk=cod*2-1;

noise=randn(1,length(bpsk)); recv=bpsk+3*sigma*noise;

for j=1:length(recv);%先判决 再将[0 1]序列送入硬判决译码器 if recv(j)>0; recv(j)=1; else

recv(j)=0; end end

de=Hviterbi(recv);% a=a+HMdist(msg,de);% end a=a+0

e(w)=10*log10(a/19200); end

plot(c); %画图 grid on; hold on; plot(e);

4. 实验心得及总结

本次实验可以算是一波三折，一开始选择的课程设计其实并不是这个课题而是A律和μ律编码，但是由于组员较多，老师建议我们换一个课题，因为A律和μ律主要还是代码劳动，不需要太多的分析和逻辑。后来我们选择了Viterbi译码这个当时觉得很难的课题，彭老师专门为我们小组讲了一次Viterbi译码的基本原理，但是由于当时还没有接触到分组码卷积码等内容,听的其实不是很明白，后来自己在课余时间将信道编码那一章完全读过一次后才初步了解了卷积码和Viterbi译码的基本原理。但是由于各种各样的原因一直没有开始着手做，直到快到期末了才开始收集资料开始做。而这时我的队友已经初步的有了一个译码程序的雏形，本来打算就在队友的程序上修改一下用算了。但是觉得费力读懂别人的程序还不如自己做一个，所以最后自己写了整个程序，从开始对MATLAB基本上完全不会到后来能写出整个程序，途中查阅了很多书籍，参考了很多例子，总之还是付出了很多。而这次课程设计中获得是我觉得这一切都是值得的：我对信道编码整个一章有了更过更深刻的理解；我对信道特性对信息传输的影响有了直观的认识；我对信道编码的作用和受信噪比影响的特性有了更好的认识；我对卷积码编码的过程更加熟悉；我对Viterbi译码中的每一步都记忆深刻。这次课程设计还使我对MATLAB的编程有了更好的理解，充分体会到了MATLAB强大的数学计算能力。

最后我想我需要感谢彭岳星老师的详细讲解并提供了详细的资料；感谢我的队友所做的一切工作，才能使得这次课程设计能够顺利完成。

13

北京邮电大学信息与通信工程学院07114班<通信原理>课程设计之卷积码编码及Viterbi译码070422吴希龙

5. 参考文献

MATLAB 7.0 从入门到精通 人民邮电出版社

通信原理(第3版) 北京邮电大学出版社 周炯槃 庞沁华 续大我 通信原理 科学出版社 黄载禄 殷蔚华 A Viterbi Algorithm with Soft-Decision Outputs and its Applications Joachim Hagenauer Peter Hoeher

14