Viterbi_Hard_K4_1513.c

/***************************************************
Channel Coding Course Work: Convolutional Codes (Modified)
Specs: (2, 1, 4) Convolutional Code
Generators: (15, 13) octal
Constraint Length: K=4
***************************************************/

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

// --- 修改 1: 状态数改为 8，表大小改为 [8][4] ---
int state_table[8][4]; 

// --- 修改 2: 状态数变量改为 8 (K=4, m=3, 2^3=8) ---
int state_num = 8; 

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

// --- 修改 3: 回溯路径数组第二维改为 8 ---
int traceback_path[message_length][8];

// 定义无穷大 (用于 Viterbi 解码)
#define INF 999999 

void statetable();
void encoder();
void modulation();
void demodulation();
void channel();
void decoder();

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文件指针 ---
    FILE *fp = NULL; 
    // 打开 csv 文件
    fp = fopen("hard_viterbi_K4_15_13.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--- Hard Decision Viterbi Simulation (K=4, g=15,13) ---\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
			}
			for (i = message_length - state_num; i<message_length; i++)
			{
				message[i] = 0;//信息后面添加与编码器状态数相等个数的0
			}

			//convolutional encoder
			encoder();

			//BPSK modulation
			modulation();

			//AWGN channel
			channel();

			//BPSK demodulation
			demodulation();

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

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

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

    return 0;
}

// --- 修改 4: 重写 statetable 适配 (15, 13) ---
void statetable()
{
    int s0, s1, s2, u;
    int next_state;
    int c1, c2;

    // 遍历所有可能的当前状态 (0 ~ 7)
    // 状态 s 对应寄存器位 (s0, s1, s2)，其中 s0 是最近存入的位
    for (int s = 0; s < 8; s++) 
    {
        // 提取二进制位
        s2 = s & 1;          // 最老的位 (D^3)
        s1 = (s >> 1) & 1;   // 中间位 (D^2)
        s0 = (s >> 2) & 1;   // 最新位 (D^1)

        // 遍历所有可能的输入 (0 或 1)
        for (u = 0; u <= 1; u++) 
        {
            // g1 = (15)8 = 1101 -> 1*u + 1*s0 + 0*s1 + 1*s2
            c1 = u ^ s0 ^ s2;
            
            // g2 = (13)8 = 1011 -> 1*u + 0*s0 + 1*s1 + 1*s2
            c2 = u ^ s1 ^ s2;
            
            // 下一个状态: 输入 u 移入最高位，其他位右移
            // next_state = (u, s0, s1)
            next_state = (u << 2) | (s >> 1);
            
            // 输出符号 (c1 高位, c2 低位)
            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; // 初始状态 0
    int input_bit;
    int output_symbol;     

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

        // 查表
        int col_next_state = (input_bit == 0) ? 0 : 2;
        int col_output     = (input_bit == 0) ? 1 : 3;

        output_symbol = state_table[current_state][col_output];
        int next_state = state_table[current_state][col_next_state];

        codeword[2 * i] = (output_symbol >> 1) & 1;
        codeword[2 * i + 1] = output_symbol & 1;

        current_state = next_state;
    }
}

void 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()
{
    int t, state, input;
    // --- 修改 5: 内部度量数组大小必须至少为 8 ---
    int path_metric[8];     
    int new_path_metric[8]; 

    // 初始化
    path_metric[0] = 0;
    for (state = 1; state < 8; state++) path_metric[state] = INF; // 循环范围改为 8

    // 前向递推
    for (t = 0; t < message_length; t++)
    {
        int r0 = re_codeword[2 * t];
        int r1 = re_codeword[2 * t + 1];

        // 初始化下一时刻度量
        for (state = 0; state < 8; state++) new_path_metric[state] = INF; // 范围改为 8

        // 遍历当前时刻的所有可能状态 (0 ~ 7)
        for (int curr_state = 0; curr_state < 8; curr_state++) // 范围改为 8
        {
            if (path_metric[curr_state] >= INF) continue;

            for (input = 0; input <= 1; input++)
            {
                int col_next = (input == 0) ? 0 : 2;
                int col_out  = (input == 0) ? 1 : 3;

                int next_state = state_table[curr_state][col_next];
                int output_sym = state_table[curr_state][col_out];

                int c0 = (output_sym >> 1) & 1; 
                int c1 = output_sym & 1;        

                int branch_metric = 0;
                if (r0 != c0) branch_metric++;
                if (r1 != c1) branch_metric++;

                int total_metric = path_metric[curr_state] + branch_metric;

                if (total_metric < new_path_metric[next_state])
                {
                    new_path_metric[next_state] = total_metric;
                    // 记录回溯路径
                    traceback_path[t][next_state] = curr_state;
                }
            }
        }

        // 更新度量
        for (state = 0; state < 8; state++) // 范围改为 8
        {
            path_metric[state] = new_path_metric[state];
        }
    }

    // 回溯
    int curr_state = 0; 
    for (t = message_length - 1; t >= 0; t--)
    {
        int prev_state = traceback_path[t][curr_state];

        // --- 修改 6: 从状态中提取输入位 ---
        // 现在的状态是 (input, s0, s1)，所以输入位在最高位 (bit 2)
        int decoded_bit = (curr_state >> 2) & 1;
        de_message[t] = decoded_bit;

        curr_state = prev_state;
    }
}