pccc_final_csv.c

/***************************************************
Channel Coding Course Work: Turbo Codes (Final CSV Version)
功能:
1. PCCC (Turbo Code) 仿真
2. 高性能优化 (OpenMP, 内存优化)
3. 结果自动保存到 turbo_ber.csv
***************************************************/

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

// --- 系统参数 ---
#define MSG_LEN 256    
#define STATE_NUM 4    
#define ITERATIONS 4   
#define TAIL_BITS 0    

#define MESSAGE_LENGTH (MSG_LEN + TAIL_BITS) 
#define CODEWORD_LENGTH (MESSAGE_LENGTH * 3) 

float code_rate = 1.0f / 3.0f;

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

// --- 全局变量 (只读) ---
int next_state[STATE_NUM][2]; 
int output_parity[STATE_NUM][2]; 
int alpha_interleaver[MESSAGE_LENGTH];

// --- 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;
}

// --- Max-Star (移到上方并内联，提高速度) ---
static inline float max_star(float a, float b) {
    // PCCC 通常使用 Max-Log-MAP (简单取最大值) 即可收敛良好
    // 如果需要极致性能，可以像 SCCC 那样加 Log-MAP 修正表，但 PCCC 对此不敏感
    return (a > b) ? a : b;
}

// --- 函数原型 ---
void statetable();
void rsc_encoder(int* input_bits, int* parity_out, int len);
// max_star 已在上方定义为 static inline，不需要原型
void siso_decoder(float* L_in_sys, float* L_in_par, float* L_in_apriori, float* L_out_extrinsic, int len, int terminate, 
                 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();

    // 2. 打开 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"); // 写入表头

    printf("--- Turbo Code Simulation (Multi-Core + CSV Output) ---\n");
    printf("Message Length: %d, Codeword Length: %d\n", MESSAGE_LENGTH, CODEWORD_LENGTH);
    printf("Detected Processors: %d\n", omp_get_num_procs());
    printf("Result will be saved to 'turbo_ber.csv'.\n");
    
    printf("\nEnter start SNR (dB) [suggest -2.0]: ");
    if(scanf("%f", &start_snr));
    printf("\nEnter finish SNR (dB) [suggest 3.0]: ");
    if(scanf("%f", &finish_snr));
    printf("\nPlease input number of frames [e.g., 10000]: ");
    if(scanf("%ld", &seq_num));

    for (current_snr = start_snr; current_snr <= finish_snr; current_snr += 0.2f) // 步长建议 0.2
    {
        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)
        {
            // 内存分配 (使用 calloc 初始化为 0)
            int *message = (int*)calloc(MESSAGE_LENGTH, sizeof(int));
            int *parity1 = (int*)malloc(MESSAGE_LENGTH * sizeof(int));
            int *parity2 = (int*)malloc(MESSAGE_LENGTH * sizeof(int));
            int *interleaved_msg = (int*)malloc(MESSAGE_LENGTH * sizeof(int));
            
            // 接收缓冲区
            float *rx_sys = (float*)malloc(MESSAGE_LENGTH * sizeof(float));
            float *rx_par1 = (float*)malloc(MESSAGE_LENGTH * sizeof(float));
            float *rx_par2 = (float*)malloc(MESSAGE_LENGTH * sizeof(float));
            
            // 解码缓冲区
            float *L_ext1 = (float*)calloc(MESSAGE_LENGTH, sizeof(float));
            float *L_ext2 = (float*)calloc(MESSAGE_LENGTH, sizeof(float));
            float *L_apriori2 = (float*)calloc(MESSAGE_LENGTH, sizeof(float));
            float *rx_sys_int = (float*)malloc(MESSAGE_LENGTH * sizeof(float));
            float *temp_buf = (float*)malloc(MESSAGE_LENGTH * sizeof(float));

            // SISO 递归数组
            float (*alpha)[STATE_NUM] = malloc((MESSAGE_LENGTH + 1) * sizeof(*alpha));
            float (*beta)[STATE_NUM] = malloc((MESSAGE_LENGTH + 1) * sizeof(*beta));
            float (*gamma)[STATE_NUM][2] = malloc(MESSAGE_LENGTH * sizeof(*gamma));

            rng_state_t rng;
            // 种子混入 SNR 以保证独立性
            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 < MESSAGE_LENGTH - TAIL_BITS; i++)
                    message[i] = rng_next(&rng) % 2;
                for (i = MESSAGE_LENGTH - TAIL_BITS; i < MESSAGE_LENGTH; i++)
                    message[i] = 0;

                // 2. 编码
                for(i=0; i<MESSAGE_LENGTH; i++) interleaved_msg[i] = message[alpha_interleaver[i]];
                rsc_encoder(message, parity1, MESSAGE_LENGTH);
                rsc_encoder(interleaved_msg, parity2, MESSAGE_LENGTH);

                // 3. 调制 + 信道 (BPSK + AWGN)
                for (i = 0; i < MESSAGE_LENGTH; i++) {
                    // System bit
                    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 = (message[i] == 0) ? 1.0f : -1.0f;
                    rx_sys[i] = (tx + g) * Lc;

                    // Parity 1
                    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 = (parity1[i] == 0) ? 1.0f : -1.0f;
                    rx_par1[i] = (tx + g) * Lc;

                    // Parity 2
                    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 = (parity2[i] == 0) ? 1.0f : -1.0f;
                    rx_par2[i] = (tx + g) * Lc;
                    
                    L_ext2[i] = 0.0f; // 重置外信息
                }

                // 4. 解码
                for (iter = 0; iter < ITERATIONS; iter++) {
                    // Decoder 1
                    siso_decoder(rx_sys, rx_par1, L_ext2, L_ext1, MESSAGE_LENGTH, 0, alpha, beta, gamma);

                    // Interleave
                    for (i = 0; i < MESSAGE_LENGTH; i++) {
                        rx_sys_int[i] = rx_sys[alpha_interleaver[i]];
                        L_apriori2[i] = L_ext1[alpha_interleaver[i]];
                    }

                    // Decoder 2
                    siso_decoder(rx_sys_int, rx_par2, L_apriori2, L_ext2, MESSAGE_LENGTH, 0, alpha, beta, gamma);

                    // De-interleave
                    for (i = 0; i < MESSAGE_LENGTH; i++) temp_buf[i] = L_ext2[i];
                    for (i = 0; i < MESSAGE_LENGTH; i++) L_ext2[alpha_interleaver[i]] = temp_buf[i];
                }

                // 5. 判决
                long local_errors = 0;
                for (i = 0; i < MESSAGE_LENGTH - TAIL_BITS; i++) {
                    float final_LLR = rx_sys[i] + L_ext1[i] + L_ext2[i];
                    int dec_bit = (final_LLR >= 0) ? 0 : 1;
                    if (message[i] != dec_bit) local_errors++;
                }
                total_bit_error += local_errors;
            }

            // 释放内存
            free(message); free(parity1); free(parity2); free(interleaved_msg);
            free(rx_sys); free(rx_par1); free(rx_par2);
            free(L_ext1); free(L_ext2); free(L_apriori2); free(rx_sys_int); free(temp_buf);
            free(alpha); free(beta); free(gamma);
        }

        double t_end = omp_get_wtime();
        double BER = (double)total_bit_error / (double)((MESSAGE_LENGTH - TAIL_BITS) * seq_num);
        
        // 打印与保存
        printf("SNR=%4.1f | Frames=%ld | Time=%.2fs | BER=%E\n", 
               current_snr, seq_num, t_end - t_start, BER);
        
        fprintf(fp, "%.2f,%.6E\n", current_snr, BER);
        fflush(fp);
    }
    
    fclose(fp);
    printf("\nSimulation finished. Data saved to turbo_ber.csv.\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);
    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;
    }
}

void rsc_encoder(int* input_bits, int* parity_out, int len)
{
    int s = 0;
    for (int i = 0; i < len; i++) {
        int u = input_bits[i];
        parity_out[i] = output_parity[s][u];
        s = next_state[s][u];
    }
}

void siso_decoder(float* L_in_sys, float* L_in_par, float* L_in_apriori, float* L_out_extrinsic, int len, int terminate,
                 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++) {
        // PCCC 通常假设终止于状态0，或不假设
        // 这里沿用之前的逻辑：terminate=0 -> 平坦分布
        beta[len][s] = (terminate && s==0) ? 0.0f : (terminate ? -INF : 0.0f);
    }

    // Gamma & Alpha
    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];
                float sign_u = (u == 0) ? 1.0f : -1.0f;
                float sign_p = (parity == 0) ? 1.0f : -1.0f;
                gamma[k][s][u] = 0.5f * (sign_u * (L_in_sys[k] + L_in_apriori[k]) + sign_p * L_in_par[k]);
            }
        }
        
        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;
        float 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 ext = (L0 - L1) - L_in_sys[k] - L_in_apriori[k];
        
        if (ext > 50.0f) ext = 50.0f;
        else if (ext < -50.0f) ext = -50.0f;
        L_out_extrinsic[k] = ext;
    }
}