Turbo_code_version1.cpp

/***************************************************
Channel Coding Course Work: Turbo Codes (Final Working Version)
Fixes:
1. Corrected De-interleaving mapping logic (CRITICAL FIX).
2. Removed strict termination assumption in decoder (Beta init).
3. Corrected data types and print formats.
***************************************************/

#define  _CRT_SECURE_NO_WARNINGS
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>

// --- System Parameters ---
#define MSG_LEN 256    // Message length
#define STATE_NUM 4    // States for RSC(1, 5/7)
#define ITERATIONS 8   // Iterations
#define TAIL_BITS 2    

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

float code_rate = 1.0f / 3.0f;

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

double N0, sgm;

// Trellis Structures
int next_state[STATE_NUM][2]; 
int output_parity[STATE_NUM][2]; 

// Interleaver Map
int alpha_interleaver[message_length];
// int deinterleaver[message_length]; // REMOVED: Not needed and caused confusion

// Data Buffers
int message[message_length];
int codeword[codeword_length];
int re_codeword[codeword_length];
int de_message[message_length];

// Modulation Buffers
double tx_symbol[codeword_length][2];
double rx_symbol[codeword_length][2];

// Decoder Internal Buffers
double rx_sys[message_length]; 
double rx_par1[message_length]; 
double rx_par2[message_length]; 

// Extrinsic Information exchange
double L_ext1[message_length]; 
double L_ext2[message_length]; 

// Function Prototypes
void statetable();
void encoder();
void modulation();
void channel();
void decoder();
void rsc_encoder(int* input_bits, int* parity_out, int len);
void siso_decoder(double* L_in_sys, double* L_in_par, double* L_in_apriori, double* L_out_extrinsic, int len, int terminate);
double max_star(double a, double b);

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

    statetable();
    srand((unsigned int)time(0));

    printf("--- Turbo Code Simulation (Final Version) ---\n");
    printf("Message Length: %d, Codeword Length: %d\n", message_length, codeword_length);
    
    printf("\nEnter start SNR (dB) [suggest -2.0]: ");
    scanf("%f", &start);
    
    printf("\nEnter finish SNR (dB) [suggest 3.0]: ");
    scanf("%f", &finish);
    
    printf("\nPlease input number of frames [e.g., 50]: ");
    scanf("%ld", &seq_num);

    for (SNR = start; SNR <= finish; SNR += 0.5f)
    {
        N0 = (1.0 / code_rate) / pow(10.0, (float)(SNR) / 10.0);
        sgm = sqrt(N0 / 2.0);

        bit_error = 0;

        for (seq = 1; seq <= seq_num; seq++)
        {
            // 1. Generate Message (Last 2 bits 0)
            for (i = 0; i < message_length - TAIL_BITS; i++)
                message[i] = rand() % 2;
            for (i = message_length - TAIL_BITS; i < message_length; i++)
                message[i] = 0; 

            // 2. Encode
            encoder();

            // 3. Modulate
            modulation();

            // 4. Channel
            channel();

            // 5. Decode
            decoder();

            // 6. Count Errors
            for (i = 0; i < message_length - TAIL_BITS; i++)
            {
                if (message[i] != de_message[i])
                    bit_error++;
            }

            progress = (double)(seq * 100) / (double)seq_num;
            if (seq % 10 == 0 || seq == seq_num) {
                BER = (double)bit_error / (double)((message_length - TAIL_BITS) * seq);
                printf("Progress=%3.0f%%, SNR=%4.1f, Errors=%6ld, BER=%E\r", progress, SNR, bit_error, BER);
            }
        }
        BER = (double)bit_error / (double)((message_length - TAIL_BITS) * seq_num);
        printf("Progress=%3.0f%%, SNR=%4.1f, Errors=%6ld, BER=%E\n", progress, SNR, bit_error, BER);
    }
    
    printf("\nSimulation finished. Press Enter to exit.");
    getchar(); getchar();
    return 0;
}

void statetable()
{
    int s, u;
    // RSC (1, 5/7) octal
    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;
        }
    }

    // Random Interleaver
    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;
    int i;
    for (i = 0; i < len; i++) {
        int u = input_bits[i];
        parity_out[i] = output_parity[s][u];
        s = next_state[s][u];
    }
}

void encoder()
{
    int i;
    int parity1[message_length];
    int parity2[message_length];
    int interleaved_msg[message_length];

    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);

    for (i = 0; i < message_length; i++) {
        codeword[3 * i] = message[i];
        codeword[3 * i + 1] = parity1[i];
        codeword[3 * i + 2] = parity2[i];
    }
}

void modulation()
{
    // BPSK: 0->+1, 1->-1
    int i;
    for (i = 0; i < codeword_length; i++)
    {
        tx_symbol[i][0] = (codeword[i] == 0) ? 1.0 : -1.0;
        tx_symbol[i][1] = 0;
    }
}

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

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

void siso_decoder(double* L_in_sys, double* L_in_par, double* L_in_apriori, double* L_out_extrinsic, int len, int terminate)
{
    static double alpha[MSG_LEN + TAIL_BITS + 1][STATE_NUM];
    static double beta[MSG_LEN + TAIL_BITS + 1][STATE_NUM];
    static double gamma[MSG_LEN + TAIL_BITS][STATE_NUM][2];

    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
    // Assuming unknown termination (equal probability)
    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;
                
                // LLR: Positive supports 0.
                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
    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
    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];
        
        // Clamp to avoid overflow/instability
        if (ext > 50.0) ext = 50.0;
        if (ext < -50.0) ext = -50.0;
        L_out_extrinsic[k] = ext;
    }
}

void decoder()
{
    int i, iter;
    double Lc = 2.0 / (sgm * sgm); 

    // Initialize inputs
    for (i = 0; i < message_length; i++) {
        rx_sys[i]  = Lc * rx_symbol[3 * i][0];
        rx_par1[i] = Lc * rx_symbol[3 * i + 1][0];
        rx_par2[i] = Lc * rx_symbol[3 * i + 2][0];
        L_ext2[i] = 0.0; 
    }

    for (iter = 0; iter < ITERATIONS; iter++) {
        // --- Decoder 1 ---
        siso_decoder(rx_sys, rx_par1, L_ext2, L_ext1, message_length, 0);

        // Interleave Ext1 -> Priori2
        double rx_sys_int[message_length];
        double L_apriori2[message_length];
        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);

        // De-interleave Ext2 -> Priori1
        // FIX: Mapping back from interleaved domain to natural domain
        // Interleaved position 'i' corresponds to Natural position 'alpha[i]'
        double temp[message_length];
        for (i = 0; i < message_length; i++) temp[i] = L_ext2[i]; // Store interleaved output
        
        for (i = 0; i < message_length; i++) {
            // Put the info back to where it belongs in the natural sequence
            L_ext2[alpha_interleaver[i]] = temp[i];
        }
    }

    // Hard Decision
    for (i = 0; i < message_length; i++) {
        double final_LLR = rx_sys[i] + L_ext1[i] + L_ext2[i];
        if (final_LLR >= 0) de_message[i] = 0;
        else de_message[i] = 1;
    }
}

void demodulation() {}