吴雨霏博士版本的Turbo编译码仿真较为经典,以下就原码进行解读。

一、仿真代码架构

turbo_sys_demo.m是程序的主体框架,Turbo编译码均在此程序展开。程序开始需要用户需要如下几个参数:

1)译码算法:选择使用0:Log-MAP,1:SOVA

    这是让用户选择使用何种Turbo译码算法。若Tubo为多次迭代译码,则选用Log-MAP算法;否则选SOVA。默认选择0;

2)选择帧长帧长=信息长度+卷尾;

    根据自己需要选择帧长,程序默认选择400(396+4);

3)输入Turbo码生成多项式g默认为[111 101], 即 [7 5];

有两点需要注意:

1) 生成多项式为八进制表示;

2)若反馈多项式为7,即RSC分量码表示,则生成多项式表示形式为:[7 5;7]

 4)选择Turbo码是否打孔:0 – punctured;1—unpunctured

LTE Turbo未打孔码默认码率为1/3,打孔后默认码率为1/2;

  • 选择每帧迭代次数: Turbo译码一般迭代5 ~ 7次,程序默认选择5次;
  • 选择终止程序的帧错误次数:程序默认选择15帧;
  • 选择系统信噪比Eb/N0程序默认为2.0;

1为输入显示图:

代码参数输入仿真运行图

二、编码程序流程梳理

参数说明:

errs – 比特错误计数

nferr – 帧错误计数

nfame – 帧计数

1)产生数据源

源码中使用sort函数产生Turbo编码所使用的交织器,仿真时仍倾向采用LTE自带的交织器,代价是限定码长(因为每一个交织器都与特定的码长对应);

2)Turbo编码

 output = function encoderm( x, g, alpha, puncture )

  • 交织器映射采用 alpha数组;
  • puncture = 1,表示不需要Turbo码打孔,默认码率为1/3;
  • puncture = 0,表示对Turbo码打孔,输出码率为1/2;

2.1  Turbo编码流程梳理

先给出编码部分总体流程图

(不管是自己设计程序还是读别人程序,强烈建议梳理流程图这一步骤必不可少。因为认识事物的客观规律也是从宏观到微观,从整体到局部。如果一开始就陷入细节,对完成目标就很困难了):

编码整理流程图 


根据生成多项式矩阵g,得到以下参数:

  • 寄存器阶数 m = 矩阵g 列数 – 1;
  • 总编码信息长度 L_total = 信息长度 + m(加入尾比特处理,迫使寄存器状态最终归0)

程序调用 rsc_encode(g,input,1)函数:

1)首先,生成系统信息比特 d_k:    

% generate the codeword
for i = 1:L_total
   if terminated<0 | (terminated>0 & i<=L_info)
      d_k = x(1,i);
   elseif terminated>0 & i>L_info
      % terminate the trellis
      d_k = rem( g(1,2:K)*state', 2 );
   end

   a_k = rem( g(1,:)*[d_k state]', 2 );
   [output_bits, state] = encode_bit(g, a_k, state);

     若归零标志 terminated >0(代码里terminated = 1)并且编码序号 i  < 信息长度 L_info

     d_k = x(1,i);

     否则,d_k = rem( g(1,2:K) * state', 2 );

2)a_k = rem( g(1,:) * [d_k state]', 2 );

3)调用函数 encode_bit,输出output为总信息长度的2倍,即800bit:

% the rate is 1/n
% k is the constraint length
% m is the amount of memory
[n,k] = size(g);
m = k-1;

% determine the next output bit
for i=1:n
   output(i) = g(i,1)*input;
   for j = 2:k
      output(i) = xor(output(i),g(i,j)*state(j-1));
   end
end

state = [input, state(1:m-1)];

该函数详见流程图解析。如果对代码还不是很理解,那就说明对卷积码编码原理理解还不到位。这里推荐大家再回顾一下卷积码编码原理:

卷积码编码原理https://blog.csdn.net/snowman898/article/details/124148068


output1 = rsc_encode(g,input,1);
 
% make a matrix with first row corresponing to info sequence
% second row corresponsing to RSC #1's check bits.
% third row corresponsing to RSC #2's check bits.
 
y(1,:) = output1(1:2:2*L_total);
y(2,:) = output1(2:2:2*L_total);


y 的第一行为系统比特,即信息比特;

y的第二行为校验比特1;

同样的方法,  在encoder.m文件中再次调用rsc_encode(g,input,-1)函数。和上次不同的是,归零标志 terminated = -1,即校验码2并不归零。

这时,交织器开始起作用,输入rsc_encode里面也不再是原始数据,而是经过交织器交织后、长度变为L_total的数据:

% interleave input to second encoder
for i = 1:L_total
   input1(1,i) = y(1,alpha(i)); 
end

output2 = rsc_encode(g, input1(1,1:L_total), -1 );

% third row corresponsing to RSC #2's check bits.

y(3,:) = output2(2:2:2*L_total);

2.2  Turbo码打孔 (速率匹配)

(1)不打孔,即 puncture = 1

     rate = 1 / ( 2 + puncture ) ,则不打孔时默认速率为1/3;

    最终输出 en_output 按照 [ sys_data ; parity_data1;  parity_data2 ] 列模式方式输出,即          en_ouput 共有3行,有 L_total 列

(2)打孔,puncture = 0

     按照上述公式计算,打孔后 Turbo码率提高至 1/2;

   for i=1:L_total
       en_output(1,n*(i-1)+1) = y(1,i);
       if rem(i,2)
      % odd check bits from RSC1
          en_output(1,n*i) = y(2,i);
       else
      % even check bits from RSC2
          en_output(1,n*i) = y(3,i);
       end 
    end  

即:输出的奇数项为 系统信息, 偶数项 为 RSC1 校验位 和 RSC2 校验位。

最终,对en_output进行极性变换输出:

% antipodal modulation: +1/-1
en_output = 2 * en_output - ones(size(en_output));

三、 Turbo译码流程梳理

对下列参数进行初始化:

nferr , errs = 0,随后根据信噪比产生信号和噪声的叠加 recv signal:

 r = en_output+sigma*randn(1,L_total*(2+puncture)); % received bits

接下来,我们仍给出译码流程图:

译码整理流程图  

3.1 logmapo 函数介绍

 L_all = logmapo(rec_s,g,L_a,ind_dec)

输入参数:

1)rec_s:scaled received bits 缩放接收比特

     rec_s = 0.5 * L_c * yk = ( 2 * a * rate * Eb/N0 ) * yk

2)g:生产码多项式

3)L_a:先验 L 值

4)ind_dec: index of decoder. Either 1 or 2. Encoder 1 is assumed to be terminated, while                                encoder 2 is open.

该函数是turbo译码的核心,因此,先给出函数流程图:

 logmapo流程图 

3.2  Trellis函数介绍

trellis函数主要是根据码生成多项式,推出以下参数:

  • 后向输出  next_output
  • 后向状态  next_state
  • 前向输出  last_output
  • 前向状态  last_state

 1)首先,确定各个初始参数:

% g = [ 1 1 1; 1 0 1];

 [n,K] = size(g);             % n = 2, K = 3
m = K - 1;                      % m = 2
max_state = 2^m;         % max_state = 4

for state=1:max_state

    state_vector = bin_state( state-1, m );  % 将整数向量state转换成 m-bit 向量形式

end

state_vector 四个状态为:[ 0, 0 ] , [ 0, 1 ] , [ 1, 0 ] , [ 1, 1 ] 

2)当输入为0时:d_k = 0

a_k = rem( g(1,:)*[0 state_vector]', 2 );

[out_0, state_0] = encode_bit(g, a_k, state_vector);

可以认为根据 g, a_k, state_vector 完成一次 Turbo编码,得到 next_out, next_state;

encode_bit函数在编码流程中已详细梳理,此次不再赘述。最终,我们给出trellis网格图:

表1  LTE Turbo Trellis 

当前状态

state_vector

输入d_k=0

输入d_k=1

next_vector

(bit形式)

[ 0, 0 ]

[ 0, 0 ]

[ 1, 0 ]

[ 0, 1 ]

[ 1, 0 ]

[ 0, 0 ]

[ 1, 0 ]

[ 1, 1 ]

[ 0, 1 ]

[ 1, 1 ]

[ 0, 1 ]

[ 1, 1 ]

next_out

(bit形式)

[ 0, 0 ]

[ 0, 0 ]

[ 1, 1 ]

[ 0, 1 ]

[ 0, 0 ]

[ 1, 1 ]

[ 1, 0 ]

[ 0, 1 ]

[ 1, 0 ]

[ 1, 1 ]

[ 0, 1 ]

[ 1, 0 ]

next_out

(双极性形式)

[ 0, 0 ]

[ -1  -1  1  1 ]

[ 0, 1 ]

[ -1  -1  1  1 ]

[ 1, 0 ]

[ -1  1  1  -1 ]

[ 1, 1 ]

[ -1  1  1  -1 ]

next_state

(整数形式)

[ 0, 0 ]

[ 1  3]

[ 0, 1 ]

[ 3  1]

[ 1, 0 ]

[ 4  2]

[ 1, 1 ]

[ 2  4]

同时,根据该表,由当前 [d_k, next_out, next_state] 推导出前一时刻的 [ last_out, last_state]:

% find out which two previous states can come to present state
last_state = zeros(max_state,2);
for bit=0:1
   for state=1:max_state
      last_state(next_state(state,bit+1), bit+1)=state;
      last_out(next_state(state, bit+1), bit*2+1:bit*2+2) ...
         = next_out(state, bit*2+1:bit*2+2);
   end 
end

表2 LTE Turbo Trellis 

当前状态

state_vector

输入d_k=0

输入d_k=1

last_vector

(bit形式)

[ 0, 0 ]

[ 0, 0 ]

[ 0, 1 ]

[ 0, 1 ]

[ 1, 1 ]

[ 1, 0 ]

[ 1, 0 ]

[ 0, 1 ]

[ 0, 0 ]

[ 1, 1 ]

[ 1, 0 ]

[ 1, 1 ]

last_out

(bit形式)

[ 0, 0 ]

[ 0, 0 ]

[ 1, 1 ]

[ 0, 1 ]

[ 0, 1 ]

[ 1, 0 ]

[ 1, 0 ]

[ 0, 0 ]

[ 1, 1 ]

[ 1, 1 ]

[ 0, 1 ]

[ 1, 0 ]

last_out

(双极性形式)

[ 0, 0 ]

[ -1  -1  1  1 ]

[ 0, 1 ]

[ -1  1  1  -1 ]

[ 1, 0 ]

[ -1  -1  1  1 ]

[ 1, 1 ]

[ -1  1  1  -1 ]

last_state

(整数形式)

[ 0, 0 ]

[ 1  2]

[ 0, 1 ]

[ 4  3]

[ 1, 0 ]

[ 2  1]

[ 1, 1 ]

[ 3  4]

p.s. 这里提供另一角度去计算 last_out:

根据当前状态和输入 [ state,d_k ] 查表 last_state, 再由 [last_stae d_k] 查找 next_output;

3.3 前向度量计算

前向度量 α 计算由双循环组成:

内循环:以所有状态为循环:state = 1 : 2^m

外循环:以所有信息数量为循环:k = 2 : L_total + 1

先计算输入 d_k = 0时:

% Log_MAP算法计算

gamma( last_state(state2,1) ) = ( -rec_s(2*k-3) + rec_s(2*k-2) * last_out(state2,2) )
                                                   - log(1+exp(L_a(k-1)));

再计算输入 d_k = 1时: 

 gamma(last_state(state2,2)) = (rec_s(2*k-3)+rec_s(2*k-2)*last_out(state2,4))
                                                   + L_a(k-1) - log(1+exp(L_a(k-1)));

log_map算法目前用的不多,由于计算量较大,主流均采用 max-log-map计算,这里只给出计算公式,详细算法不再推导;

p.s. 未完,见后续(2)

       一定要看哈,文末有彩蛋(* ̄︶ ̄)

吴雨霏博士论文集及MATLAB原版程序https://download.csdn.net/download/qq_36756847/12036978?utm_medium=distribute.pc_relevant_download.none-task-download-2~default~BlogCommendFromBaidu~Rate-15-12036978-download-4229931.dl_default&depth_1-utm_source=distribute.pc_relevant_download.none-task-download-2~default~BlogCommendFromBaidu~Rate-15-12036978-download-4229931.dl_default&dest=https%3A%2F%2Fdownload.csdn.net%2Fdownload%2Fqq_36756847%2F12036978&spm=1003.2020.3001.6616.17

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐