turbo.cpp

/***************************************************
Channel Coding Course Work: Turbo Codes (Final Working Version)
Turbo 码仿真程序 (C++)
包含功能：
1. RSC 编码器 (递归系统卷积码)
2. 随机交织器 (Random Interleaver)
3. BPSK 调制与 AWGN 信道
4. 迭代解码 (Iterative Decoding)
5. SISO 解码器 (Max-Log-MAP 算法)
***************************************************/

#define  _CRT_SECURE_NO_WARNINGS
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>

// --- 系统参数定义 ---
#define MSG_LEN 1024    // 信息比特长度 (每帧数据量)
#define STATE_NUM 4      // 状态数量 (RSC生成多项式为 1, 5/7，约束长度 K=3，状态数=2^(K-1)=4)
#define ITERATIONS 8     // Turbo 解码的迭代次数
#define TAIL_BITS 4      // 归零尾比特 (用于让编码器回到全0状态)

// 总消息长度 = 数据长度 + 尾比特
#define message_length (MSG_LEN + TAIL_BITS) 
// 码字长度 = 系统比特 + 校验比特1 + 校验比特2 (码率 R=1/3)
#define codeword_length (message_length * 3) 

float code_rate = 1.0f / 3.0f; // 码率

// --- 全局常量与变量 ---
#define pi 3.141592653589793
#define INF 1e9        // 定义无穷大，用于初始化对数概率

double N0, sgm; // N0: 噪声功率谱密度, sgm: 噪声标准差 (sigma)

// --- 网格图 (Trellis) 结构 ---
// 用于描述卷积码的状态转移
int next_state[STATE_NUM][2];    // [当前状态][输入比特] -> 下一状态
int output_parity[STATE_NUM][2]; // [当前状态][输入比特] -> 输出校验位

// --- 交织器映射表 ---
int alpha_interleaver[message_length]; // 存储交织后的索引位置

// --- 数据缓冲区 ---
int message[message_length];        // 原始消息比特 (u)
int codeword[codeword_length];      // 编码后的码字
int re_codeword[codeword_length];   // (未使用的保留变量)
int de_message[message_length];     // 解码后的硬判决结果

// --- 调制/信道缓冲区 ---
double tx_symbol[codeword_length][2]; // 发送符号 (BPSK映射后)
double rx_symbol[codeword_length][2]; // 接收符号 (经过信道加噪后)

// --- 解码器内部缓冲区 (接收到的LLR) ---
double rx_sys[message_length];   // 接收到的系统比特信息 (Systematic)
double rx_par1[message_length];  // 接收到的校验比特1信息 (Parity 1)
double rx_par2[message_length];  // 接收到的校验比特2信息 (Parity 2)

// --- 外部信息 (Extrinsic Information) ---
// 用于在两个分量解码器之间交换信息
double L_ext1[message_length]; // 解码器1 输出的外部信息
double L_ext2[message_length]; // 解码器2 输出的外部信息

// --- 函数声明 ---
void statetable(); // 初始化网格图状态表和交织器
void encoder();    // Turbo 编码主函数
void modulation(); // BPSK 调制
void channel();    // AWGN 信道模拟
void decoder();    // Turbo 解码主函数 (迭代控制)
void rsc_encoder(int* input_bits, int* parity_out, int len); // RSC 分量编码器
// SISO (软输入软输出) 解码核心函数
void siso_decoder(double* L_in_sys, double* L_in_par, double* L_in_apriori, double* L_out_extrinsic, int len, int terminate);
double max_star(double a, double b); // Max-Log-MAP 运算核心

int main()
{
    int i;
    float SNR, start, finish, step = 0.2;
    long int bit_error, seq, seq_num;
    double BER;
    double progress;

    // --- 新增：定义文件指针 ---
    FILE *fp = NULL;

    statetable();
    srand((unsigned int)time(0));

    // --- 新增：打开 CSV 文件并写入表头 ---
    fp = fopen("turbo_ber.csv", "w");
    if (fp == NULL) {
        printf("Error: Could not open file turbo_ber.csv for writing.\n");
        return 1;
    }
    fprintf(fp, "SNR,BER\n"); // 写入 CSV 表头

    // 初始化状态表和随机数种子
    statetable();
    srand((unsigned int)time(0));

    printf("--- Turbo Code Simulation (Final Version) ---\n");
    printf("Message Length: %d, Codeword Length: %d\n", message_length, codeword_length);
    
    // 用户输入仿真参数
    printf("\nEnter start SNR (dB) [suggest -2.0]: ");
    scanf("%f", &start);
    
    printf("\nEnter finish SNR (dB) [suggest 3.0]: ");
    scanf("%f", &finish);
    
    printf("\nPlease input number of frames [e.g., 50]: ");
    scanf("%ld", &seq_num);

    // --- SNR 循环 ---
    for (SNR = start; SNR <= finish; SNR += step)
    {
        // 计算噪声参数
        // Eb/N0 = SNR_dB, Es/N0 = Eb/N0 * Rate
        // N0 = 1 / (Rate * 10^(SNR/10))
        N0 = (1.0 / code_rate) / pow(10.0, (float)(SNR) / 10.0);
        sgm = sqrt(N0 / 2.0); // 高斯白噪声的标准差

        bit_error = 0; // 重置误码计数

        // --- 帧/序列循环 ---
        for (seq = 1; seq <= seq_num; seq++)
        {
            // 1. 生成随机消息 (最后 TAIL_BITS 位假装不设置)
            for (i = 0; i < message_length; i++)
                message[i] = rand() % 2;
            for (i = message_length - TAIL_BITS; i < message_length; i++)
                message[i] = 0; //先置尾比特为0

            // 2. 编码 (Turbo Encoder)
            encoder();

            // 3. 调制 (BPSK)
            modulation();

            // 4. 过信道 (Add Gaussian Noise)
            channel();

            // 5. 解码 (Iterative Decoder)
            decoder();

            // 6. 统计误码 (Error Counting)
            // 只统计有效信息位，不统计尾比特
            for (i = 0; i < message_length - TAIL_BITS; i++)
            {
                if (message[i] != de_message[i])
                    bit_error++;
            }

            // 打印进度
            progress = (double)(seq * 100) / (double)seq_num;
            if (seq % 10 == 0 || seq == seq_num) {
                BER = (double)bit_error / (double)((message_length - TAIL_BITS) * seq);
                printf("Progress=%3.0f%%, SNR=%4.1f, Errors=%6ld, BER=%E\r", progress, SNR, bit_error, BER);
            }
        }
        // 计算该 SNR 点下的最终 BER
        BER = (double)bit_error / (double)((message_length - TAIL_BITS) * seq_num);
        
        // --- 新增：将结果写入文件 ---
        fprintf(fp, "%f,%E\n", SNR, BER);
        
        printf("Progress=%3.0f%%, SNR=%4.1f, Errors=%6ld, BER=%E\n", progress, SNR, bit_error, BER);
    }
    
    printf("\nSimulation finished. Press Enter to exit.");
    getchar(); getchar();

    // --- 新增：关闭文件 ---
    fclose(fp);

    return 0;
}

// --- 初始化函数 ---
void statetable()
{
    int s, u;
    // 构建 RSC (1, 5/7) 的状态转移表
    // 生成多项式八进制表示: Feedforward=5 (101), Feedback=7 (111)
    for (s = 0; s < STATE_NUM; s++) {
        for (u = 0; u < 2; u++) {
            // 获取当前状态的二进制位 (s = m1 m0)
            int m1 = (s >> 1) & 1; 
            int m0 = s & 1;        
            
            // 反馈多项式逻辑 (Feedback=7 -> 1 + D + D^2)
            int a_k = (u + m1 + m0) % 2; // 输入 + 反馈
            
            // 前馈多项式逻辑 (Feedforward=5 -> 1 + D^2)
            // 输出 parity = a_k + m0 (对应D^2)
            int out = (a_k + m0) % 2;
            
            // 下一状态移位 (Next State = a_k m1)
            int next_s = (a_k << 1) | m1;
            
            next_state[s][u] = next_s;
            output_parity[s][u] = out;
        }
    }

    // 生成随机交织器 (Random Interleaver)
    // 算法：Knuth shuffle (洗牌算法)
    int i, j, temp;
    for (i = 0; i < message_length; i++) alpha_interleaver[i] = i; // 初始化为顺序
    for (i = message_length - 1; i > 0; i--) {
        j = rand() % (i + 1); // 随机选一个位置
        // 交换
        temp = alpha_interleaver[i];
        alpha_interleaver[i] = alpha_interleaver[j];
        alpha_interleaver[j] = temp;
    }
}

// --- 单个 RSC 编码器 ---
void rsc_encoder(int* input_bits, int* parity_out, int len)
{
    int s = 0; // 初始状态为 0
    int i;
    for (i = 0; i < len; i++) {
        int u = input_bits[i];               // 当前输入比特
        parity_out[i] = output_parity[s][u]; // 查表得校验位
        s = next_state[s][u];                // 状态更新
    }
}

// --- Turbo 编码器 (Parallel Concatenation) ---
// --- 修改后的 Turbo 编码器 ---
void encoder()
{
    int i, k;
    int parity1[message_length];
    int parity2[message_length];
    int interleaved_msg[message_length];
    
    // ==========================================
    // 新增：计算尾比特以终止编码器1 (Trellis Termination)
    // ==========================================
    int s = 0; // 模拟编码器1的状态
    int u, m1, m0, a_k;

    // 1. 先跑完前面的数据位，获取当前状态
    for (i = 0; i < message_length - TAIL_BITS; i++) {
        u = message[i];
        // 只有这里需要模拟状态跳转，逻辑同 statetable
        m1 = (s >> 1) & 1; 
        m0 = s & 1;       
        a_k = (u + m1 + m0) % 2; // 反馈
        s = (a_k << 1) | m1;     // 下一状态
    }

    // 2. 计算剩余的 4 个尾比特
    // 目标：使反馈 a_k = 0，即 (u + m1 + m0) % 2 = 0 => u = (m1 + m0) % 2
    for (i = message_length - TAIL_BITS; i < message_length; i++) {
        m1 = (s >> 1) & 1; 
        m0 = s & 1;
        
        // 计算能让反馈归零的输入比特 u
        // 对于反馈多项式 7 (111): d_k = u + m1 + m0. 要让 d_k=0，则 u = m1 + m0
        message[i] = (m1 + m0) % 2; 

        // 更新状态 (此时 a_k 必定为 0)
        u = message[i];
        a_k = (u + m1 + m0) % 2; 
        s = (a_k << 1) | m1;
    }
    // ==========================================

    // 执行交织：将消息（包含刚才计算好的尾比特）顺序打乱
    for(i=0; i<message_length; i++) interleaved_msg[i] = message[alpha_interleaver[i]];

    // 编码器1：输入经过归零处理的序列
    rsc_encoder(message, parity1, message_length);
    
    // 编码器2：输入交织后的序列 
    // 注意：交织后的序列通常无法让编码器2也归零，这是Turbo码的标准特性
    rsc_encoder(interleaved_msg, parity2, message_length);

    // 组装码字
    for (i = 0; i < message_length; i++) {
        codeword[3 * i] = message[i];
        codeword[3 * i + 1] = parity1[i];
        codeword[3 * i + 2] = parity2[i];
    }
}

// --- BPSK 调制 ---
void modulation()
{
    // 映射规则: 0 -> +1.0, 1 -> -1.0
    int i;
    for (i = 0; i < codeword_length; i++)
    {
        tx_symbol[i][0] = (codeword[i] == 0) ? 1.0 : -1.0;
        tx_symbol[i][1] = 0; // 虚部为0
    }
}

// --- AWGN 信道 ---
void channel()
{
    int i, j;
    double u, r, g;
    // Box-Muller 变换生成高斯噪声
    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; // 接收信号 = 发送 + 噪声
        }
    }
}

// --- Max-Log-MAP 运算核心 ---
// 用于近似计算 log(e^a + e^b) ≈ max(a, b)
// 真正的 Log-MAP 会加上修正项 log(1 + e^-|a-b|)
double max_star(double a, double b) {
    return (a > b) ? a : b;
}

// --- SISO 解码器 (Soft-Input Soft-Output) ---
// 实现 BCJR (MAP) 算法的对数域版本
// L_in_sys: 系统位 LLR
// L_in_par: 校验位 LLR
// L_in_apriori: 先验信息 (来自另一个解码器的外信息)
// L_out_extrinsic: 输出的外部信息
// len: 序列长度
// terminate: 是否已知终止状态 (1: 已知, 0: 未知)
void siso_decoder(double* L_in_sys, double* L_in_par, double* L_in_apriori, double* L_out_extrinsic, int len, int terminate)
{
    // 定义前向概率 Alpha, 后向概率 Beta, 转移概率 Gamma
    // 维度说明：[时间 t][状态 s]
    static double alpha[MSG_LEN + TAIL_BITS + 1][STATE_NUM];
    static double beta[MSG_LEN + TAIL_BITS + 1][STATE_NUM];
    static double gamma[MSG_LEN + TAIL_BITS][STATE_NUM][2]; // [时间][状态][输入比特0/1]

    int k, s, u;

    // 1. 初始化 Alpha (前向度量)
    // t=0时刻，编码器必定在状态0
    for (s = 0; s < STATE_NUM; s++) {
        if (s == 0) alpha[0][s] = 0.0;
        else alpha[0][s] = -INF;
    }

    // 2. 初始化 Beta (后向度量)
    // 已知 trellis 终止 (terminate=1)，则 t=len 时必定在状态0
    for (s = 0; s < STATE_NUM; s++) {
        if (terminate) beta[len][s] = (s == 0) ? 0.0 : -INF;
        else beta[len][s] = 0.0;
    }

    // 3. 计算 Gamma (分支度量 Branch Metric)
    // Gamma(s -> s', u) = 0.5 * [ u * L_a(u) + u * L_c * y_s + p * L_c * y_p ]
    for (k = 0; k < len; k++) {
        for (s = 0; s < STATE_NUM; s++) {
            for (u = 0; u < 2; u++) {
                int parity = output_parity[s][u];
                // BPSK 符号: 0 -> +1, 1 -> -1
                double sign_u = (u == 0) ? 1.0 : -1.0;
                double sign_p = (parity == 0) ? 1.0 : -1.0;
                
                // 计算分支度量值
                double metric = 0.5 * (sign_u * (L_in_sys[k] + L_in_apriori[k]) + sign_p * L_in_par[k]);
                gamma[k][s][u] = metric;
            }
        }
    }

    // 4. Alpha 递归 (Forward Recursion)
    // Alpha_k(s) = max* (Alpha_{k-1}(prev_s) + Gamma(prev_s -> s))
    for (k = 1; k <= len; k++) {
        for (s = 0; s < STATE_NUM; s++) {
            alpha[k][s] = -INF;
            int prev_s, input_bit;
            // 遍历所有可能的前一状态
            for(prev_s = 0; prev_s < STATE_NUM; prev_s++) {
                for(input_bit = 0; input_bit < 2; input_bit++) {
                    if (next_state[prev_s][input_bit] == s) {
                        alpha[k][s] = max_star(alpha[k][s], alpha[k-1][prev_s] + gamma[k-1][prev_s][input_bit]);
                    }
                }
            }
        }
    }

    // 5. Beta 递归 (Backward Recursion)
    // Beta_k(s) = max* (Beta_{k+1}(next_s) + Gamma(s -> next_s))
    for (k = len - 1; k >= 0; k--) {
        for (s = 0; s < STATE_NUM; s++) {
            beta[k][s] = -INF;
            for (u = 0; u < 2; u++) {
                int n_s = next_state[s][u]; // 下一状态
                beta[k][s] = max_star(beta[k][s], beta[k+1][n_s] + gamma[k][s][u]);
            }
        }
    }

    // 6. 软输出与外部信息计算 (Extrinsic Calculation)
    // LLR = log( P(u=0)/P(u=1) )
    for (k = 0; k < len; k++) {
        double L0 = -INF; // u=0 的联合概率对数
        double L1 = -INF; // u=1 的联合概率对数

        for (s = 0; s < STATE_NUM; s++) {
            // 对所有输入为0的分支进行 max* 累加
            int ns0 = next_state[s][0];
            L0 = max_star(L0, alpha[k][s] + gamma[k][s][0] + beta[k+1][ns0]);
            
            // 对所有输入为1的分支进行 max* 累加
            int ns1 = next_state[s][1];
            L1 = max_star(L1, alpha[k][s] + gamma[k][s][1] + beta[k+1][ns1]);
        }
        
        double L_total = L0 - L1; // 后验对数似然比 (A Posteriori LLR)
        
        // 减去输入的系统信息和先验信息，得到纯外部信息
        // L_ext = L_total - L_channel_sys - L_apriori
        double ext = L_total - L_in_sys[k] - L_in_apriori[k];
        
        // 限幅处理：防止数值溢出导致计算不稳定
        if (ext > 50.0) ext = 50.0;
        if (ext < -50.0) ext = -50.0;
        L_out_extrinsic[k] = ext;
    }
}

// --- Turbo 解码器主控 ---
void decoder()
{
    int i, iter;
    // 信道可靠性因子 Lc = 4 * Es / N0 = 2 / sigma^2
    double Lc = 2.0 / (sgm * sgm); 

    // 初始化输入 LLR
    for (i = 0; i < message_length; i++) {
        rx_sys[i]  = Lc * rx_symbol[3 * i][0];     // 系统位 LLR
        rx_par1[i] = Lc * rx_symbol[3 * i + 1][0]; // 校验位1 LLR
        rx_par2[i] = Lc * rx_symbol[3 * i + 2][0]; // 校验位2 LLR
        L_ext2[i] = 0.0; // 初始先验信息为 0 (等概)
    }

    // 开始迭代解码
    for (iter = 0; iter < ITERATIONS; iter++) {
        // --- 解码器 1 (自然顺序) ---
        // 输入：系统位，校验位1，来自DEC2的先验信息(L_ext2)
        // 输出：L_ext1 (传给DEC2)
        siso_decoder(rx_sys, rx_par1, L_ext2, L_ext1, message_length, 1);// 终止状态已知

        // --- 交织处理 ---
        // 为了让解码器2使用，需要将系统信息和外信息进行交织
        double rx_sys_int[message_length];
        double L_apriori2[message_length];
        for (i = 0; i < message_length; i++) {
            rx_sys_int[i] = rx_sys[alpha_interleaver[i]]; // 系统位交织
            L_apriori2[i] = L_ext1[alpha_interleaver[i]]; // 外信息1交织 -> 成为DEC2的先验
        }

        // --- 解码器 2 (交织顺序) ---
        // 输入：交织后的系统位，校验位2，来自DEC1的先验信息
        // 输出：L_ext2 (交织域)
        siso_decoder(rx_sys_int, rx_par2, L_apriori2, L_ext2, message_length, 0);

        // --- 解交织处理 ---
        // 解码器2输出的 L_ext2 是在交织顺序下的，需要映射回自然顺序传给解码器1
        double temp[message_length];
        for (i = 0; i < message_length; i++) temp[i] = L_ext2[i]; // 暂存交织域输出
        
        for (i = 0; i < message_length; i++) {
            // 解交织逻辑：alpha[i] 是原序列第i位在交织后的位置
            // 所以，交织序列的第 alpha[i] 位，对应原序列的第 i 位
            L_ext2[alpha_interleaver[i]] = temp[i];
        }
    }

    // --- 硬判决 (Hard Decision) ---
    // 最终判决基于 L_total = L_sys + L_ext1 + L_ext2
    for (i = 0; i < message_length; i++) {
        double final_LLR = rx_sys[i] + L_ext1[i] + L_ext2[i];
        if (final_LLR >= 0) de_message[i] = 0; // LLR > 0 判决为 0
        else de_message[i] = 1;                // LLR < 0 判决为 1
    }
}

void demodulation() {} // 占位函数，本程序中解调在 decoder 初始化部分完成