#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
// BCH(7,4)编码参数 (t=1, 本原多项式x³+x+1)
#define N_BCH 7
#define K_BCH 4
#define T_BCH 1 // 纠错能力
#define GF_M 3 // 有限域GF(2^3)
// GF(2^3)有限域定义（本原多项式x³+x+1）
// 指数表: alpha^i -> 多项式表示 (i=0~6)
const int alpha_to[8] = {1, 2, 4, 3, 6, 7, 5, 0}; // alpha^0=1, alpha^1=2, ..., alpha^6=5
// 对数表: 多项式表示 -> alpha^i (值0对应索引7)
const int index_of[8] = {7, 0, 1, 3, 2, 6, 4, 5}; // 1=alpha^0, 2=alpha^1, ..., 5=alpha^6
// 有限域加法（异或）
static inline int gf_add(int a, int b) {
return a ^ b;
}
// 有限域乘法（基于指数表和对数表）
static inline int gf_mul(int a, int b) {
if (a == 0 || b == 0) return 0;
// alpha^i * alpha^j = alpha^((i+j) mod 7)
return alpha_to[(index_of[a] + index_of[b]) % 7];
}
// 有限域乘法逆元（a^(-1) = alpha^(-i) = alpha^(7-i mod 7)）
static inline int gf_inv(int a) {
if (a == 0) return 0; // 0没有逆元
return alpha_to[(7 - index_of[a]) % 7];
}
// 高斯随机数生成器（Box-Muller变换）
double gaussrand(double mean, double stddev) {
static double V1, V2, S;
static int phase = 0;
double X;
if (phase == 0) {
do {
double U1 = (double)rand() / RAND_MAX;
double U2 = (double)rand() / RAND_MAX;
V1 = 2.0 * U1 - 1.0;
V2 = 2.0 * U2 - 1.0;
S = V1 * V1 + V2 * V2;
} while (S >= 1.0 || S == 0.0);
X = V1 * sqrt(-2.0 * log(S) / S);
} else {
X = V2 * sqrt(-2.0 * log(S) / S);
}
phase = 1 - phase;
return mean + X * stddev;
}
// BCH(7,4)编码函数（保持不变）
void bch_encode(const int* input, int* output, int num_blocks) {
int G[K_BCH][N_BCH] = {
{1, 0, 0, 0, 1, 1, 0},
{0, 1, 0, 0, 1, 0, 1},
{0, 0, 1, 0, 0, 1, 1},
{0, 0, 0, 1, 1, 1, 1}
};
for (int i = 0; i < num_blocks; i++) {
for (int j = 0; j < N_BCH; j++) {
output[i * N_BCH + j] = 0;
for (int k = 0; k < K_BCH; k++) {
output[i * N_BCH + j] ^= input[i * K_BCH + k] * G[k][j];
}
}
}
}
// 计算校正子（关键修复：基于GF(2^3)的正确计算）
void calculate_syndromes(const int* received, int* syndromes) {
// BCH(7,4)纠错t=1，需要2t=2个校正子（S1, S3）
// 校正子定义：S_j = r(alpha^j) = sum_{i=0}^{6} received[i] * (alpha^j)^i
syndromes[0] = 0; // S1 = r(alpha^1)
syndromes[1] = 0; // S3 = r(alpha^3)
for (int i = 0; i < N_BCH; i++) { // i: 接收位位置（对应x^i）
if (received[i]) {
// 计算(alpha^1)^i = alpha^i，累加到S1
int alpha_pow_i = alpha_to[i % 7]; // alpha^i
syndromes[0] = gf_add(syndromes[0], alpha_pow_i);
// 计算(alpha^3)^i = alpha^(3i)，累加到S3
int alpha_pow_3i = alpha_to[(3 * i) % 7]; // alpha^(3i)
syndromes[1] = gf_add(syndromes[1], alpha_pow_3i);
}
}
}
// Berlekamp-Massey算法（关键修复：多项式更新和逆元计算）
void berlekamp_massey(const int* syndromes, int* sigma, int* deg) {
int L = 0; // 错误位置多项式当前阶数
int m = 1; // 上次更新后的步数
int b = 1; // 上次更新的偏差值
int sigma_old[3] = {1, 0, 0}; // 前一次的多项式
// 初始化错误位置多项式：sigma(x) = 1
sigma[0] = 1;
sigma[1] = 0;
sigma[2] = 0;
// 迭代处理2t个校正子
for (int n = 0; n < 2 * T_BCH; n++) {
// 计算偏差d = S_{n+1} + sum_{i=1}^L sigma_i * S_{n+1-i}
int d = syndromes[n];
for (int i = 1; i <= L; i++) {
if (n - i >= 0) {
d = gf_add(d, gf_mul(sigma[i], syndromes[n - i]));
}
}
if (d == 0) { // 偏差为0，仅更新步数
m++;
continue;
}
// 保存当前多项式用于后续更新
int t[3];
memcpy(t, sigma, 3 * sizeof(int));
// 计算更新量：d * b^{-1} * x^m * sigma_old(x)
int b_inv = gf_inv(b); // b的逆元
int c = gf_mul(d, b_inv); // 缩放因子：d / b
for (int i = 0; i <= L; i++) {
if (i + m < 3) { // 防止越界
sigma[i + m] = gf_add(sigma[i + m], gf_mul(c, sigma_old[i]));
}
}
// 若当前阶数不足，更新阶数
if (2 * L <= n) {
L = n + 1 - L;
memcpy(sigma_old, t, 3 * sizeof(int)); // 保存当前多项式
b = d; // 更新偏差基准
m = 1; // 重置步数
} else {
m++;
}
}
*deg = L; // 返回多项式阶数
}
// 钱搜索算法（关键修复：根的验证和错误位置映射）
int chien_search(const int* sigma, int deg, int* error_positions) {
int num_errors = 0;
// 错误位置多项式：sigma(x) = 1 + sigma_1 x + sigma_2 x^2
// 根为alpha^{-k}，对应错误位置k（即x^k项出错）
for (int k = 0; k < N_BCH; k++) { // k: 可能的错误位置（0~6）
// 计算sigma(alpha^{-k}) = 1 + sigma_1*alpha^{-k} + sigma_2*alpha^{-2k}
int x = alpha_to[(7 - k) % 7]; // alpha^{-k} = alpha^(7-k mod7)
int val = sigma[0]; // 常数项1
if (deg >= 1) {
val = gf_add(val, gf_mul(sigma[1], x)); // + sigma_1 * x
}
if (deg >= 2) {
int x2 = gf_mul(x, x); // x^2 = alpha^{-2k}
val = gf_add(val, gf_mul(sigma[2], x2)); // + sigma_2 * x^2
}
if (val == 0) { // alpha^{-k}是根，错误位置为k
error_positions[num_errors++] = k;
if (num_errors >= T_BCH) break; // 达到最大纠错能力
}
}
return num_errors;
}
// BCH解码函数（整合校正子、BM算法和钱搜索）
void bch_decode(const int* input, int* output, int num_blocks) {
for (int i = 0; i < num_blocks; i++) {
int received[N_BCH];
memcpy(received, &input[i * N_BCH], N_BCH * sizeof(int));
// 1. 计算校正子
int syndromes[2 * T_BCH];
calculate_syndromes(received, syndromes);
// 2. 若校正子全为0，无错误
if (syndromes[0] == 0 && syndromes[1] == 0) {
for (int j = 0; j < K_BCH; j++) {
output[i * K_BCH + j] = received[j];
}
continue;
}
// 3. BM算法求错误位置多项式
int sigma[3]; // 错误位置多项式（最高2阶）
int deg; // 多项式阶数
berlekamp_massey(syndromes, sigma, °);
// 4. 钱搜索找错误位置
int error_positions[T_BCH];
int num_errors = chien_search(sigma, deg, error_positions);
// 5. 纠正错误
for (int j = 0; j < num_errors; j++) {
int pos = error_positions[j];
if (pos >= 0 && pos < N_BCH) {
received[pos] ^= 1; // 翻转错误位
}
}
// 6. 提取信息位（前4位）
for (int j = 0; j < K_BCH; j++) {
output[i * K_BCH + j] = received[j];
}
}
}
// 计算误码率
double calculate_ber(const int* original, const int* received, int length) {
if (length <= 0) return 0.0;
int errors = 0;
for (int i = 0; i < length; i++) {
if (original[i] != received[i]) errors++;
}
return (double)errors / length;
}
int main() {
srand((unsigned int)time(NULL));
// 参数初始化（确保数据长度为K_BCH的整数倍）
const int base_length = 40000000;
const int num_blocks = base_length / K_BCH;
const int n = num_blocks * K_BCH;
const double SNR_dB_start = 0.0;
const double SNR_dB_end = 13.0;
const double SNR_dB_step = 1.0;
const int num_snr = (int)((SNR_dB_end - SNR_dB_start) / SNR_dB_step) + 1;
// 分配BER结果数组
double* BER_encoded = (double*)malloc(num_snr * sizeof(double));
double* BER_theory = (double*)malloc(num_snr * sizeof(double));
double* BER_test = (double*)malloc(num_snr * sizeof(double));
if (!BER_encoded || !BER_theory || !BER_test) {
fprintf(stderr, "内存分配失败\n");
exit(1);
}
// 生成随机数据
int* data = (int*)malloc(n * sizeof(int));
if (!data) {
fprintf(stderr, "数据内存分配失败\n");
exit(1);
}
for (int i = 0; i < n; i++) {
data[i] = rand() % 2;
}
// BCH编码
int* encoded_bits = (int*)malloc(num_blocks * N_BCH * sizeof(int));
if (!encoded_bits) {
fprintf(stderr, "编码内存分配失败\n");
exit(1);
}
bch_encode(data, encoded_bits, num_blocks);
const int encoded_length = num_blocks * N_BCH;
// BPSK调制
double* BPSK_encoded = (double*)malloc(encoded_length * sizeof(double));
double* BPSK_test = (double*)malloc(n * sizeof(double));
if (!BPSK_encoded || !BPSK_test) {
fprintf(stderr, "调制内存分配失败\n");
exit(1);
}
for (int i = 0; i < encoded_length; i++) {
BPSK_encoded[i] = encoded_bits[i] ? 1.0 : -1.0;
}
for (int i = 0; i < n; i++) {
BPSK_test[i] = data[i] ? 1.0 : -1.0;
}
printf("开始仿真...数据长度: %d bits, SNR范围: %.1f dB 到 %.1f dB\n",
n, SNR_dB_start, SNR_dB_end);
printf("编码方案: BCH(7,4), 纠错能力: %d bit\n", T_BCH);
printf("解码算法: Berlekamp-Massey算法\n");
// 主仿真循环
for (int z = 0; z < num_snr; z++) {
double snr_db = SNR_dB_start + z * SNR_dB_step;
double snr_lin = pow(10.0, snr_db / 10.0);
// 噪声标准差（确保Eb/N0公平比较）
double noise_std = 1.0 / sqrt(2.0 * snr_lin);
double noise_std_encoded = 1.0 / sqrt(2.0 * snr_lin * (double)K_BCH / N_BCH);
// 编码系统处理
double* noisy_encoded = (double*)malloc(encoded_length * sizeof(double));
int* demod_encoded = (int*)malloc(encoded_length * sizeof(int));
int* decoded_bits = (int*)malloc(n * sizeof(int));
if (!noisy_encoded || !demod_encoded || !decoded_bits) {
fprintf(stderr, "SNR=%.1f dB内存分配失败\n", snr_db);
exit(1);
}
// 添加噪声并解调
for (int i = 0; i < encoded_length; i++) {
noisy_encoded[i] = BPSK_encoded[i] + gaussrand(0.0, noise_std_encoded);
}
for (int i = 0; i < encoded_length; i++) {
demod_encoded[i] = (noisy_encoded[i] > 0.0) ? 1 : 0;
}
// BM算法译码
bch_decode(demod_encoded, decoded_bits, num_blocks);
BER_encoded[z] = calculate_ber(data, decoded_bits, n);
// 释放内存
free(noisy_encoded);
free(demod_encoded);
free(decoded_bits);
// 未编码系统处理
double* noisy_test = (double*)malloc(n * sizeof(double));
int* demod_test = (int*)malloc(n * sizeof(int));
if (!noisy_test || !demod_test) {
fprintf(stderr, "SNR=%.1f dB未编码内存分配失败\n", snr_db);
exit(1);
}
// 添加噪声并解调
for (int i = 0; i < n; i++) {
noisy_test[i] = BPSK_test[i] + gaussrand(0.0, noise_std);
}
for (int i = 0; i < n; i++) {
demod_test[i] = (noisy_test[i] > 0.0) ? 1 : 0;
}
BER_test[z] = calculate_ber(data, demod_test, n);
BER_theory[z] = 0.5 * erfc(sqrt(snr_lin));
// 释放内存
free(noisy_test);
free(demod_test);
// 打印进度
printf("SNR = %4.1f dB: 编码BER = %8.3e, 未编码BER = %8.3e, 理论BER = %8.3e\n",
snr_db, BER_encoded[z], BER_test[z], BER_theory[z]);
}
// 保存结果
FILE* fp = fopen("ber_results.txt", "w");
if (fp) {
fprintf(fp, "SNR_dB BER_theory BER_encoded BER_test\n");
for (int z = 0; z < num_snr; z++) {
double snr_db = SNR_dB_start + z * SNR_dB_step;
fprintf(fp, "%.1f %.6e %.6e %.6e\n",
snr_db, BER_theory[z], BER_encoded[z], BER_test[z]);
}
fclose(fp);
printf("\n结果已保存到 ber_results.txt\n");
// 分析编码增益
int gain_start = -1;
for (int z = 0; z < num_snr; z++) {
if (BER_encoded[z] < BER_test[z]) {
gain_start = z;
break;
}
}
if (gain_start != -1) {
double snr = SNR_dB_start + gain_start * SNR_dB_step;
printf("编码增益从 %.1f dB 开始: 编码BER = %.2e, 未编码BER = %.2e\n",
snr, BER_encoded[gain_start], BER_test[gain_start]);
}
} else {
printf("保存结果失败\n");
}
// 释放所有内存
free(BER_encoded);
free(BER_theory);
free(BER_test);
free(data);
free(encoded_bits);
free(BPSK_encoded);
free(BPSK_test);
return 0;
}
最新发布
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
