Turbo码在卷积码的基础上引入了一个交织序列，在生成的信息码字中引入了一定的随机性，其编码结构如下图：
对于信息向量
$$u=[u_1, u_2,...,u_{t^{'}},...,u_k]$$
分别经过两个卷积码编码器得到输出，其中$p^{2}$是第二个编码器输入经过了交织器乱序后的信息向量
$$p^{(1)}=[p^{(1)}1, p^{(1)}_2,...,p^{(1)}{t^{'}},...,p^{(1)}_k]$$
$$p^{(2)}=[p^{(2)}1, p^{(2)}_2,...,p^{(2)}{t^{'}},...,p^{(2)}_k]$$
对于1/3码率的Turbo码，其码字向量即为：
$$c=[u_1p^{(1)}_1p^{(2)}_1 u_2p^{(1)}_2p^{(2)}_2...]$$
显然，对于这种编码方式，其码率为1/3，在Turbo码中，引入了打孔矩阵，为
$$\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}$$
通过在第偶数个信息符号时去除$p^{(2)}$，在第奇数个信息符号时去除$p^{(1)}$，得到了打孔后的码字向量
$$c=[u_1p^{(1)}_1 u_2p^{(2)}_2...]$$
则此时的码率增加至$\frac{1}{2}$
编码过程如下(使用$(1,5/7)_8$卷积码)：
void encoder()
{
int s1 = 0, s2 = 0; //两个寄存器状态
int feedback = 0;
//第一个编码器
for (int i = 0; i < message_length - state_num; i++)
{
feedback = (s1 + s2) % 2;
p1[i] = (message[i] + feedback + s2) % 2;
s2 = s1;
s1 = (feedback + message[i]) % 2;
}
if (s1 == 0 && s2 == 0) //判断归零的输入
{
message[message_length - 2] = 0;
message[message_length - 1] = 0;
}
else if (s1 == 0 && s2 == 1)
{
message[message_length - 2] = 1;
message[message_length - 1] = 0;
}
else if (s1 == 1 && s2 == 0)
{
message[message_length - 2] = 1;
message[message_length - 1] = 1;
}
else if (s1 == 1 && s2 == 1)
{
message[message_length - 2] = 0;
message[message_length - 1] = 1;
}
for (int i = message_length - state_num; i < message_length; i++)
{
feedback = (s1 + s2) % 2;
p1[i] = (message[i] + feedback + s2) % 2;
s2 = s1;
s1 = (feedback + message[i]) % 2;
}
//第二个编码器
s1 = 0; s2 = 0;
feedback = 0;
for (int i = 0; i < message_length - state_num; i++)
{
feedback = (s1 + s2) % 2;
p2[i] = (message_mix[i] + feedback + s2) % 2;
s2 = s1;
s1 = (feedback + message_mix[i]) % 2;
}
if (s1 == 0 && s2 == 0) //判断归零的输入
{
message_mix[message_length - 2] = 0;
message_mix[message_length - 1] = 0;
}
else if (s1 == 0 && s2 == 1)
{
message_mix[message_length - 2] = 1;
message_mix[message_length - 1] = 0;
}
else if (s1 == 1 && s2 == 0)
{
message_mix[message_length - 2] = 1;
message_mix[message_length - 1] = 1;
}
else if (s1 == 1 && s2 == 1)
{
message_mix[message_length - 2] = 0;
message_mix[message_length - 1] = 1;
}
for (int i = message_length - state_num; i < message_length; i++)
{
feedback = (s1 + s2) % 2;
p2[i] = (message_mix[i] + feedback + s2) % 2;
s2 = s1;
s1 = (feedback + message_mix[i]) % 2;
}
//将上述结果结合得到最终的打孔后编码输出
for (int i = 0; i < message_length; i++)
{
codeword[i * 2] = message[i];
if (i % 2 == puncturing)
codeword[i * 2 + 1] = p1[i];
else codeword[i * 2 + 1] = p2[i]; //若当前为第偶数个信息符号，保留第一个编码器的输出，否则保留第二个
}
// convolution encoder, the input is message[] and the output is codeword[]
}
在译码时，流程如下：
在第一次迭代时，取$P_a(u_{t^{'}} = 0) = P_a(u_{t^{'}} = 1) = 0.5$，对于被打孔点，取信道观察概率为$P_{ch}(p^{(1)} = 0) = P_{ch}(p^{(1)} = 1) = 0.5$，那么在进行了一次BCJR译码后，可以得到后验概率$P_p(u_{t^{'}})$，和作为下一次迭代的先验概率
$$P_e(u_{t^{'}})=P(A|B)=\frac{P(AB)}{P(B)}$$
其中$P(AB)$为后验概率，$P(B)$为这一次迭代的先验概率，将这个概率作为下一次迭代的先验概率，最后可以得到一个趋近信道容量的编码。其仿真结果如下：
完整代码：
/***************************************************
Channel Coding Course Work: conolutional codes
This program template has given the message generator, BPSK modulation, AWGN channel model and BPSK demodulation,
you should first determine the encoder structure, then define the message and codeword length, generate the state table, write the convolutional encoder and decoder.
If you have any question, please contact me via e-mail: wuchy28@mail2.sysu.edu.cn
***************************************************/
#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <algorithm>
#include <random> // std::default_random_engine
#include <chrono> // std::chrono::system_clock
//设计时考虑了以下两个长度都为偶数
//否则decoder里面的aa等矩阵维度会出现小数
#define puncturing 0 // 若为0，则偶数位置保留第一个编码器的输出，若为1，偶数位保留第二个编码器的输出
#define message_length 1000 // the length of message
#define codeword_length 2000 // the length of codeword
#define INF 0x7fffff
#define iteration 10//迭代次数
float code_rate = (float)message_length / (float)codeword_length;
// channel coefficient
#define pi 3.1415926
long double N0, sgm;
unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
int state_table[10][10]; // state table, the size should be defined yourself
int state_num; // the number of the state of encoder structure
int message[message_length], codeword[codeword_length]; // message and codeword
int re_codeword[codeword_length]; // the received codeword
int de_message[message_length]; // the decoding message
int mix[message_length], reverse_mix[message_length]; //随机的乱序矩阵
int message_mix[message_length]; //打乱后的信息
int p1[message_length], p2[message_length]; //两个卷积码编码器的输出
long double pa1[message_length][2], pa2[message_length][2], pa1t[message_length][2], pa2t[message_length][2]; //两个BCJR译码器的输出(乱序后),一个临时变量
long double fix1[codeword_length][2], fix2[codeword_length][2];
long double tx_symbol[codeword_length][2]; // the transmitted symbols
long double rx_symbol[codeword_length][2]; // the received symbols
long double Pa[message_length][2];
void statetable();
void encoder();
void modulation();
void demodulation();
void channel();
void decoder(long double Pa[message_length][2], long double rx_symbol[codeword_length][2], int bcjr_num);
void fix();
int main()
{
int i;
float SNR, start, finish;
long int bit_error, seq, seq_num;
long double BER;
long double progress;
// generate state table
statetable();
for (int i = 0; i < message_length - state_num; i++)
mix[i] = i;
std::shuffle(mix, mix + message_length - state_num, std::default_random_engine(seed));
for (int i = 0; i < message_length - state_num; i++)
for (int j = 0; j < message_length - state_num; j++)
if (mix[i] == j)
reverse_mix[j] = i;
// random seed
srand((int)time(0));
for (int i = 0; i < message_length; i++)
for (int j = 0; j < 2; j++)
Pa[i][j] = 0.5;
// input the SNR and frame number
printf("\nEnter start SNR: ");
scanf("%f", &start);
printf("\nEnter finish SNR: ");
scanf("%f", &finish);
printf("\nPlease input the number of message: ");
scanf("%d", &seq_num);
for (SNR = start; SNR <= finish; SNR++)
{
// channel noise
N0 = (1.0 / code_rate) / pow(10.0, (float)(SNR) / 10.0);
//原始的信息每个比特能量为1（长度为k），经过编码后的信息长度变为n
//分母对应的SNR是以dBw为单位的
sgm = sqrt(N0 / 2); //噪声功率谱密度？
bit_error = 0;
for (seq = 1; seq <= seq_num; seq++)
{
// generate binary message randomly
/****************
Pay attention that message is appended by 0 whose number is equal to the state of encoder structure.
****************/
for (i = 0; i < message_length - state_num; i++)
message[i] = rand() % 2; //随机产生01
for (i = 0; i < message_length - state_num; i++)
message_mix[i] = message[mix[i]]; //进行乱序
for (i = message_length - state_num; i < message_length; i++)
message[i] = message_mix[i] = 0;
//对于每一个SNR，产生的序列都不同
// convolutional encoder
encoder();
// BPSK modulation
modulation();
// AWGN channel
channel();
// BPSK demodulation, it's needed in hard-decision Viterbi decoder
// demodulation();
fix(); //用于将接受序列转化为两个独立的正常顺序序列
// convolutional decoder
decoder(Pa, fix1, 1);
decoder(Pa, fix2, 2);
for (int i = 1; i <= iteration; i++)
{
for (int j = 0; j < message_length; j++)
for (int k = 0; k < 2; k++)
{
pa1t[j][k] = pa1[j][k];
pa2t[j][k] = pa2[j][k];
}
decoder(pa2t, fix1, 1);
decoder(pa1t, fix2, 2);
}
// calculate the number of bit error
for (i = 0; i < message_length - state_num; i++)
{
if (message_mix[i] != de_message[i])
bit_error++;
}
progress = (long double)(seq * 100) / (long double)seq_num;
// calculate the BER
BER = (long double)bit_error / (long double)(message_length * seq);
// print the intermediate result
printf("Progress=%2.1f, SNR=%2.1f, Bit Errors=%2.1d, BER=%E\r", progress, SNR, bit_error, BER);
}
// calculate the BER
BER = (long double)bit_error / (long double)(message_length * seq_num);
// print the final result
printf("Progress=%2.1f, SNR=%2.1f, Bit Errors=%2.1d, BER=%E\n", progress, SNR, bit_error, BER);
}
system("pause");
}
void statetable()
{
state_num = 2;
}
void encoder()
{
int s1 = 0, s2 = 0; //两个寄存器状态
int feedback = 0;
//第一个编码器
for (int i = 0; i < message_length - state_num; i++)
{
feedback = (s1 + s2) % 2;
p1[i] = (message[i] + feedback + s2) % 2;
s2 = s1;
s1 = (feedback + message[i]) % 2;
}
if (s1 == 0 && s2 == 0) //判断归零的输入
{
message[message_length - 2] = 0;
message[message_length - 1] = 0;
}
else if (s1 == 0 && s2 == 1)
{
message[message_length - 2] = 1;
message[message_length - 1] = 0;
}
else if (s1 == 1 && s2 == 0)
{
message[message_length - 2] = 1;
message[message_length - 1] = 1;
}
else if (s1 == 1 && s2 == 1)
{
message[message_length - 2] = 0;
message[message_length - 1] = 1;
}
for (int i = message_length - state_num; i < message_length; i++)
{
feedback = (s1 + s2) % 2;
p1[i] = (message[i] + feedback + s2) % 2;
s2 = s1;
s1 = (feedback + message[i]) % 2;
}
//第二个编码器
s1 = 0; s2 = 0;
feedback = 0;
for (int i = 0; i < message_length - state_num; i++)
{
feedback = (s1 + s2) % 2;
p2[i] = (message_mix[i] + feedback + s2) % 2;
s2 = s1;
s1 = (feedback + message_mix[i]) % 2;
}
if (s1 == 0 && s2 == 0) //判断归零的输入
{
message_mix[message_length - 2] = 0;
message_mix[message_length - 1] = 0;
}
else if (s1 == 0 && s2 == 1)
{
message_mix[message_length - 2] = 1;
message_mix[message_length - 1] = 0;
}
else if (s1 == 1 && s2 == 0)
{
message_mix[message_length - 2] = 1;
message_mix[message_length - 1] = 1;
}
else if (s1 == 1 && s2 == 1)
{
message_mix[message_length - 2] = 0;
message_mix[message_length - 1] = 1;
}
for (int i = message_length - state_num; i < message_length; i++)
{
feedback = (s1 + s2) % 2;
p2[i] = (message_mix[i] + feedback + s2) % 2;
s2 = s1;
s1 = (feedback + message_mix[i]) % 2;
}
//将上述结果结合得到最终的打孔后编码输出
for (int i = 0; i < message_length; i++)
{
codeword[i * 2] = message[i];
if (i % 2 == puncturing)
codeword[i * 2 + 1] = p1[i];
else codeword[i * 2 + 1] = p2[i]; //若当前为第偶数个信息符号，保留第一个编码器的输出，否则保留第二个
}
// convolution encoder, the input is message[] and the output is codeword[]
}
void modulation()
{
// BPSK modulation
int i;
// 0 is mapped to (1,0) and 1 is mapped tp (-1,0)
// tx_symbol[i][0]是实轴上的符号
for (i = 0; i < codeword_length; i++)
{
tx_symbol[i][0] = -1 * (2 * codeword[i] - 1);
tx_symbol[i][1] = 0;
}
}
void channel()
{
// AWGN channel
int i, j;
long double u, r, g;
for (i = 0; i < codeword_length; i++)
{
for (j = 0; j < 2; j++)
{
u = (float)rand() / (float)RAND_MAX;
if (u == 1.0)
u = 0.999999;
r = sgm * sqrt(2.0 * log(1.0 / (1.0 - u)));
u = (float)rand() / (float)RAND_MAX;
if (u == 1.0)
u = 0.999999;
g = (float)r * cos(2 * pi * u);
rx_symbol[i][j] = tx_symbol[i][j] + g;
}
}
}
void demodulation()
{
int i;
long double d1, d2;
for (i = 0; i < codeword_length; i++)
{
//解调，通过到星座点的距离来判断
d1 = (rx_symbol[i][0] - 1) * (rx_symbol[i][0] - 1) + rx_symbol[i][1] * rx_symbol[i][1];
d2 = (rx_symbol[i][0] + 1) * (rx_symbol[i][0] + 1) + rx_symbol[i][1] * rx_symbol[i][1];
if (d1 < d2)
re_codeword[i] = 0;
else
re_codeword[i] = 1;
}
}
void fix()
{
for (int i = 0; i < message_length; i++)
{
for (int j = 0; j < 2; j++)
{
fix1[i * 2][j] = rx_symbol[i * 2][j];
if (i < message_length - state_num)
fix2[i * 2][j] = rx_symbol[mix[i] * 2][j];
else fix2[i * 2][j] = rx_symbol[i * 2][j];
}
if (i % 2 == puncturing)
{
for (int j = 0; j < 2; j++)
{
fix1[i * 2 + 1][j] = rx_symbol[i * 2 + 1][j];
fix2[i * 2 + 1][j] = 0;
}
}
else
{
for (int j = 0; j < 2; j++)
{
fix2[i * 2 + 1][j] = rx_symbol[i * 2 + 1][j];
fix1[i * 2 + 1][j] = 0;
}
}
}
}
// BCJR decode
// re_symbol[codeword_length][2] long double type
void decoder(long double Pa[message_length][2], long double rx_symbol[codeword_length][2], int bcjr_num)
{
long double sum = 0;
long double Pch[codeword_length][2]; //用于存储Pch，0维度存到0的概率，1维度存到1的概率
long double d1, d2;
long double aa[codeword_length / 2], ac[codeword_length / 2], ba[codeword_length / 2], bc[codeword_length / 2];
long double cb[codeword_length / 2], cd[codeword_length / 2], db[codeword_length / 2], dd[codeword_length / 2]; //转移概率
//先算出转移概率，在算出A
long double A[codeword_length / 2 + 1][4]; //起始状态概率
A[0][0] = 1;
A[0][1] = 0;
A[0][2] = 0;
A[0][3] = 0;
for (int i = 0; i < codeword_length; i = i + 2)
{
d1 = (rx_symbol[i][0] - 1) * (rx_symbol[i][0] - 1) + rx_symbol[i][1] * rx_symbol[i][1];
Pch[i][0] = (1 / sqrt(pi * N0)) * exp(-d1 / N0);
d2 = (rx_symbol[i][0] + 1) * (rx_symbol[i][0] + 1) + rx_symbol[i][1] * rx_symbol[i][1];
Pch[i][1] = (1 / sqrt(pi * N0)) * exp(-d2 / N0);
long double sum = Pch[i][0] + Pch[i][1];
if (sum == 0) sum = 0.000001;
sum = 1 / sum;
Pch[i][0] *= sum;
Pch[i][1] *= sum;
d1 = (rx_symbol[i + 1][0] - 1) * (rx_symbol[i + 1][0] - 1) + rx_symbol[i + 1][1] * rx_symbol[i + 1][1];
Pch[i + 1][0] = (1 / sqrt(pi * N0)) * exp(-d1 / N0);
d2 = (rx_symbol[i + 1][0] + 1) * (rx_symbol[i + 1][0] + 1) + rx_symbol[i + 1][1] * rx_symbol[i + 1][1];
Pch[i + 1][1] = (1 / sqrt(pi * N0)) * exp(-d2 / N0);
if (rx_symbol[i + 1][0] == rx_symbol[i + 1][1] && rx_symbol[i + 1][0] == 0) //若当前点为打孔点，信道观察概率为0.5
Pch[i + 1][0] = Pch[i + 1][1] = 0.5;
sum = Pch[i + 1][0] + Pch[i + 1][1];
if (sum == 0) sum = 0.000001;
sum = 1 / sum;
Pch[i + 1][0] *= sum;
Pch[i + 1][1] *= sum;
bc[i / 2] = aa[i / 2] = Pa[i / 2][0] * Pch[i][0] * Pch[i + 1][0];
ba[i / 2] = ac[i / 2] = Pa[i / 2][1] * Pch[i][1] * Pch[i + 1][1];
cb[i / 2] = dd[i / 2] = Pa[i / 2][1] * Pch[i][1] * Pch[i + 1][0];
cd[i / 2] = db[i / 2] = Pa[i / 2][0] * Pch[i][0] * Pch[i + 1][1];
int k = i / 2;
if (k == message_length - 1)
break;
//前向迭代更新A节点
A[k + 1][0] = A[k][0] * aa[k] + A[k][1] * ba[k];
A[k + 1][1] = A[k][2] * cb[k] + A[k][3] * db[k];
A[k + 1][2] = A[k][1] * bc[k] + A[k][0] * ac[k];
A[k + 1][3] = A[k][2] * cd[k] + A[k][3] * dd[k];
//归一化
if (k == message_length - 2)
{
sum = A[k + 1][0] + A[k + 1][1];
if (sum == 0) sum = 0.000001;
A[k + 1][0] = A[k + 1][0] / sum;
A[k + 1][1] = A[k + 1][1] / sum;
A[k + 1][2] = 0;
A[k + 1][3] = 0;
A[k + 2][0] = 1;
A[k + 2][1] = A[k + 2][2] = A[k + 2][3] = 0;
}
else
{
sum = A[k + 1][0] + A[k + 1][1] + A[k + 1][2] + A[k + 1][3];
if (sum == 0) sum = 0.000001;
A[k + 1][0] = A[k + 1][0] / sum;
A[k + 1][1] = A[k + 1][1] / sum;
A[k + 1][2] = A[k + 1][2] / sum;
A[k + 1][3] = A[k + 1][3] / sum;
}
}
long double B[codeword_length / 2 + 1][4]; //终止状态概率
B[codeword_length / 2][0] = 1;
B[codeword_length / 2][1] = 0;
B[codeword_length / 2][2] = 0;
B[codeword_length / 2][3] = 0;
for (int i = message_length; i > 0; i--)
{
B[i - 1][0] = B[i][0] * aa[i - 1] + B[i][2] * ac[i - 1];
B[i - 1][1] = B[i][0] * ba[i - 1] + B[i][2] * bc[i - 1];
B[i - 1][2] = B[i][1] * cb[i - 1] + B[i][3] * cd[i - 1];
B[i - 1][3] = B[i][1] * db[i - 1] + B[i][3] * dd[i - 1];
sum = B[i - 1][0] + B[i - 1][1] + B[i - 1][2] + B[i - 1][3];
if (sum == 0) sum = 0.000001;
B[i - 1][0] = B[i - 1][0] / sum;
B[i - 1][1] = B[i - 1][1] / sum;
B[i - 1][2] = B[i - 1][2] / sum;
B[i - 1][3] = B[i - 1][3] / sum;
if (i == 2)
{
sum = B[i - 1][0] + B[i - 1][2];
if (sum == 0) sum = 0.000001;
B[i - 1][0] = B[i - 1][0] / sum;
B[i - 1][2] = B[i - 1][2] / sum;
B[i - 1][1] = B[i - 1][3] = 0;
B[0][0] = 1;
B[0][1] = B[0][2] = B[0][3] = 0;
break;
}
else
{
sum = B[i - 1][0] + B[i - 1][1] + B[i - 1][2] + B[i - 1][3];
if (sum == 0) sum = 0.000001;
B[i - 1][0] = B[i - 1][0] / sum;
B[i - 1][1] = B[i - 1][1] / sum;
B[i - 1][2] = B[i - 1][2] / sum;
B[i - 1][3] = B[i - 1][3] / sum;
}
}
long double p0, p1, pe0, pe1;
//计算后验概率
for (int i = 0; i < message_length - state_num; i++)
{
p0 = A[i][0] * aa[i] * B[i + 1][0] + A[i][1] * bc[i] * B[i + 1][2] + A[i][2] * cd[i] * B[i + 1][3] + A[i][3] * db[i] * B[i + 1][1];
p1 = A[i][0] * ac[i] * B[i + 1][2] + A[i][1] * ba[i] * B[i + 1][0] + A[i][2] * cb[i] * B[i + 1][1] + A[i][3] * dd[i] * B[i + 1][3];
long double sum = p0 + p1;
p0 /= sum;
p1 /= sum;
if (Pa[i][0] <= 0.000001) Pa[i][0] = 0.000001;
if (Pa[i][1] <= 0.000001) Pa[i][1] = 0.000001;
pe0 = p0 / Pa[i][0];
pe1 = p1 / Pa[i][1];
if (bcjr_num == 1)
{
pa1[i][0] = pe0 / (pe0 + pe1); //得到下一次迭代的先验概率,这是一个条件概率P(p0|pa[i][0])，为已知当前迭代的先验信息时，输出为0的概率
pa1[i][1] = 1 - pa1[i][0];
}
else if (bcjr_num == 2)
{
pa2[i][0] = pe0 / (pe0 + pe1);
pa2[i][1] = 1 - pa2[i][0];
}
if (p0 > p1) de_message[i] = 0;
else de_message[i] = 1;
pa1[message_length - 1][0] = pa1[message_length - 2][0] = pa1[message_length - 1][1] = pa1[message_length - 2][1] = 0.5;
pa2[message_length - 1][0] = pa2[message_length - 2][0] = pa2[message_length - 1][1] = pa2[message_length - 2][1] = 0.5;
}
long double temp[message_length][2];
for (int i = 0; i < message_length - state_num; i++) //对两个BCJR译码器的输出进行乱序或反乱序
for (int j = 0; j < 2; j++)
{
if (bcjr_num == 1)
temp[i][j] = pa1[i][j];
else if (bcjr_num == 2)
temp[i][j] = pa2[i][j];
}
for (int i = 0; i < message_length - state_num; i++)
for (int j = 0; j < 2; j++)
{
if (bcjr_num == 1)
pa1[i][j] = temp[mix[i]][j];
else if (bcjr_num == 2)
pa2[i][j] = temp[reverse_mix[i]][j];
}
return;
}
Comments NOTHING