【信道编码/Channel Coding】卷积码和Viterbi译码及其MATLAB实现

简介：

---

这是本专栏信道编码/Channel Coding的最后一站，想对信道编码有一个系统性的认识可以看本专栏的 信道编码的整体框架 一文。而在本篇文章中，将介绍卷积码的基本原理和Viterbi译码的过程，以及其MATLAB实现。为什么是最后一站呢，明明还有Turbo码，LDPC码和Polar码，这是因为这是我本科的总结，那些码由于实现起来需要很多细节，有很多研究价值，我希望在我的研究生生涯能够完善。

目录

简介：

一、卷积码的原理

1.1 卷积码的符号表示

1.2 运作原理

二、 举例：（2,1,2）卷积编码器

2.1 状态转移表

2.2 状态转移图

2.3 网格图

三、Viterbi译码

3.1 译码过程（以上述（2，1，2）卷积码为例）（假设用硬判决）

3.2 硬判决和软判决

四、MATLAB实现

---

一、卷积码的原理

---

1.1 卷积码的符号表示

和分组码稍有不同，卷积码的符号表示有三个维度，（n，k，L）中，n和k和分组码不同，n是每次输出的码字长度，k是每次输入的信息比特长度，而L就是记忆深度，卷积编码器的图表示如下：

1.2 运作原理

• 每过一个时间单位，比如从 t=0 到 t=1， k个比特从输入端输入k位的移位寄存器，而每一个移位寄存器里的k个比特输入到这个寄存器之后的一个寄存器。总共有L个k位移位寄存器。
• 紧接着上面的，有L个k位寄存器意味着有k*L个比特。根据设计的不同，每一个比特对n个模二加运算器有独特的影响。
• 模二加运算器的输出结果作为输出，所以一次输出n个比特。
• 移位寄存器初始状态都是0

二、 举例：（2,1,2）卷积编码器

---

 可以看到，（2,1,2）卷积编码器是一个，每一时刻只有一个输入，但是会输出两位（两个模二加运算器）的卷积编码器，里面的寄存器都是1位的，共有两个寄存器（两个记忆深度）。而模二加运算器的数学表达如下：

 假如输入序列为 x=10101，则其规律会是这样：

2.1 状态转移表

我们会发现，寄存器的状态就只有00,01,10,11四种可能，而输入只有0/1两种可能，所以输出一定也是伪随机的，也就是说输入输出的关系只有2x4=8种可能，我们将其列表：

输入	当前状态	下个状态	输出1	输出2
0	00	00	0	0
0	01	00	1	1
0	10	01	1	0
0	11	01	0	1
1	00	10	1	1
1	01	10	0	0
1	10	11	0	1
1	11	11	1	0

这就是状态转移表，通过查表可以非常简单的看出不同寄存器状态对应的输入输出之间的关系。这种表结构常用于编程。

2.2 状态转移图

我们只需要稍微改一下，将表转化成图就能得到状态转移图，如下：

 箭头表示从某一状态到每一状态，0/10表示在这个状态下，输入0得到输出10.

2.3 网格图

 网格图也非常直观简单，每行是一种寄存器的状态，所以有四行分别是00,01,10,11，而箭头仍然表示从什么状态到什么状态，0/00仍然表示在这个状态下输入0则输出11。但是不同的是，网格图能体现时刻的变化导致状态的转移路线的变化。

三、Viterbi译码

---

Viterbi译码实际上并不是非常复杂，其基于最大似然估计思想（ML）。我们先通过一个简单的例子迅速的理解一下Viterbi译码的流程。

3.1 译码过程（以上述（2，1，2）卷积码为例）（假设用硬判决）

• 将接收到的码两两分组（毕竟编码的时候一个输入对应两个输出）
• 利用和编码器一样的网格图，此时t=0，从状态00出发，画出输入为0和输入为1的两条路线,。并将输入为0的输出 和 输入为1的输出 分别和 t=0 要翻译的码进行异或（计算汉明距离），得到的值称为路径度量（PM，Path Metric），我们将其作为像路径的重量一样累积起来。注意，这时候应该有两个路径了、
• t=1，t=2，...的操作和上一步一样，从每个路径在上一时刻的状态出发，比较输入为0输入为1的输出和此时刻要翻译的码，计算汉明距离，再计算路径度量（累积的汉明距离）。所以在t=n时刻，理应有  个路径
• 剪枝（Pruning）：由于路径过于多，所以在路径达到某个数量的时候，译码器仅保留路径度量（累积汉明距离）最小的两个路径，然后继续译码。直至译码结束。

举例：

假设接收端接收到信号，并将其量化判决为011100（体现了硬判决，下面会讲）

（1）第一步是把信息两两分组：01 11 00，作为t=0，t=1，t=2时刻的输入

（2）第二步，译码，按照上面所述的流程，t=0时，收到的码为00；此时译码器状态为00，输入为0时输出为00，和收到的码的汉明距离是0，故PM=0；输入为1时输出为11，和收到的码汉明距离是2，故PM=2.

（3）继续译码，从上一步的两条路径出发，会衍生出四条路径

（4）进行剪枝（看系统的设计，不一定是这个时候剪枝），只保留最小PM的两个路径，继续译码

（5）回溯，从PM最短的路径开始回溯，译码为110

【注意】通常接收端会在一帧的信息比特发完末尾加入0 0，这样接收端的译码器就会回归00状态

3.2 硬判决和软判决

大家都注意到了，上面栗子的小标题提示了这是硬判决，那么什么是硬判决什么是软判决呢？

• 硬判决（Hard Decision）是指在接收端，信号被接收后，先进行种种滤波，均衡等操作，然后直接被判决器判为0或1，根据得到的比特串进行如上面所述的Viterbi译码。注意，这时候我们的PM（路径度量）是以汉明距离为度量。由它的二元特性判断，它比较适合BPSK等二元系统。
• 软判决（Soft Decision）是指在接收端，信号被接收后，进行完种种滤波，均衡等操作，在采样量化后，先不进行判决（也就是先不把多电平状态的量化信号变成0101的数字信号），而是直接经过软判决Viterbi译码器，这时候译码器会以欧氏距离为PM（路径度量），判断某一电平状态距离欧氏距离最近的符号。软判决通常会比硬判决性能更好，因为软判决减少了因量化器判决导致的量化噪声。

四、MATLAB实现

---

%%
clear all;clc;close all;
% Generation of Bits stream
nob=100; % The amount of bit
bit_stream=[randi([0 1],nob,1);[0;0]];

codeword=zeros(2*(nob+2),1);
D1=0; % Initial state of registers
D2=0;
for i_b=1:1:nob+2
    y1=mod(bit_stream(i_b)+D1+D2,2);
    y2=mod(bit_stream(i_b)+D2,2);
    D2=D1;
    D1=bit_stream(i_b);
    codeword(2*i_b-1)=y1;
    codeword(2*i_b)=y2;
end

%%
% Viterbi Decoding
received_codeword=codeword;
state_table=[0 0;0 1;1 0;1 1];  %find(ismember(state_table,[1 1],'row')==1)
Transfer_table={[0 0],[0 0],[1 1],[1 0]; % Output with input 0/1;Next State
    [1 1],[0 0],[0 0],[1 0];
    [1 0],[0 1],[0 1],[1 1];
    [0 1],[0 1],[1 0],[1 1]};
D=[0 0]; % State of registers
branch=cell(8,4); % weight, estimated bit stream, last state
branch{1,1}=0;branch{1,2}=[];branch{1,3}=[0 0];branch{1,4}=1;
for i=2:1:8
    branch{i,1}=0;
    branch{i,4}=0;
end

for i_b=1:2:length(received_codeword)
    code_buffer=received_codeword(i_b:i_b+1,1).';
    a_b=find([branch{:,4}]==0); % Available Branches
    n_b=8-length(a_b);
    for i_s=1:1:n_b
        i_a_b1=i_s;
        i_a_b2=a_b(i_s);
        weight=branch{i_s,1};  % Accumulated weight of branch
        est_b_s=branch{i_s,2}; % Estimated Bit stream
        last_state=branch{i_s,3}; % Last state of this branch
        index=find(ismember(state_table,last_state,'row')==1); % Find the index of Last state
        
        % Look up table
        weight1=weight+sum(xor(code_buffer,Transfer_table{index,1})); % Hamming Distance
        weight2=weight+sum(xor(code_buffer,Transfer_table{index,3}));
        est_b_s1=[est_b_s,0];
        est_b_s2=[est_b_s,1];
        last_state1=Transfer_table{index,2};
        last_state2=Transfer_table{index,4};
        
        % Update the branch list
        branch{i_a_b1,1}=weight1;
        branch{i_a_b1,2}=est_b_s1;
        branch{i_a_b1,3}=last_state1;
        branch{i_a_b1,4}=1;
        branch{i_a_b2,1}=weight2;
        branch{i_a_b2,2}=est_b_s2;
        branch{i_a_b2,3}=last_state2;
        branch{i_a_b2,4}=1;
    end
    
    % Pruning
    if n_b==4
        weight_list=[branch{:,1}];
        [~,i_tmw]=sort(weight_list); % Get the index of the two minimum weight
        buf11=branch{i_tmw(1),1};
        buf12=branch{i_tmw(1),2};
        buf13=branch{i_tmw(1),3};
        buf21=branch{i_tmw(2),1};
        buf22=branch{i_tmw(2),2};
        buf23=branch{i_tmw(2),3};
        
        branch{1,1}=buf11;
        branch{1,2}=buf12;
        branch{1,3}=buf13;
        branch{1,4}=1;
        branch{2,1}=buf21;
        branch{2,2}=buf22;
        branch{2,3}=buf23;
        branch{2,4}=1;
        for i=3:1:8
            branch{i,4}=0;
        end
    end
end

%%
n_survival_branch=length(find([branch{:,4}]==1));
w_survival_branch=[branch{1:n_survival_branch,1}];
[~,i_weight]=sort(w_survival_branch);
i_minimum_weight=i_weight(1);
estimated_bit_stream=branch{i_minimum_weight,2};
BER=sum(abs(estimated_bit_stream-bit_stream'))/(nob+2)