Viterbi_Sort_RSC.c

/***************************************************
Channel Coding Course Work: Soft Decision Viterbi (RSC 1, 5/7)
Description:
1. Soft Decision: Uses Euclidean distance metric.
2. Code: RSC (1, 5/7)_8.
3. Termination: Tail-biting logic preserved.
***************************************************/

#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) 
#define TERMINATION_BITS 4             

float code_rate = 0.5; 

#define pi 3.1415926
#define INF 1.0e9 // 软判决度量值可能较大，使用较大的浮点数作为无穷大

// --- 全局变量 ---
double N0, sgm;
int state_table[4][4]; 
int state_num = 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]; // 接收到的软信息 (浮点数)

// 路径记录
int traceback_path[message_length][4];

// --- 函数声明 ---
void statetable_RSC();
void encoder();
void modulation();
void channel();
// void demodulation(); // 软判决不需要此步骤
void decoder();

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

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

    statetable_RSC(); 
    srand((int)time(0));

    printf("\n--- Soft Decision Viterbi (RSC 1, 5/7) ---\n");
    printf("Enter start SNR: ");
    scanf("%f", &start);
    printf("Enter finish SNR: ");
    scanf("%f", &finish);
    printf("Please input the number of message: ");
    scanf("%ld", &seq_num);

    for (SNR = start; SNR <= finish; SNR += step)
    {
        N0 = (1.0 / code_rate) / pow(10.0, (float)(SNR) / 10.0);
        sgm = sqrt(N0 / 2);
        bit_error = 0;

        for (seq = 1; seq <= seq_num; seq++)
        {
            // 1. 生成数据 (含归零位预留)
            for (i = 0; i < message_length - TERMINATION_BITS; i++) {
                message[i] = rand() % 2;
            }
            for (; i < message_length; i++) {
                message[i] = 0; 
            }

            // 2. 编码 (RSC + Termination)
            encoder();

            // 3. 调制与信道
            modulation();
            channel();
            // demodulation(); // 软判决跳过此步，直接把 rx_symbol 送给 decoder

            // 4. 软判决译码
            decoder();

            // 5. 误码统计
            for (i = 0; i < message_length; i++) {
                if (message[i] != de_message[i]) bit_error++;
            }

            if (seq % 50 == 0 || seq == seq_num) {
                progress = (double)(seq * 100) / (double)seq_num;
                BER = (double)bit_error / (double)(message_length * seq);
                printf("SNR=%4.1f Progress=%5.1f%% BER=%E\r", SNR, progress, BER);
            }
        }

        BER = (double)bit_error / (double)(message_length * seq_num);
        fprintf(fp, "%f,%E\n", SNR, BER);
        printf("\nFile saved: SNR=%4.1f BER=%E\n", SNR, BER);
    }

    fclose(fp);
    return 0;
}

void statetable_RSC()
{
    // RSC (1, 5/7) 状态表
    // s0=D^1, s1=D^2
    int s, u, s0, s1, a_in, next_state, output;

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

        for (u = 0; u <= 1; u++) {
            a_in = u ^ s0 ^ s1;
            next_state = (a_in << 1) + s0;

            int out1 = u;          // Systematic
            int out2 = a_in ^ s1;  // Parity (1+D^2)

            output = (out1 << 1) + out2;

            state_table[s][u * 2]     = next_state; 
            state_table[s][u * 2 + 1] = output;     
        }
    }
}

void encoder()
{
    int i, input_bit, output_symbol;
    int current_state = 0;
    int data_len = message_length - TERMINATION_BITS;

    for (i = 0; i < message_length; i++)
    {
        // 归零逻辑 (Trellis Termination)
        if (i < data_len) {
            input_bit = message[i];
        } 
        else {
            int s1 = current_state & 1;
            int s0 = (current_state >> 1) & 1;
            input_bit = s0 ^ s1; // 强制归零输入
            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()
{
    for (int i = 0; i < codeword_length; i++) {
        // 0 -> +1.0, 1 -> -1.0
        tx_symbol[i][0] = -1.0 * (2.0 * codeword[i] - 1.0);
        tx_symbol[i][1] = 0;
    }
}

void channel()
{
    for (int i = 0; i < codeword_length; i++) {
        double u1 = (double)rand() / RAND_MAX; if (u1 >= 1.0) u1 = 0.999999;
        double u2 = (double)rand() / RAND_MAX; if (u2 >= 1.0) u2 = 0.999999;
        double g = sgm * sqrt(-2.0 * log(1.0 - u1)) * cos(2 * pi * u2);
        
        rx_symbol[i][0] = tx_symbol[i][0] + g;
        rx_symbol[i][1] = 0;
    }
}

void decoder()
{
    int t, state, input;
    // --- 修改点1: 路径度量改为 double ---
    double path_metric[4];
    double new_path_metric[4];

    // 初始化
    path_metric[0] = 0.0;
    for (state = 1; state < 4; state++) path_metric[state] = INF;

    // --- 前向递推 ---
    for (t = 0; t < message_length; t++)
    {
        // --- 修改点2: 读取软信息 (浮点数) ---
        double r0 = rx_symbol[2 * t][0];
        double r1 = rx_symbol[2 * t + 1][0];

        for (state = 0; state < 4; state++) new_path_metric[state] = INF;

        for (int curr_state = 0; curr_state < 4; curr_state++)
        {
            if (path_metric[curr_state] >= INF) continue;

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

                int next_state = state_table[curr_state][col_next];
                int output_sym = state_table[curr_state][col_out];

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

                // --- 修改点3: 计算欧氏距离平方 (Euclidean Distance Squared) ---
                // 期望符号: bit 0 -> +1.0, bit 1 -> -1.0
                double expected_0 = (c0 == 0) ? 1.0 : -1.0;
                double expected_1 = (c1 == 0) ? 1.0 : -1.0;

                // 分支度量 = (r0 - e0)^2 + (r1 - e1)^2
                // 展开后等价于最小化 (r - e)^2
                double dist0 = r0 - expected_0;
                double dist1 = r1 - expected_1;
                double branch_metric = (dist0 * dist0) + (dist1 * dist1);

                double total_metric = path_metric[curr_state] + branch_metric;

                // 更新路径 (找最小距离)
                if (total_metric < new_path_metric[next_state]) {
                    new_path_metric[next_state] = total_metric;
                    traceback_path[t][next_state] = curr_state;
                }
            }
        }
        for (state = 0; state < 4; state++) path_metric[state] = new_path_metric[state];
    }

    // --- 回溯 (不变) ---
    // 由于有 Trellis Termination，直接从状态 0 开始
    int curr_state = 0; 
    
    for (t = message_length - 1; t >= 0; t--)
    {
        int prev_state = traceback_path[t][curr_state];
        
        // 查表反推输入
        if (state_table[prev_state][0] == curr_state) {
            de_message[t] = 0;
        } 
        else if (state_table[prev_state][2] == curr_state) {
            de_message[t] = 1;
        }
        else {
            de_message[t] = 0; 
        }

        curr_state = prev_state;
    }
}