import math
import numpy as np
import matplotlib.pyplot as plt
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(occurrences.astype(int)) * idx))
indexes = np.random.permutation(samples.shape[0])
return samples[indexes]
def xesmd(apps, idx):
'''
Estimates symbolwise equivocation from reference symbols indices and a posteriori probabilities.
'''
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)
n = 100_000
SNR_dBs = np.arange(5,22)
np.arange(0,2.5,0.01).shape
M = 8
plots_dict = np.empty(shape=(17,25000))
for idx, shaping in enumerate(np.arange(0,2.5,0.0001)):
alphabet = np.arange(-7,8,2)
alphabet = alphabet / np.sqrt(np.mean(alphabet**2))
scaling = (alphabet[1]-alphabet[0])/2
# scaling = 0.21821789023599236
# print('Scaling: ', scaling)
denominator = np.sum(np.exp(-shaping*np.square(alphabet)))
probs = 1/denominator * np.exp(-shaping*np.square(alphabet))
power = np.sum(np.square(np.arange(-7,8,2)*(scaling))*probs)
# print('Power: ', power)
indices = sampler(probs, n).astype(int)
norm_factor = np.sqrt(np.sum(np.square(alphabet) * probs))
alphabet = alphabet / norm_factor
symbols = alphabet[indices]
scaling = (alphabet[1]-alphabet[0])/2
# print('Scaling: ', scaling)
H = np.sum(-np.log2(probs)*probs)
# print('Entropy: ',H)
power = np.mean(symbols**2)
# print('Power: ', power)
# plt.rcParams['figure.figsize'] = [4, 4]
# plt.hist(symbols, bins=100)
# plt.show()
for jdx, snrdB in enumerate(SNR_dBs):
sigma2 = 1/(10**(snrdB/10))
sigma2 = sigma2
y = AWGN_channel(symbols, sigma2)
apps = AWGNdemapper(y, alphabet, sigma2)
xe = xesmd(apps, indices)
plots_dict[jdx][idx] = (2*(H - (xe) / np.log(2)))
print(np.max(plots_dict, 1).tolist())
mb = np.max(plots_dict, 1)
mb = [2.0193916630846918, 2.2654022278277393, 2.5231306887622242, 2.7943201472250596, 3.0806359602235873, 3.3809640408051758, 3.69383842699387, 4.010358520455199, 4.318829495451576, 4.61704902872525, 4.895880050688773, 5.151007989533582, 5.376996748960784, 5.569021543421744, 5.722730806499707, 5.841398311333901, 5.9200053405287045]
from numpy import genfromtxt
my_data = genfromtxt('/home/ddeandres/Projects/internship_pcs/documentation/figs/mb_8_ask.csv', delimiter=',')
plt.plot(my_data[:,0], my_data[:,1])
# Plot
plt.rcParams['figure.figsize'] = [8, 6]
for value in plots_dict:
# print(key, value)
plt.plot(SNR_dBs, mb)
plt.plot(my_data[:,0], 2*my_data[:,1])
plt.plot(SNR_dBs, np.log2(1+10**(SNR_dBs/10)), color='black', label='Capacity')
plt.grid()
解释代码,并进行注释,请注意,可以对代码进行规范化,但不要修改原本的结构。另外,将原本的信噪比SNR_dBs = np.arange(5,22)修改为仅包含5,10,15,20.增加星座点可视化,可参考下列代码:
3. 星座图可视化
plt.figure(figsize=(4, 4))
# 计算二维网格点(实部+虚部,本例为实数调制)
alph = alphabet_norm.detach().numpy()
a = np.array([(d, c) for c in np.flip(alph) for d in alph])
# 计算联合概率
joint_probs = (probs.reshape(-1, 1) * probs.reshape(1, -1)).flatten().detach().numpy()
# 绘制散点图(点大小表示概率)
plt.scatter(a[:, 0], a[:, 1], s=joint_probs * 2000, alpha=0.6)
plt.title(“Constellation Diagram”)
plt.xlabel(“In-phase”)
plt.ylabel(“Quadrature”)
# 保存为TikZ格式(用于LaTeX文档)
tikzplotlib_save(f"aref_pcs_{SNR_dB}dB.tex")
plt.show()