<think>我们首先对代码进行规范化,并添加注释。同时,我们将按照要求修改SNR_dBs为仅包含5,10,15,20,并增加星座点可视化。
注意:原代码中有些地方有语法错误(如字符串使用中文引号,缩进不一致等),我们将修正这些错误。
步骤:
1. 修正语法错误,调整代码结构。
2. 添加注释。
3. 修改SNR_dBs为[5,10,15,20]。
4. 为16QAM、64QAM(均匀和概率整形)、256QAM添加星座图可视化。
5. 注意:原代码中有些函数(如`my_bisection`)在后续未使用,但根据要求保留。
由于代码较长,我们将分部分进行。
注意:原代码中`tikzplotlib`用于生成LaTeX的TikZ代码,我们保留,但如果没有安装,可以注释掉或者用普通保存。
首先,我们导入必要的库,并定义函数。</think>### 解释与注释代码
以下是对代码的详细解释和注释,同时进行了规范化处理(修正了语法错误、调整了缩进、优化了变量命名),但保留了原始结构:
```python
import math
import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
import torch.optim as optim
# import tikzplotlib # 注释掉,因为可能未安装
def my_bisection(f, a, b, P, tol):
"""
二分法求根函数
f: 目标函数
a, b: 搜索区间
P: 目标值
tol: 容差
"""
# 检查区间是否包含根
if (f(a) > P and f(b) > P) or (f(a) < P and f(b) < P):
raise Exception("The scalars a and b do not bound P")
# 计算中点
m = (a + b) / 2
# 停止条件
if np.abs(f(m) - P) < tol:
return m
# 递归搜索
elif f(m) > P and f(m) < f(a):
return my_bisection(f, m, b, P, tol)
elif f(m) < P and f(m) > f(b):
return my_bisection(f, a, m, P, tol)
def AWGN_channel(x, sigma2):
"""添加高斯白噪声"""
noise = np.sqrt(sigma2) * np.random.randn(x.size)
return x + noise
def AWGNdemapper(y, const, varN):
"""计算后验概率"""
apps = np.exp(-np.abs(np.transpose([y]) - const)**2 / (2 * varN))
return apps / np.transpose([np.sum(apps, 1)])
def sampler(prob, n):
"""根据概率分布生成样本"""
samples = np.empty(0)
for idx, p in enumerate(prob):
occurrences = np.round(n * p)
samples = np.concatenate((samples, np.ones(int(occurrences)) * idx))
indexes = np.random.permutation(samples.shape[0])
return samples[indexes]
def xesmd(apps, idx):
"""计算符号平均互信息"""
eq = -np.log(np.take_along_axis(apps, idx[:, None], axis=1) / np.transpose([np.sum(apps, 1)]))
eq[eq == np.inf] = 1000 # 处理无穷大值
return np.mean(eq)
def power_for_shaping(shaping):
"""计算给定整形参数下的功率"""
alphabet = np.arange(-7, 8, 2)
scaling = (alphabet[1] - alphabet[0]) / 2
# 计算概率分布
denominator = np.sum(np.exp(-shaping * np.square(alphabet)))
probs = np.exp(-shaping * np.square(alphabet)) / denominator
# 计算功率
power = np.sum(np.square(np.arange(-7, 8, 2) * scaling) * probs)
return power
# ================= 主程序 =================
if __name__ == "__main__":
# 设置参数
n = 100_000 # 样本数
SNR_dBs = np.array([5, 10, 15, 20]) # 修改为指定SNR点
# ===== 16QAM仿真 =====
M = 4
alphabet = np.arange(-3, 4, 2) # 星座点 [-3, -1, 1, 3]
alphabet = alphabet / np.sqrt(np.mean(alphabet**2)) # 功率归一化
indices = np.random.choice(np.arange(M), n)
symbols = alphabet[indices]
# 星座图可视化
plt.figure(figsize=(4, 4))
plt.scatter(symbols, np.zeros_like(symbols), alpha=0.3)
plt.title("16QAM Constellation")
plt.xlabel("In-phase")
plt.grid(True)
plt.show()
mi_16 = []
for snrdB in SNR_dBs:
sigma2 = 1 / (10**(snrdB/10)) # 噪声方差
y = AWGN_channel(symbols, sigma2) # 通过AWGN信道
apps = AWGNdemapper(y, alphabet, sigma2) # 计算后验概率
xe = xesmd(apps, indices) # 计算互信息
mi_16.append(2 * (2 - xe / np.log(2))) # 16QAM每符号2比特
# ===== 64QAM均匀分布 =====
M = 8
alphabet = np.arange(-7, 8, 2) # [-7, -5, -3, -1, 1, 3, 5, 7]
alphabet = alphabet / np.sqrt(np.mean(alphabet**2)) # 功率归一化
indices = np.random.choice(np.arange(M), n)
symbols = alphabet[indices]
# 星座图可视化
plt.figure(figsize=(4, 4))
plt.scatter(symbols, np.zeros_like(symbols), alpha=0.3)
plt.title("64QAM Uniform Constellation")
plt.xlabel("In-phase")
plt.grid(True)
plt.show()
mi_64 = []
for snrdB in SNR_dBs:
sigma2 = 1 / (10**(snrdB/10))
y = AWGN_channel(symbols, sigma2)
apps = AWGNdemapper(y, alphabet, sigma2)
xe = xesmd(apps, indices)
mi_64.append(2 * (3 - xe / np.log(2))) # 64QAM每符号3比特
# ===== 概率整形64QAM =====
shaping = 2.5
alphabet = np.arange(-7, 8, 2)
scaling = 0.21821789023599236 # 固定缩放因子
# 计算概率分布
denominator = np.sum(np.exp(-shaping * np.square(alphabet)))
probs = np.exp(-shaping * np.square(alphabet)) / denominator
# 生成样本
indices = sampler(probs, n).astype(int)
# 功率归一化
norm_factor = np.sqrt(np.sum(np.square(alphabet) * probs))
alphabet = alphabet / norm_factor
symbols = alphabet[indices]
# 计算熵和功率
H = np.sum(-np.log2(probs) * probs)
power = np.mean(symbols**2)
# 星座图可视化(概率整形)
plt.figure(figsize=(4, 4))
plt.hist(symbols, bins=100, density=True)
plt.title("64QAM Probabilistic Shaping")
plt.xlabel("Symbol Value")
plt.ylabel("Probability Density")
plt.show()
mi_mb_64 = []
for snrdB in SNR_dBs:
sigma2 = 1 / (10**(snrdB/10))
y = AWGN_channel(symbols, sigma2)
apps = AWGNdemapper(y, alphabet, sigma2)
xe = xesmd(apps, indices)
mi_mb_64.append(H - xe / np.log(2)) # 使用实际熵
# ===== 256QAM仿真 =====
M = 16
alphabet = np.arange(-15, 16, 2) # 16个点
alphabet = alphabet / np.sqrt(np.mean(alphabet**2)) # 功率归一化
indices = np.random.choice(np.arange(M), n)
symbols = alphabet[indices]
# 星座图可视化
plt.figure(figsize=(4, 4))
plt.scatter(symbols, np.zeros_like(symbols), alpha=0.3)
plt.title("256QAM Constellation")
plt.xlabel("In-phase")
plt.grid(True)
plt.show()
mi_256 = []
for snrdB in SNR_dBs:
sigma2 = 1 / (10**(snrdB/10))
y = AWGN_channel(symbols, sigma2)
apps = AWGNdemapper(y, alphabet, sigma2)
xe = xesmd(apps, indices)
mi_256.append(2 * (4 - xe / np.log(2))) # 256QAM每符号4比特
# ===== 结果绘图 =====
plt.figure(figsize=(8, 6))
plt.plot(SNR_dBs, mi_16, label='16QAM')
plt.plot(SNR_dBs, mi_64, label='64QAM Uniform')
plt.plot(SNR_dBs, mi_mb_64, label='64QAM Prob Shaped')
plt.plot(SNR_dBs, mi_256, label='256QAM')
plt.plot(SNR_dBs, np.log2(1 + 10**(SNR_dBs/10)),
color='black', label='AWGN Capacity')
plt.legend()
plt.grid()
plt.title('AIR of QAM Constellations')
plt.ylabel('Bits per Channel Use')
plt.xlabel('SNR (dB)')
plt.savefig("QAM_MI.png")
plt.show()
# 容量间隙分析
plt.figure(figsize=(8, 6))
plt.plot(SNR_dBs, np.log2(1 + 10**(SNR_dBs/10)),
color='black', label=r'$C(P/\sigma^2)$')
plt.plot(SNR_dBs, np.log2(1 + 10**(SNR_dBs/10)) - 0.5 * np.log2((np.pi * np.e)/6),
linestyle='dashed', color='C1',
label=r'$C(P/\sigma^2) - \frac{1}{2}\log_2\frac{\pi e}{6}$')
plt.grid()
plt.ylabel('Bits per Channel Use')
plt.xlabel('SNR (dB)')
plt.xlim([5, 20])
plt.ylim([1, 8])
plt.title('AWGN Channel Capacity Gap')
plt.legend()
plt.savefig("SNR_GAP.png")
plt.show()
# 测试二分法
f = lambda x: -np.log10(x)
P = 0.4
p1 = my_bisection(f, 0.2, 0.8, P, 0.1)
print("P1 =", p1)
p01 = my_bisection(f, 0.2, 0.8, P, 0.01)
print("P01 =", p01)
print("f(r1) =", f(p1))
print("f(r01) =", f(p01))
```
### 主要修改说明:
1. **SNR点修改**:
```python
SNR_dBs = np.array([5, 10, 15, 20]) # 仅保留指定点
```
2. **星座图可视化**:
为每种调制方式添加了星座图可视化:
```python
# 16QAM星座图
plt.figure(figsize=(4, 4))
plt.scatter(symbols, np.zeros_like(symbols), alpha=0.3)
plt.title("16QAM Constellation")
plt.xlabel("In-phase")
plt.grid(True)
plt.show()
```
3. **概率整形可视化**:
使用直方图展示概率分布:
```python
plt.hist(symbols, bins=100, density=True)
plt.title("64QAM Probabilistic Shaping")
```
4. **错误修复**:
- 修正了`sampler`函数中的整数转换问题
- 修复了字符串中的中文引号问题
- 规范了缩进和空格使用
- 添加了完整的函数文档字符串
5. **图像保存**:
使用`plt.savefig`替代`tikzplotlib`保存为PNG格式
6. **代码结构优化**:
- 将主程序放在`if __name__ == "__main__":`中
- 添加清晰的章节注释
- 使用更规范的变量命名
### 代码功能总结:
1. 实现了AWGN信道中的QAM调制仿真
2. 比较了16QAM、64QAM(均匀和概率整形)、256QAM的性能
3. 计算并绘制了互信息随SNR变化曲线
4. 分析了AWGN信道容量与QAM性能的差距
5. 包含星座图可视化展示调制方案