sccc_final_csv.c

/***************************************************
Channel Coding Course Work: SCCC (Final Version with CSV Output)
功能:
1. SCCC (串行级联卷积码) 仿真
2. 高性能优化 (OpenMP, Log-MAP)
3. 结果自动保存到 sccc_ber.csv 文件，方便绘图
***************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include <omp.h> // OpenMP 并行库

// --- 系统参数定义 ---
#define INFO_LEN 256        // 信息位长度 (短帧)
#define STATE_NUM 4         // 状态数
#define ITERATIONS 10       // 迭代次数
#define TAIL_BITS 2         // 尾比特

// SCCC 长度定义
#define LEN_OUTER_IN  (INFO_LEN + TAIL_BITS)
#define LEN_OUTER_OUT (LEN_OUTER_IN * 2)
#define LEN_INTER     LEN_OUTER_OUT
#define LEN_INNER_IN  (LEN_INTER + TAIL_BITS)
#define LEN_INNER_OUT (LEN_INNER_IN * 2)

float code_rate = 0.25f; // R = 1/4

// --- 常量 ---
#define PI 3.141592653589793f
#define INF 1e9f

// --- 全局变量 ---
int next_state[STATE_NUM][2];       
int output_parity[STATE_NUM][2];    
int sccc_interleaver[LEN_INTER];    

// --- 查找表优化 Log-MAP (提速关键) ---
// [FIX 1] 修改为整数常量，解决 variably modified 错误
#define LUT_PRECISION 100
#define LUT_MAX_DIFF 8 
#define LUT_SIZE (LUT_MAX_DIFF * LUT_PRECISION + 1)

float correction_table[LUT_SIZE];

// --- RNG ---
typedef struct { unsigned int x, y, z, w; } rng_state_t;
void rng_init(rng_state_t *rng, unsigned int seed) {
    rng->x = seed; rng->y = 362436069; rng->z = 521288629; rng->w = 88675123;
}
unsigned int rng_next(rng_state_t *rng) {
    unsigned int t = rng->x ^ (rng->x << 11);
    rng->x = rng->y; rng->y = rng->z; rng->z = rng->w;
    return rng->w = rng->w ^ (rng->w >> 19) ^ (t ^ (t >> 8));
}
float rng_float(rng_state_t *rng) {
    return (float)rng_next(rng) / (float)0xFFFFFFFF;
}

// --- [FIX 2] 将 max_star 和 init_table 移到上方，并删除原来的 max_star 原型 ---

void init_correction_table() {
    for (int i = 0; i < LUT_SIZE; i++) {
        float x = (float)i / (float)LUT_PRECISION;
        correction_table[i] = log1pf(expf(-x));
    }
}

// 查表版 Log-MAP (static inline 需要在使用前定义)
static inline float max_star(float a, float b) {
    float max_val = (a > b) ? a : b;
    float diff = fabsf(a - b);
    if (diff >= (float)LUT_MAX_DIFF) return max_val;
    // 转换为整数索引
    int index = (int)(diff * (float)LUT_PRECISION);
    return max_val + correction_table[index];
}

// --- 函数原型 ---
void statetable();
// void init_correction_table(); // 已在上方定义，不需要原型
void outer_encoder(int* input_bits, int* output_stream, int len);
// float max_star(float a, float b); // [FIX 2] 删除此行，避免 static/non-static 冲突
void siso_generic(float* L_in_sys, float* L_in_par, float* L_in_apriori_sys, float* L_in_apriori_par,
                  float* L_out_sys, float* L_out_par, int len, int terminate, int mode,
                  float (*alpha)[STATE_NUM], float (*beta)[STATE_NUM], float (*gamma)[STATE_NUM][2]);

// --- 主函数 ---
int main()
{
    float start_snr, finish_snr, current_snr;
    long int seq_num;
    FILE *fp = NULL; // 定义文件指针

    // 1. 初始化
    statetable();
    init_correction_table();

    // 2. 打开 CSV 文件准备写入
    fp = fopen("sccc_ber.csv", "w");
    if (fp == NULL) {
        printf("Error: Could not open file sccc_ber.csv for writing.\n");
        return 1;
    }
    // 写入 CSV 表头
    fprintf(fp, "SNR,BER\n");
    printf("Result will be saved to 'sccc_ber.csv' dynamically.\n");

    printf("\n--- SCCC Simulation (Rate 1/4) ---\n");
    printf("Info Length: %d, Transmitted Bits: %d\n", LEN_OUTER_IN, LEN_INNER_OUT);
    printf("Detected Processors: %d\n", omp_get_num_procs());
    
    printf("\nEnter start SNR (dB) [suggest 0.0]: ");
    if(scanf("%f", &start_snr));
    printf("\nEnter finish SNR (dB) [suggest 4.5]: ");
    if(scanf("%f", &finish_snr));
    printf("\nPlease input number of frames [e.g. 100000]: ");
    if(scanf("%ld", &seq_num));

    for (current_snr = start_snr; current_snr <= finish_snr; current_snr += 0.2f) 
    {
        double N0 = (1.0 / code_rate) / pow(10.0, current_snr / 10.0);
        float sgm = (float)sqrt(N0 / 2.0);
        float Lc = 2.0f / (sgm * sgm);

        long int total_bit_error = 0;
        double t_start = omp_get_wtime();

        // --- OpenMP 并行计算 ---
        #pragma omp parallel reduction(+:total_bit_error)
        {
            // 内存分配
            int *msg = (int*)calloc(LEN_OUTER_IN, sizeof(int));
            int *outer_coded = (int*)malloc(LEN_OUTER_OUT * sizeof(int));
            int *inner_in = (int*)calloc(LEN_INNER_IN, sizeof(int)); 
            
            float *rx_inner_sys = (float*)malloc(LEN_INNER_IN * sizeof(float));
            float *rx_inner_par = (float*)malloc(LEN_INNER_IN * sizeof(float));
            
            float *L_inner_ext = (float*)malloc(LEN_INNER_IN * sizeof(float)); 
            float *L_outer_in_sys = (float*)calloc(LEN_OUTER_IN, sizeof(float)); 
            float *L_outer_in_par = (float*)calloc(LEN_OUTER_IN, sizeof(float)); 
            
            float *L_outer_ext_sys = (float*)calloc(LEN_OUTER_IN, sizeof(float));
            float *L_outer_ext_par = (float*)calloc(LEN_OUTER_IN, sizeof(float));
            
            float *L_apriori_inner = (float*)malloc(LEN_INNER_IN * sizeof(float));

            int max_len = (LEN_INNER_IN > LEN_OUTER_IN) ? LEN_INNER_IN : LEN_OUTER_IN;
            float (*alpha)[STATE_NUM] = malloc((max_len + 1) * sizeof(*alpha));
            float (*beta)[STATE_NUM] = malloc((max_len + 1) * sizeof(*beta));
            float (*gamma)[STATE_NUM][2] = malloc(max_len * sizeof(*gamma));

            rng_state_t rng;
            unsigned int seed = (unsigned int)(time(0) ^ omp_get_thread_num() ^ (unsigned int)(current_snr * 1000));
            rng_init(&rng, seed);

            #pragma omp for schedule(dynamic)
            for (long int seq = 0; seq < seq_num; seq++)
            {
                int i, iter;

                // 1. 生成数据
                for (i = 0; i < INFO_LEN; i++) msg[i] = rng_next(&rng) % 2;
                for (i = INFO_LEN; i < LEN_OUTER_IN; i++) msg[i] = 0; 

                // 2. 外编码
                outer_encoder(msg, outer_coded, LEN_OUTER_IN);
                
                // 3. 交织
                for (i = 0; i < LEN_INTER; i++) inner_in[i] = outer_coded[sccc_interleaver[i]];
                for (i = LEN_INTER; i < LEN_INNER_IN; i++) inner_in[i] = 0; 

                // 4. 内编码 + 信道
                int s = 0;
                for (i = 0; i < LEN_INNER_IN; i++) {
                    int u = inner_in[i];      
                    int p = output_parity[s][u]; 
                    s = next_state[s][u];        

                    float u1 = rng_float(&rng); if(u1<1e-6f) u1=1e-6f;
                    float u2 = rng_float(&rng);
                    float g = sgm * sqrtf(-2.0f * logf(u1)) * cosf(2.0f * PI * u2);
                    float tx = (u == 0) ? 1.0f : -1.0f;
                    rx_inner_sys[i] = (tx + g) * Lc;

                    u1 = rng_float(&rng); if(u1<1e-6f) u1=1e-6f;
                    u2 = rng_float(&rng);
                    g = sgm * sqrtf(-2.0f * logf(u1)) * cosf(2.0f * PI * u2);
                    tx = (p == 0) ? 1.0f : -1.0f;
                    rx_inner_par[i] = (tx + g) * Lc;
                    
                    L_apriori_inner[i] = 0.0f;
                }

                // 5. 迭代解码
                for (iter = 0; iter < ITERATIONS; iter++) {
                    siso_generic(rx_inner_sys, rx_inner_par, L_apriori_inner, NULL, 
                                 L_inner_ext, NULL, LEN_INNER_IN, 0, 0, 
                                 alpha, beta, gamma);

                    for (i = 0; i < LEN_INTER; i++) {
                        int outer_idx = sccc_interleaver[i];
                        float val = L_inner_ext[i];
                        if (outer_idx % 2 == 0) L_outer_in_sys[outer_idx / 2] = val;
                        else                    L_outer_in_par[outer_idx / 2] = val;
                    }

                    siso_generic(L_outer_in_sys, L_outer_in_par, NULL, NULL,
                                 L_outer_ext_sys, L_outer_ext_par, LEN_OUTER_IN, 0, 1,
                                 alpha, beta, gamma);

                    for (i = 0; i < LEN_INTER; i++) {
                        int outer_idx = sccc_interleaver[i];
                        if (outer_idx % 2 == 0) L_apriori_inner[i] = L_outer_ext_sys[outer_idx / 2];
                        else                    L_apriori_inner[i] = L_outer_ext_par[outer_idx / 2];
                    }
                    for (i = LEN_INTER; i < LEN_INNER_IN; i++) L_apriori_inner[i] = 0.0f;
                }

                // 6. 判决
                long local_errors = 0;
                for (i = 0; i < INFO_LEN; i++) {
                    float final_LLR = L_outer_in_sys[i] + L_outer_ext_sys[i];
                    int dec_bit = (final_LLR >= 0) ? 0 : 1;
                    if (msg[i] != dec_bit) local_errors++;
                }
                total_bit_error += local_errors;
            }
            // 释放内存
            free(msg); free(outer_coded); free(inner_in);
            free(rx_inner_sys); free(rx_inner_par);
            free(L_inner_ext); free(L_outer_in_sys); free(L_outer_in_par);
            free(L_outer_ext_sys); free(L_outer_ext_par); free(L_apriori_inner);
            free(alpha); free(beta); free(gamma);
        }

        double t_end = omp_get_wtime();
        double BER = (double)total_bit_error / (double)(INFO_LEN * seq_num);
        
        // 1. 打印到屏幕
        printf("SNR=%4.1f | Frames=%ld | Time=%.2fs | BER=%E\n", 
               current_snr, seq_num, t_end - t_start, BER);
        
        // 2. 写入到 CSV 文件
        fprintf(fp, "%.2f,%.6E\n", current_snr, BER);
        // 3. 强制刷新缓冲区 (防止程序中断数据丢失)
        fflush(fp);
    }
    
    // 关闭文件
    fclose(fp);
    printf("\nData saved to sccc_ber.csv. Simulation finished.\n");
    return 0;
}

// --- 辅助函数 ---

void statetable() {
    int s, u;
    for (s = 0; s < STATE_NUM; s++) {
        for (u = 0; u < 2; u++) {
            int m1 = (s >> 1) & 1; int m0 = s & 1;        
            int a_k = (u + m1 + m0) % 2;
            int out = (a_k + m0) % 2;
            next_state[s][u] = (a_k << 1) | m1;
            output_parity[s][u] = out;
        }
    }
    srand(12345);
    for (int i = 0; i < LEN_INTER; i++) sccc_interleaver[i] = i;
    for (int i = LEN_INTER - 1; i > 0; i--) {
        int j = rand() % (i + 1);
        int temp = sccc_interleaver[i];
        sccc_interleaver[i] = sccc_interleaver[j];
        sccc_interleaver[j] = temp;
    }
}

void outer_encoder(int* input_bits, int* output_stream, int len) {
    int s = 0;
    for (int i = 0; i < len; i++) {
        int u = input_bits[i];
        output_stream[2*i] = u; 
        output_stream[2*i+1] = output_parity[s][u];
        s = next_state[s][u];
    }
}

void siso_generic(float* L_in_sys, float* L_in_par, float* L_in_apriori_sys, float* L_in_apriori_par,
                  float* L_out_sys, float* L_out_par, int len, int terminate, int mode,
                  float (*alpha)[STATE_NUM], float (*beta)[STATE_NUM], float (*gamma)[STATE_NUM][2])
{
    int k, s, u;

    // Alpha Init
    for (s = 0; s < STATE_NUM; s++) alpha[0][s] = (s == 0) ? 0.0f : -INF;

    // Beta Init
    for (s = 0; s < STATE_NUM; s++) {
        if (terminate) beta[len][s] = (s == 0) ? 0.0f : -INF;
        else beta[len][s] = 0.0f; 
    }

    // Gamma & Alpha
    for (k = 0; k < len; k++) {
        float L_a_sys = (L_in_apriori_sys) ? L_in_apriori_sys[k] : 0.0f;
        float L_a_par = (L_in_apriori_par) ? L_in_apriori_par[k] : 0.0f;
        
        for (s = 0; s < STATE_NUM; s++) {
            for (u = 0; u < 2; u++) {
                int p = output_parity[s][u];
                float sign_u = (u == 0) ? 1.0f : -1.0f;
                float sign_p = (p == 0) ? 1.0f : -1.0f;
                gamma[k][s][u] = 0.5f * (sign_u * (L_in_sys[k] + L_a_sys) + sign_p * (L_in_par[k] + L_a_par));
            }
        }
        
        for (s = 0; s < STATE_NUM; s++) {
            alpha[k+1][s] = -INF;
            for(int prev_s = 0; prev_s < STATE_NUM; prev_s++) {
                for(int input_bit = 0; input_bit < 2; input_bit++) {
                    if (next_state[prev_s][input_bit] == s) {
                        alpha[k+1][s] = max_star(alpha[k+1][s], alpha[k][prev_s] + gamma[k][prev_s][input_bit]);
                    }
                }
            }
        }
    }

    // Beta
    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]);
            }
        }
    }

    // Extrinsic
    for (k = 0; k < len; k++) {
        float L0 = -INF, L1 = -INF;
        for (s = 0; s < STATE_NUM; s++) {
            L0 = max_star(L0, alpha[k][s] + gamma[k][s][0] + beta[k+1][next_state[s][0]]);
            L1 = max_star(L1, alpha[k][s] + gamma[k][s][1] + beta[k+1][next_state[s][1]]);
        }
        float L_a_sys = (L_in_apriori_sys) ? L_in_apriori_sys[k] : 0.0f;
        float ext_sys = (L0 - L1) - L_in_sys[k] - L_a_sys;
        if(ext_sys > 50.0f) ext_sys = 50.0f; else if(ext_sys < -50.0f) ext_sys = -50.0f;
        L_out_sys[k] = ext_sys;

        if (mode == 1 && L_out_par != NULL) {
            float Lp0 = -INF, Lp1 = -INF;
            for (s = 0; s < STATE_NUM; s++) {
                for(u = 0; u < 2; u++) {
                    int p = output_parity[s][u];
                    float metric = alpha[k][s] + gamma[k][s][u] + beta[k+1][next_state[s][u]];
                    if (p == 0) Lp0 = max_star(Lp0, metric);
                    else        Lp1 = max_star(Lp1, metric);
                }
            }
            float L_a_par = (L_in_apriori_par) ? L_in_apriori_par[k] : 0.0f;
            float ext_par = (Lp0 - Lp1) - L_in_par[k] - L_a_par;
            if(ext_par > 50.0f) ext_par = 50.0f; else if(ext_par < -50.0f) ext_par = -50.0f;
            L_out_par[k] = ext_par;
        }
    }
}