fast_turbo.cpp

#define  _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <vector>
#include <random>
#include <omp.h> // OpenMP 头文件

// --- System Parameters ---
#define MSG_LEN 16384   
#define STATE_NUM 4    
#define ITERATIONS 8   
#define TAIL_BITS 4      // 保持你要求的 4 位尾比特

#define message_length (MSG_LEN + TAIL_BITS) 
#define codeword_length (message_length * 3) 

float code_rate = 1.0f / 3.0f;

// --- Global Constants ---
#define pi 3.141592653589793
#define INF 1e9        

// 只保留只读的全局查表变量，其他全部移入局部变量以支持并行
int next_state[STATE_NUM][2]; 
int output_parity[STATE_NUM][2]; 
int alpha_interleaver[message_length];

// --- 线程安全的随机数生成器 ---
// 替代 rand()，rand() 在多线程中会锁住导致性能下降
int get_rand_bit() {
    // 每个线程拥有独立的随机数引擎
    static thread_local std::mt19937 gen(std::random_device{}());
    std::uniform_int_distribution<> dis(0, 1);
    return dis(gen);
}

// 线程安全的高斯噪声生成 (Box-Muller)
void get_gaussian_noise(double sgm, double* n1) {
    static thread_local std::mt19937 gen(std::random_device{}());
    static thread_local std::uniform_real_distribution<double> dis(0.0, 1.0);

    double u1 = dis(gen);
    if (u1 >= 1.0) u1 = 0.999999;
    double u2 = dis(gen);
    if (u2 >= 1.0) u2 = 0.999999;
    
    double r = sgm * sqrt(2.0 * log(1.0 / (1.0 - u1)));
    *n1 = r * cos(2 * pi * u2); 
}

// Function Prototypes
void statetable();
void rsc_encoder(const int* input_bits, int* parity_out, int len);
double max_star(double a, double b);

// --- 核心 SISO 解码器 (并行改造版) ---
// 增加了 alpha_buf, beta_buf, gamma_buf 参数，避免在函数内使用 static
void siso_decoder(const double* L_in_sys, const double* L_in_par, const double* L_in_apriori, 
                  double* L_out_extrinsic, int len, int terminate,
                  std::vector<double>& alpha_buf, std::vector<double>& beta_buf, std::vector<double>& gamma_buf)
{
    // 利用传入的 vector 内存作为数组使用
    // alpha 映射到 alpha_buf
    double (*alpha)[STATE_NUM] = (double (*)[STATE_NUM])alpha_buf.data();
    // beta 映射到 beta_buf
    double (*beta)[STATE_NUM]  = (double (*)[STATE_NUM])beta_buf.data();
    // gamma 映射到 gamma_buf
    double (*gamma)[STATE_NUM][2] = (double (*)[STATE_NUM][2])gamma_buf.data();

    int k, s, u;

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

    // 2. Init Beta
    // 根据 terminate 参数决定边界条件
    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 Calculation
    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];
                double sign_u = (u == 0) ? 1.0 : -1.0;
                double sign_p = (parity == 0) ? 1.0 : -1.0;
                // Max-Log-MAP 分支度量
                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 Recursion (Forward)
    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 Recursion (Backward)
    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
    for (k = 0; k < len; k++) {
        double L0 = -INF;
        double L1 = -INF;

        for (s = 0; s < STATE_NUM; s++) {
            int ns0 = next_state[s][0];
            L0 = max_star(L0, alpha[k][s] + gamma[k][s][0] + beta[k+1][ns0]);
            
            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;
        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;
    }
}

double max_star(double a, double b) {
    return (a > b) ? a : b;
}

void rsc_encoder(const 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 statetable()
{
    int s, u;
    // 构建 RSC (1, 5/7)
    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;
            int next_s = (a_k << 1) | m1;
            next_state[s][u] = next_s;
            output_parity[s][u] = out;
        }
    }
    // Knuth shuffle 生成交织器
    srand((unsigned int)time(0)); // 仅用于初始化交织器
    for (int i = 0; i < message_length; i++) alpha_interleaver[i] = i;
    for (int i = message_length - 1; i > 0; i--) {
        int j = rand() % (i + 1);
        int temp = alpha_interleaver[i];
        alpha_interleaver[i] = alpha_interleaver[j];
        alpha_interleaver[j] = temp;
    }
}

int main()
{
    float SNR, start, finish, step = 0.2f;
    long int seq_num;
    double N0, sgm;
    FILE *fp = NULL;

    statetable(); // 初始化

    // OpenMP 核心数检测
    int num_procs = omp_get_num_procs();
    printf("Detected Processors: %d\n", num_procs);
    omp_set_num_threads(num_procs);

    fp = fopen("turbo_ber_final.csv", "w");
    if (fp) fprintf(fp, "SNR,BER\n");
    else { printf("File open error!\n"); return -1; }

    printf("--- Turbo Code Parallel Simulation (Corrected Logic) ---\n");
    printf("Enter start SNR (dB): "); scanf("%f", &start);
    printf("Enter finish SNR (dB): "); scanf("%f", &finish);
    printf("Enter number of frames: "); 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.0);
        long int total_bit_error = 0;

        // --- OpenMP 并行区域 ---
        // 每个线程处理一部分帧，最后汇总 total_bit_error
        #pragma omp parallel for reduction(+:total_bit_error) schedule(dynamic)
        for (long int seq = 1; seq <= seq_num; seq++)
        {
            // --- 线程局部变量 (移出全局，确保安全) ---
            int message[message_length];
            int codeword[codeword_length];
            int parity1[message_length];
            int parity2[message_length];
            int interleaved_msg[message_length];
            
            double rx_sys[message_length]; 
            double rx_par1[message_length]; 
            double rx_par2[message_length]; 
            double L_ext1[message_length]; 
            double L_ext2[message_length]; 
            double rx_sys_int[message_length];
            double L_apriori2[message_length];
            double temp_ext[message_length];

            // 预分配 SISO 计算所需的内存，避免 malloc/free 开销
            std::vector<double> alpha_buf((message_length + 1) * STATE_NUM);
            std::vector<double> beta_buf((message_length + 1) * STATE_NUM);
            std::vector<double> gamma_buf(message_length * STATE_NUM * 2);

            int i;

            // 1. 生成数据与计算尾比特 (归零逻辑)
            // 先随机生成数据位
            for (i = 0; i < message_length - TAIL_BITS; i++)
                message[i] = get_rand_bit();
            
            // 归零计算：模拟 Encoder 1 的状态
            int s = 0; 
            int u, m1, m0, a_k;
            // 跑完前段数据
            for (i = 0; i < message_length - TAIL_BITS; i++) {
                u = message[i];
                m1 = (s >> 1) & 1; m0 = s & 1;       
                a_k = (u + m1 + m0) % 2; 
                s = (a_k << 1) | m1;
            }
            // 计算尾比特让状态归零
            for (i = message_length - TAIL_BITS; i < message_length; i++) {
                m1 = (s >> 1) & 1; m0 = s & 1;
                message[i] = (m1 + m0) % 2; // 你的归零公式
                
                // Update State
                u = message[i];
                a_k = (u + m1 + m0) % 2; 
                s = (a_k << 1) | m1;
            }

            // 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. 调制与过信道 (合并计算，提升缓存命中率)
            double n;
            double Lc = 2.0 / (sgm * sgm);
            
            // 处理系统位 (Systematic)
            for (i = 0; i < message_length; i++) {
                double tx = (message[i] == 0) ? 1.0 : -1.0;
                get_gaussian_noise(sgm, &n);
                rx_sys[i] = Lc * (tx + n);
                L_ext2[i] = 0.0; // 初始化外信息
            }
            // 处理校验位 1
            for (i = 0; i < message_length; i++) {
                double tx = (parity1[i] == 0) ? 1.0 : -1.0;
                get_gaussian_noise(sgm, &n);
                rx_par1[i] = Lc * (tx + n);
            }
            // 处理校验位 2
            for (i = 0; i < message_length; i++) {
                double tx = (parity2[i] == 0) ? 1.0 : -1.0;
                get_gaussian_noise(sgm, &n);
                rx_par2[i] = Lc * (tx + n);
            }

            // 4. 迭代解码
            for (int iter = 0; iter < ITERATIONS; iter++) {
                // --- DEC 1 (自然顺序) ---
                // 注意：terminate = 1 (因为我们做了归零处理)
                siso_decoder(rx_sys, rx_par1, L_ext2, L_ext1, message_length, 1, alpha_buf, beta_buf, gamma_buf);

                // 交织
                for (i = 0; i < message_length; i++) {
                    rx_sys_int[i] = rx_sys[alpha_interleaver[i]];
                    L_apriori2[i] = L_ext1[alpha_interleaver[i]];
                }

                // --- DEC 2 (交织顺序) ---
                // 注意：terminate = 0 (必须为0，交织后状态未知)
                siso_decoder(rx_sys_int, rx_par2, L_apriori2, L_ext2, message_length, 0, alpha_buf, beta_buf, gamma_buf);

                // 解交织
                for (i = 0; i < message_length; i++) temp_ext[i] = L_ext2[i];
                for (i = 0; i < message_length; i++) L_ext2[alpha_interleaver[i]] = temp_ext[i];
            }

            // 5. 误码统计
            // 不统计尾比特
            for (i = 0; i < message_length - TAIL_BITS; i++) {
                double final_LLR = rx_sys[i] + L_ext1[i] + L_ext2[i];
                int dec_bit = (final_LLR >= 0) ? 0 : 1;
                if (message[i] != dec_bit)
                    total_bit_error++;
            }
        }
        // --- End Parallel Region ---

        double final_BER = (double)total_bit_error / (double)((message_length - TAIL_BITS) * seq_num);
        if (fp) fprintf(fp, "%f,%E\n", SNR, final_BER);
        
        printf("SNR=%4.1f, Errors=%6ld, BER=%E\n", SNR, total_bit_error, final_BER);
    }
    
    if (fp) fclose(fp);
    printf("Simulation Finished.\n");
    getchar(); getchar();
    return 0;
}