BCJR_RSC.c

/***************************************************
Channel Coding Course Work: RSC (Recursive Systematic Convolutional) Code
Systematic Code with Trellis Termination
***************************************************/

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

#define message_length 65536 // 消息长度
#define codeword_length (message_length * 2) // 码字长度
float code_rate = (float)message_length / (float)codeword_length;

// Math constants
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

double N0, sgm;

int state_table[4][4]; // 状态表

// [修正] 这里 state_num 保持为4，但在逻辑上它对应尾比特长度
// 对于K=3 (Memory=2) 的码，理论上只需要2位，但按照要求保留4位
int tail_length = 4; 

int message[message_length];
int codeword[codeword_length];
int re_codeword[codeword_length];
int de_message[message_length];

double tx_symbol[codeword_length][2];
double rx_symbol[codeword_length][2];

// BCJR Variables
double alpha[message_length + 1][4];
double beta[message_length + 1][4];

void statetable();
void encoder();
void modulation();
void channel();
void bcjr_decoder();

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

    FILE *fp = NULL;
    fp = fopen("bcjr_rsc_soft.csv", "w");
    if (fp == NULL) {
        printf("Error: Could not open file for writing.\n");
        return 1;
    }
    fprintf(fp, "SNR,BER\n");

    // 生成状态表 (RSC结构)
    statetable();

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

    printf("\n--- Soft Decision BCJR (MAP) Simulation for RSC ---\n");
    printf("Enter start SNR: ");
    scanf("%f", &start);
    printf("Enter finish SNR: ");
    scanf("%f", &finish);
    printf("Please input the number of message frames: ");
    scanf("%ld", &seq_num);

    for (SNR = start; SNR <= finish; SNR += step)
    {
        // SNR = Es/N0 (BPSK), Es=1. Code Rate R=1/2.
        // Eb/N0 = Es/N0 * (1/R) -> N0 = 1/(R * 10^(SNR/10))
        N0 = (1.0 / code_rate) / pow(10.0, (double)(SNR) / 10.0);
        sgm = sqrt(N0 / 2.0);

        bit_error = 0;

        for (seq = 1; seq <= seq_num; seq++)
        {
            // 1. 生成随机信息 (前 65532 位)
            for (i = 0; i < message_length - tail_length; i++)
            {
                message[i] = rand() % 2;
            }
            
            // [修正] 后4位先置0，但在encoder中会根据状态回溯进行修改
            for (i = message_length - tail_length; i < message_length; i++)
            {
                message[i] = 0; 
            }

            // 2. 编码 (含网格终止计算)
            encoder();

            // 3. BPSK 调制
            modulation();

            // 4. AWGN 信道
            channel();

            // 5. BCJR 解码
            bcjr_decoder();

            // 6. 计算误码 (注意：尾比特通常不计入误码，但这里全长计算)
            // 严格来说应该只统计有效信息位，但为保持代码结构简单，统计全长
            for (i = 0; i < message_length; i++)
            {
                if (message[i] != de_message[i])
                    bit_error++;
            }

            progress = (double)(seq * 100) / (double)seq_num;
            BER = (double)bit_error / (double)(message_length * seq);
            printf("Progress=%2.1f, SNR=%2.1f, Bit Errors=%ld, BER=%E\r", progress, SNR, bit_error, BER);
        }

        BER = (double)bit_error / (double)(message_length * seq_num);
        fprintf(fp, "%f,%E\n", SNR, BER);
        printf("SNR=%2.1f, Bit Errors=%ld, BER=%E\n", SNR, bit_error, BER);
    }

    fclose(fp);
    printf("\nResult saved to bcjr_rsc_soft.csv\n");

    return 0;
}

// [修正] 生成 RSC (1, 5/7) 状态表
// 状态定义 s = (D1, D2) -> (s0, s1)
// 输入 u
// 反馈输入 a = u + s0 + s1 (mod 2)
// 下一状态 (a, s0)
// 输出1 (系统): u
// 输出2 (校验): a + s1 (即 1+D^2)
void statetable()
{
    int s, u;
    int s0, s1; // s0是高位(D1), s1是低位(D2)
    int feedback_input;
    int next_state;
    int out_sys, out_par, output;

    for (s = 0; s < 4; s++)
    {
        s1 = s & 1;       // D2
        s0 = (s >> 1) & 1;// D1

        for (u = 0; u <= 1; u++)
        {
            // 反馈多项式 7 (111): 1 + D + D^2
            // 这里的加法是模2加 (异或)
            feedback_input = u ^ s0 ^ s1;

            // 下一状态: 新位移入高位，原高位移至低位
            // next_state = (feedback_input, s0)
            next_state = (feedback_input << 1) | s0;

            // 输出1: 系统位 (Systematic)
            out_sys = u;

            // 输出2: 前馈多项式 5 (101): 1 + D^2
            // 物理上对应 tap 在输入端(a)和D2(s1)端
            out_par = feedback_input ^ s1;

            output = (out_sys << 1) | out_par;

            // 存入表
            // state_table[s][0/1] -> next_state
            // state_table[s][2/3] -> output
            // 这里我们按照 encoder/decoder 的习惯索引：
            // 列 0: 输入0的下一状态
            // 列 1: 输入0的输出
            // 列 2: 输入1的下一状态
            // 列 3: 输入1的输出
            state_table[s][u * 2] = next_state;
            state_table[s][u * 2 + 1] = output;
        }
    }
}

// [修正] 编码器：实现 RSC 编码及网格终止
void encoder()
{
    int i;
    int current_state = 0;
    int input_bit;
    int output_symbol;

    for (i = 0; i < message_length; i++)
    {
        // 区分数据阶段和终止阶段
        if (i < message_length - tail_length)
        {
            // 正常数据位
            input_bit = message[i];
        }
        else
        {
            // [关键修正] 网格终止 (Trellis Termination)
            // 此时我们要强制让寄存器回归0状态。
            // 在 RSC 中，意味着进入寄存器的反馈值 (feedback_input) 必须为 0。
            // feedback_input = u ^ s0 ^ s1 = 0
            // 因此，必须令 u = s0 ^ s1
            int s1 = current_state & 1;
            int s0 = (current_state >> 1) & 1;
            input_bit = s0 ^ s1;
            
            // 修改 message 数组记录真实的输入，以便计算误码率
            message[i] = input_bit;
        }

        int col_next = (input_bit == 0) ? 0 : 2;
        int col_out  = (input_bit == 0) ? 1 : 3;

        output_symbol = state_table[current_state][col_out];
        int next_state = state_table[current_state][col_next];

        codeword[2 * i] = (output_symbol >> 1) & 1;
        codeword[2 * i + 1] = output_symbol & 1;

        current_state = next_state;
    }
}

void modulation()
{
    int i;
    for (i = 0; i < codeword_length; i++)
    {
        // 0 -> +1, 1 -> -1
        tx_symbol[i][0] = -1.0 * (2.0 * codeword[i] - 1.0);
        tx_symbol[i][1] = 0.0;
    }
}

void channel()
{
    int i, j;
    double u1, u2, r, g;

    for (i = 0; i < codeword_length; i++)
    {
        for (j = 0; j < 2; j++)
        {
            // Box-Muller transform
            do {
                u1 = (double)rand() / RAND_MAX;
            } while (u1 <= 1e-9); // 避免 log(0)
            
            u2 = (double)rand() / RAND_MAX;
            
            r = sgm * sqrt(-2.0 * log(u1));
            g = r * cos(2.0 * M_PI * u2);

            rx_symbol[i][j] = tx_symbol[i][j] + g;
        }
    }
}

// BCJR Decoder 保持原有逻辑框架，但在计算 Gamma 时会自动利用新的 state_table
void bcjr_decoder()
{
    int t, state, input;
    
    // 1. 初始化
    for (t = 0; t <= message_length; t++) {
        for (state = 0; state < 4; state++) {
            alpha[t][state] = 0.0;
            beta[t][state] = 0.0;
        }
    }

    alpha[0][0] = 1.0;
    // 由于执行了网格终止，最后时刻必然回到状态 0
    beta[message_length][0] = 1.0; 

    // 2. 前向递归 (Forward Recursion)
    for (t = 0; t < message_length; t++) {
        double sum_alpha = 0.0;
        
        for (int cur_s = 0; cur_s < 4; cur_s++) {
            if (alpha[t][cur_s] < 1e-50) continue; // 忽略极小概率

            for (input = 0; input <= 1; input++) {
                int col_next = (input == 0) ? 0 : 2;
                int col_out  = (input == 0) ? 1 : 3;

                int next_s = state_table[cur_s][col_next];
                int output_sym = state_table[cur_s][col_out];

                int c0 = (output_sym >> 1) & 1;
                int c1 = output_sym & 1;

                double tx0 = (c0 == 0) ? 1.0 : -1.0;
                double tx1 = (c1 == 0) ? 1.0 : -1.0;

                double r0 = rx_symbol[2 * t][0];
                double r1 = rx_symbol[2 * t + 1][0];

                double dist_sq = (r0 - tx0)*(r0 - tx0) + (r1 - tx1)*(r1 - tx1);
                double gamma = exp(-dist_sq / (2.0 * sgm * sgm));

                alpha[t + 1][next_s] += alpha[t][cur_s] * gamma;
            }
        }

        // 归一化
        for (state = 0; state < 4; state++) sum_alpha += alpha[t + 1][state];
        if (sum_alpha > 0) {
            for (state = 0; state < 4; state++) alpha[t + 1][state] /= sum_alpha;
        }
    }

    // 3. 后向递归 (Backward Recursion)
    for (t = message_length - 1; t >= 0; t--) {
        double sum_beta = 0.0;

        for (int cur_s = 0; cur_s < 4; cur_s++) {
            for (input = 0; input <= 1; input++) {
                int col_next = (input == 0) ? 0 : 2;
                int col_out  = (input == 0) ? 1 : 3;

                int next_s = state_table[cur_s][col_next];
                int output_sym = state_table[cur_s][col_out];

                int c0 = (output_sym >> 1) & 1;
                int c1 = output_sym & 1;
                double tx0 = (c0 == 0) ? 1.0 : -1.0;
                double tx1 = (c1 == 0) ? 1.0 : -1.0;
                double r0 = rx_symbol[2 * t][0];
                double r1 = rx_symbol[2 * t + 1][0];
                double dist_sq = (r0 - tx0)*(r0 - tx0) + (r1 - tx1)*(r1 - tx1);
                double gamma = exp(-dist_sq / (2.0 * sgm * sgm));

                beta[t][cur_s] += beta[t + 1][next_s] * gamma;
            }
        }

        // 归一化
        for (state = 0; state < 4; state++) sum_beta += beta[t][state];
        if (sum_beta > 0) {
            for (state = 0; state < 4; state++) beta[t][state] /= sum_beta;
        }
    }

    // 4. LLR 计算与判决
    for (t = 0; t < message_length; t++) {
        double prob0 = 0.0;
        double prob1 = 0.0;

        for (int cur_s = 0; cur_s < 4; cur_s++) {
            // 这里只需判断 alpha 是否有值，避免无效计算
            if (alpha[t][cur_s] < 1e-50) continue; 

            for (input = 0; input <= 1; input++) {
                int col_next = (input == 0) ? 0 : 2;
                int col_out  = (input == 0) ? 1 : 3;
                int next_s = state_table[cur_s][col_next];
                int output_sym = state_table[cur_s][col_out];

                int c0 = (output_sym >> 1) & 1;
                int c1 = output_sym & 1;
                double tx0 = (c0 == 0) ? 1.0 : -1.0;
                double tx1 = (c1 == 0) ? 1.0 : -1.0;
                double r0 = rx_symbol[2 * t][0];
                double r1 = rx_symbol[2 * t + 1][0];
                double dist_sq = (r0 - tx0)*(r0 - tx0) + (r1 - tx1)*(r1 - tx1);
                double gamma = exp(-dist_sq / (2.0 * sgm * sgm));

                double metric = alpha[t][cur_s] * gamma * beta[t + 1][next_s];

                if (input == 0) prob0 += metric;
                else            prob1 += metric;
            }
        }

        if (prob1 > prob0) de_message[t] = 1;
        else               de_message[t] = 0;
    }
}