BCJR_K3_75.c

/***************************************************
Channel Coding Course Work: conolutional codes
This program template has given the message generator, BPSK modulation, AWGN channel model and BPSK demodulation,
you should first determine the encoder structure, then define the message and codeword length, generate the state table, write the convolutional encoder and decoder.
 
If you have any question, please contact me via e-mail: liangjw59@mail2.sysu.edu.cn, or, yangzhj59@mail2.sysu.edu.cn.
***************************************************/

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

#define message_length 65536 //the length of message
#define codeword_length (message_length * 2) //the length of codeword
float code_rate = (float)message_length / (float)codeword_length;

// channel coefficient
#define pi 3.1415926
double N0, sgm;

int state_table[4][4];//state table, the size should be defined yourself

//加一些冗余
int state_num = 4;//the number of the state of encoder structure

int message[message_length], codeword[codeword_length];//message and codeword
int re_codeword[codeword_length];//the received codeword
int de_message[message_length];//the decoding message

double tx_symbol[codeword_length][2];//the transmitted symbols
double rx_symbol[codeword_length][2];//the received symbols

int traceback_path[message_length][4];

// --- BCJR 新增全局变量 ---
// alpha[t][s] 表示时刻 t 处于状态 s 的前向概率
// beta[t][s]  表示时刻 t 处于状态 s 的后向概率
// 长度设为 message_length + 1，因为包含时刻 0 到 T
double alpha[message_length + 1][4];
double beta[message_length + 1][4];
// -----------------------

void statetable();
void encoder();
void modulation();
void demodulation();
void channel();
void decoder();
void bcjr_decoder(); // --- 新增 BCJR 解码函数声明 ---

int main()
{
	int i;
	float SNR, start, finish, step = 0.5;//SNR代表信噪比, start代表起始信噪比, finish代表结束信噪比
	long int bit_error, seq, seq_num;//bit_error代表误码数, seq代表当前帧数, seq_num代表总帧数
	double BER;//BER代表比特误码率
	double progress;//progress是用来显示进度的变量

    // --- 定义csv文件指针 (修改文件名为 bcjr_soft.csv) ---
    FILE *fp = NULL; 
    fp = fopen("bcjr_soft.csv", "w"); 
    if (fp == NULL) {
        printf("Error: Could not open file for writing.\n");
        return 1;
    }
    // 写入表头
    fprintf(fp, "SNR,BER\n"); 
    // -------------------------

	//generate state table 生成状态转移表
	statetable();

	//random seed
	srand((int)time(0));

	//input the SNR and frame number
    printf("\n--- Soft Decision BCJR (MAP) Simulation ---\n");
	printf("Enter start SNR: ");
	scanf("%f", &start);
	printf("Enter finish SNR: ");
	scanf("%f", &finish);
	printf("Please input the number of message: ");
	scanf("%ld", &seq_num);

	for (SNR = start; SNR <= finish; SNR+=step)//外层循环，遍历不同的信噪比
	{
		//channel noise
		N0 = (1.0 / code_rate) / pow(10.0, (float)(SNR) / 10.0);
		sgm = sqrt(N0 / 2);
		
		bit_error = 0;

		for (seq = 1; seq<=seq_num; seq++)//内层循环，同一信噪比下，进行seq_num次传输
		{
			//generate binary message randomly
			/*
			Pay attention that message is appended by 0 whose number is equal to the state of encoder structure.
			*/
			for (i = 0; i<message_length - state_num; i++)
			{
				message[i] = rand() % 2;//生成随机信息0或1
			}
            // --- 注意：BCJR算法非常依赖尾比特归零 ---
            // 必须确保最后几位是0，以便 Beta递归从状态0开始
			for (i = message_length - state_num; i<message_length; i++)
			{
				message[i] = 0;//信息后面添加与编码器状态数相等个数的0
			}

			//convolutional encoder
			//卷积编码器，输入是message[]，输出是codeword[]
			encoder();

			//BPSK modulation
			//BPSK调制
			modulation();

			//AWGN channel
			//高斯白噪声信道
			channel();

			//BPSK demodulation, it's needed in hard-decision Viterbi decoder
            // --- BCJR 不需要硬判决解调，它直接使用 rx_symbol (软信息) ---
			// demodulation(); 

			//convolutional decoder
            // --- 使用 BCJR 解码替代原来的 Viterbi 解码 ---
			bcjr_decoder();
            // decoder();

			//calculate the number of bit error
			for (i = 0; i<message_length; i++)
			{
				if (message[i] != de_message[i])
					bit_error++;
			}

			progress = (double)(seq * 100) / (double)seq_num;//计算当前信噪比下进度

			//calculate the intermediate BER
			BER = (double)bit_error / (double)(message_length*seq);

			//print the intermediate result
			printf("Progress=%2.1f, SNR=%2.1f, Bit Errors=%ld, BER=%E\r", progress, SNR, bit_error, BER);
		}

		//calculate the final BER
		BER = (double)bit_error / (double)(message_length*seq_num);

        // --- 将结果写入文件 ---
        // 格式: SNR, BER
        fprintf(fp, "%f,%E\n", SNR, BER);
        // 同时在屏幕打印结果
        printf("SNR=%2.1f, Bit Errors=%ld, BER=%E\n", SNR, bit_error, BER);

		//print the final result
		//printf("Progress=%2.1f, SNR=%2.1f, Bit Errors=%ld, BER=%E\n", progress, SNR, bit_error, BER);
	}

	//system("pause");

    // --- 关闭csv文件 ---
    fclose(fp);
    printf("Result saved to bcjr_soft.csv\n");

    return 0;
}

// 生成状态表
void statetable()
{
    int s0, s1, u;
    //int current_state;
    int next_state;
    int c1, c2;

    // 遍历所有可能的当前状态 (00, 01, 10, 11)
    for (int s = 0; s < 4; s++) 
    {
        // 提取二进制位: s = (s0, s1) 也就是 高位是s0, 低位是s1
        s1 = s & 1;      // 取最低位
        s0 = (s >> 1) & 1; // 取次低位

        // 遍历所有可能的输入 (0 或 1)
        for (u = 0; u <= 1; u++) 
        {
            // 电路逻辑 (保持与你硬判决一致)
            c1 = u ^ s0 ^ s1;
            c2 = u ^ s1;
            
            // 下一个状态
            next_state = (u << 1) + s0;
            
            // 输出符号
            int output = (c1 << 1) + c2;

            // --- 填表 ---
            state_table[s][u * 2]     = next_state;
            state_table[s][u * 2 + 1] = output;
        }
    }
}

void encoder()
{
    int i;
    int current_state = 0; // 初始状态必须是 00
    int input_bit;
    int output_symbol;     // 十进制的输出 (0~3)

    for (i = 0; i < message_length; i++)
    {
        input_bit = message[i];

        // 1. 根据输入位，决定查哪两列
        int col_next_state = (input_bit == 0) ? 0 : 2;
        int col_output     = (input_bit == 0) ? 1 : 3;

        // 2. 查表获取输出符号和下一状态
        output_symbol = state_table[current_state][col_output];
        int next_state = state_table[current_state][col_next_state];

        // 3. 将十进制输出符号转换为两个二进制位存入 codeword
        codeword[2 * i] = (output_symbol >> 1) & 1;
        codeword[2 * i + 1] = output_symbol & 1;

        // 4. 更新状态
        current_state = next_state;
    }
}

void modulation()
{
	//BPSK modulation
	int i;

	//0 is mapped to (1,0) and 1 is mapped tp (-1,0)
	for (i = 0; i<codeword_length; i++)
	{
		tx_symbol[i][0] = -1 * (2 * codeword[i] - 1);
		tx_symbol[i][1]=0;
	}
}

void channel()
{
	//AWGN 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;
		}
	}
}

void demodulation()
{
	int i;
	double d1, d2;
	for (i = 0; i<codeword_length; i++)
	{
		d1 = (rx_symbol[i][0] - 1)*(rx_symbol[i][0] - 1) + rx_symbol[i][1] * rx_symbol[i][1];
		d2 = (rx_symbol[i][0] + 1)*(rx_symbol[i][0] + 1) + rx_symbol[i][1] * rx_symbol[i][1];
		if (d1<d2)
			re_codeword[i] = 0;
		else
			re_codeword[i] = 1;
	}
}

void decoder()
{

}


// ----------------------------------------------------------------
// BCJR (MAP) Decoder Implementation
// ----------------------------------------------------------------
void bcjr_decoder()
{
    int t, state, input;
    
    // 1. 初始化 (Initialization)
    // 清空 alpha 和 beta 数组
    for (t = 0; t <= message_length; t++) {
        for (state = 0; state < 4; state++) {
            alpha[t][state] = 0.0;
            beta[t][state] = 0.0;
        }
    }

    // 初始状态概率：编码器从状态0开始，所以 alpha[0][0] = 1，其他为0
    alpha[0][0] = 1.0;

    // 结束状态概率：由于有尾比特(Zero Padding)，编码器必然回到状态0
    // 所以 beta[message_length][0] = 1，其他为0
    beta[message_length][0] = 1.0; 

    // 2. 前向递归 (Forward Recursion) - 计算 Alpha
    for (t = 0; t < message_length; t++) {
        double sum_alpha = 0.0; // 用于归一化
        
        for (int cur_s = 0; cur_s < 4; cur_s++) {
            // 如果当前状态概率已经为0，则无需计算其分支
            if (alpha[t][cur_s] == 0.0) continue; 

            for (input = 0; input <= 1; input++) {
                // 查表获取下一状态和输出比特
                int col_next = (input == 0) ? 0 : 2;
                int col_out  = (input == 0) ? 1 : 3;

                int next_s = state_table[cur_s][col_next];
                int output_sym = state_table[cur_s][col_out];

                // 解析理论输出的两个比特
                int c0 = (output_sym >> 1) & 1; // 高位
                int c1 = output_sym & 1;        // 低位

                // --- 计算分支度量 Gamma ---
                // Gamma = P(y|x) * P(u)
                // P(u) 假定等概为 0.5 (可以忽略常数，因为会归一化)
                // P(y|x) 对于高斯信道正比于 exp(-d^2 / 2sigma^2)
                
                // 将理论比特映射为 BPSK 符号: 0->1.0, 1->-1.0 (与 modulation 函数一致)
                double tx0 = (c0 == 0) ? 1.0 : -1.0;
                double tx1 = (c1 == 0) ? 1.0 : -1.0;

                // 获取当前时刻接收到的软信息 (rx_symbol)
                double r0 = rx_symbol[2 * t][0];
                double r1 = rx_symbol[2 * t + 1][0];

                // 计算欧氏距离的平方 (Euclidean Distance Squared)
                double dist_sq = (r0 - tx0)*(r0 - tx0) + (r1 - tx1)*(r1 - tx1);

                // 计算 Gamma
                double gamma = exp(-dist_sq / (2.0 * sgm * sgm));

                // Alpha 更新公式: alpha[t+1][next] += alpha[t][current] * gamma
                alpha[t + 1][next_s] += alpha[t][cur_s] * gamma;
            }
        }

        // --- 归一化 Alpha ---
        // 防止浮点数下溢 (Underflow)，将当前时刻所有状态的概率和设为1
        for (state = 0; state < 4; state++) sum_alpha += alpha[t + 1][state];
        if (sum_alpha > 0) {
            for (state = 0; state < 4; state++) alpha[t + 1][state] /= sum_alpha;
        }
    }

    // 3. 后向递归 (Backward Recursion) - 计算 Beta
    for (t = message_length - 1; t >= 0; t--) {
        double sum_beta = 0.0; // 用于归一化

        for (int cur_s = 0; cur_s < 4; cur_s++) {
            for (input = 0; input <= 1; input++) {
                // 这里我们要看：从 cur_s 输入 input 会去哪？(next_s)
                // 然后利用 beta[t+1][next_s] 来更新 beta[t][cur_s]
                
                int col_next = (input == 0) ? 0 : 2;
                int col_out  = (input == 0) ? 1 : 3;

                int next_s = state_table[cur_s][col_next];
                int output_sym = state_table[cur_s][col_out];

                // 计算 Gamma (同前向递归)
                int c0 = (output_sym >> 1) & 1;
                int c1 = output_sym & 1;
                double tx0 = (c0 == 0) ? 1.0 : -1.0;
                double tx1 = (c1 == 0) ? 1.0 : -1.0;
                double r0 = rx_symbol[2 * t][0];
                double r1 = rx_symbol[2 * t + 1][0];
                double dist_sq = (r0 - tx0)*(r0 - tx0) + (r1 - tx1)*(r1 - tx1);
                double gamma = exp(-dist_sq / (2.0 * sgm * sgm));

                // Beta 更新公式: beta[t][current] += beta[t+1][next] * gamma
                beta[t][cur_s] += beta[t + 1][next_s] * gamma;
            }
        }

        // --- 归一化 Beta ---
        for (state = 0; state < 4; state++) sum_beta += beta[t][state];
        if (sum_beta > 0) {
            for (state = 0; state < 4; state++) beta[t][state] /= sum_beta;
        }
    }

    // 4. 计算对数似然比 (LLR) 并判决
    // P(u_t = 0 | y) vs P(u_t = 1 | y)
    for (t = 0; t < message_length; t++) {
        double prob0 = 0.0; // 输入为0的联合概率和
        double prob1 = 0.0; // 输入为1的联合概率和

        for (int cur_s = 0; cur_s < 4; cur_s++) {
            if (alpha[t][cur_s] == 0.0) continue;

            for (input = 0; input <= 1; input++) {
                int col_next = (input == 0) ? 0 : 2;
                int col_out  = (input == 0) ? 1 : 3;
                int next_s = state_table[cur_s][col_next];
                int output_sym = state_table[cur_s][col_out];

                // 重新计算 Gamma
                int c0 = (output_sym >> 1) & 1;
                int c1 = output_sym & 1;
                double tx0 = (c0 == 0) ? 1.0 : -1.0;
                double tx1 = (c1 == 0) ? 1.0 : -1.0;
                double r0 = rx_symbol[2 * t][0];
                double r1 = rx_symbol[2 * t + 1][0];
                double dist_sq = (r0 - tx0)*(r0 - tx0) + (r1 - tx1)*(r1 - tx1);
                double gamma = exp(-dist_sq / (2.0 * sgm * sgm));

                // 这一条路径的概率度量 = alpha * gamma * beta
                double metric = alpha[t][cur_s] * gamma * beta[t + 1][next_s];

                if (input == 0) prob0 += metric;
                else            prob1 += metric;
            }
        }

        // MAP 判决
        if (prob1 > prob0) de_message[t] = 1;
        else               de_message[t] = 0;
    }
}