turbo_opt.c

/***************************************************
Channel Coding Course Work: Turbo Codes (Optimized Version)
Optimizations:
1. OpenMP Multi-threading (Uses all CPU cores).
2. Thread-safe Random Number Generator (XorShift).
3. Memory locality optimization (Removed globals).
4. Float precision (Faster SIMD on M2).
***************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include <omp.h> // Header for OpenMP

// --- System Parameters ---
#define MSG_LEN 256    
#define STATE_NUM 4    
#define ITERATIONS 8   
#define TAIL_BITS 2    

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

float code_rate = 1.0f / 3.0f;

// --- Constants ---
#define PI 3.141592653589793f
#define INF 1e9f        

// Global Tables (Read-only after init, so thread-safe)
int next_state[STATE_NUM][2]; 
int output_parity[STATE_NUM][2]; 
int alpha_interleaver[MESSAGE_LENGTH];

// --- XorShift RNG (Fast & Thread-safe) ---
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));
}

// Generates uniform float [0, 1)
float rng_float(rng_state_t *rng) {
    return (float)rng_next(rng) / (float)0xFFFFFFFF;
}

// --- Function Prototypes ---
void statetable();
void rsc_encoder(int* input_bits, int* parity_out, int len);
float max_star(float a, float b);
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]);

// --- Main ---
int main()
{
    float start_snr, finish_snr, current_snr;
    long int seq_num;

    statetable();

    printf("--- Turbo Code Simulation (Multi-Core Optimized) ---\n");
    printf("Message Length: %d, Codeword Length: %d\n", MESSAGE_LENGTH, CODEWORD_LENGTH);
    printf("Detected Processors: %d\n", omp_get_num_procs());
    
    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.1f)
    {
        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;
        
        // Start Timer
        double t_start = omp_get_wtime();

        // --- PARALLEL REGION STARTS HERE ---
        #pragma omp parallel reduction(+:total_bit_error)
        {
            // Thread-local variables (Allocated on stack or heap per thread)
            // Using Heap (malloc) to prevent Stack Overflow for large arrays
            int *message = (int*)malloc(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));
            int *de_message = (int*)malloc(MESSAGE_LENGTH * sizeof(int));
            
            float *tx_sym_0 = (float*)malloc(CODEWORD_LENGTH * sizeof(float)); // Only store dimension 0 (BPSK)
            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*)malloc(MESSAGE_LENGTH * sizeof(float));
            float *L_ext2 = (float*)malloc(MESSAGE_LENGTH * sizeof(float));
            float *L_apriori2 = (float*)malloc(MESSAGE_LENGTH * sizeof(float));
            float *rx_sys_int = (float*)malloc(MESSAGE_LENGTH * sizeof(float));
            float *temp_buf = (float*)malloc(MESSAGE_LENGTH * sizeof(float));

            // Decoder buffers (Thread local to avoid static race condition)
            // Using flat memory block for 2D/3D arrays to simplify malloc
            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));

            // Initialize Thread-safe RNG
            rng_state_t rng;
            rng_init(&rng, (unsigned int)(time(0) ^ omp_get_thread_num() ^ (unsigned int)t_start));

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

                // 1. Generate Message
                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. Encode
                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. Modulation (BPSK) & Channel (AWGN)
                // Merged for speed: Calculate Rx immediately, don't store codeword
                for (i = 0; i < MESSAGE_LENGTH; i++) {
                    // System bit
                    float u1 = rng_float(&rng); if(u1<1e-6) u1=1e-6;
                    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-6) u1=1e-6;
                    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-6) u1=1e-6;
                    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;
                    
                    // Init Extrinsic
                    L_ext2[i] = 0.0f;
                }

                // 4. Decode
                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 L_ext2 for next iteration
                    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. Hard Decision & Error Count
                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;
            }

            // Cleanup Thread Memory
            free(message); free(parity1); free(parity2); free(interleaved_msg); free(de_message);
            free(tx_sym_0); 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);
        } // End Parallel Region

        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);
    }
    
    printf("\nSimulation finished.\n");
    return 0;
}

// --- Helper Functions ---

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); // Fixed seed for interleaver consistency
    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];
    }
}

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

// Optimized SISO: No static vars, passed buffers
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;

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

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

    // 3. Gamma Calculation & Alpha Recursion
    // Combined loops slightly where possible, but split here for clarity
    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]);
            }
        }
        
        // Alpha step (k+1 depends on k)
        for (s = 0; s < STATE_NUM; s++) {
            alpha[k+1][s] = -INF;
            int prev_s, input_bit;
            // Unrolling loops for specific 4-state trellis could be faster, but keeping generic
            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+1][s] = max_star(alpha[k+1][s], alpha[k][prev_s] + gamma[k][prev_s][input_bit]);
                    }
                }
            }
        }
    }

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

    // 5. Extrinsic Calculation
    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];
        
        // Clamp
        if (ext > 50.0f) ext = 50.0f;
        else if (ext < -50.0f) ext = -50.0f;
        L_out_extrinsic[k] = ext;
    }
}