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):
# approximates a root, R, of f bounded
# by a and b to within tolerance
# | f(m) | < tol with m the midpoint
# between a and b Recursive implementation
# check if a and b bound a root 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") # get midpoint m = (a + b)/2 if np.abs(f(m)-P) < tol: # stopping condition, report m as root return m elif f(m) > P and f(m) < f(a): # case where m is an improvement on a. # Make recursive call with a = m return my_bisection(f, m, b, P, tol) elif f(m) < P and f(m) > f(b): # case where m is an improvement on b. # Make recursive call with b = m 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(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,25)
16 QAM
M = 4
alphabet = np.arange(-3,4,2)
alphabet = alphabet / np.sqrt(np.mean(alphabet**2))
indices = np.random.choice(np.arange(M), n)
symbols = alphabet[indices]
mi_16 = []
for snrdB in SNR_dBs:
sigma2 = 1/(10**(snrdB/10))
sigma2 = sigma2
y = AWGN_channel(symbols, sigma2)
apps = AWGNdemapper(y, alphabet, sigma2)
xe = xesmd(apps, indices)
mi_16.append(2*(2 - xe / np.log(2)))
64 QAM
M = 8
alphabet = np.arange(-7,8,2)
alphabet = alphabet / np.sqrt(np.mean(alphabet**2))
indices = np.random.choice(np.arange(M), n)
symbols = alphabet[indices]
print('Power: ',np.mean(np.square(symbols)))
mi_64 = []
for snrdB in SNR_dBs:
sigma2 = 1/(10**(snrdB/10))
sigma2 = sigma2
y = AWGN_channel(symbols, sigma2)
apps = AWGNdemapper(y, alphabet, sigma2)
xe = xesmd(apps, indices)
mi_64.append(2*(3 - (xe) / np.log(2)))
def power_for_shaping(shaping):
alphabet = np.arange(-7,8,2)
# alphabet = alphabet / np.sqrt(np.mean(alphabet**2))
scaling = (alphabet[1]-alphabet[0])/2
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) return power
power = []
for i in np.arange(0,1,0.1):
power.append(power_for_shaping(i))
plt.plot(np.arange(0,1,0.1), power)
M = 8
shaping = 2.5
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(-shapingnp.square(alphabet)))
probs = 1/denominator * np.exp(-shapingnp.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(symbols2)
print('Power: ', power)
plt.rcParams[‘figure.figsize’] = [4, 4]
plt.hist(symbols, bins=100)
plt.show()
mi_mb_64 = []
for snrdB in SNR_dBs:
sigma2 = 1/(10(snrdB/10))
sigma2 = sigma2
y = AWGN_channel(symbols, sigma2)
apps = AWGNdemapper(y, alphabet, sigma2)
xe = xesmd(apps, indices)
mi_mb_64.append((H - (xe) / np.log(2)))
256 QAM
M = 16
alphabet = np.arange(-15,16,2)
alphabet = alphabet / np.sqrt(np.mean(alphabet**2))
indices = np.random.choice(np.arange(M), n)
symbols = alphabet[indices]
print('Power: ',np.mean(np.square(symbols)))
mi_256 = []
for snrdB in SNR_dBs:
sigma2 = 1/(10**(snrdB/10))
sigma2 = sigma2
y = AWGN_channel(symbols, sigma2)
apps = AWGNdemapper(y, alphabet, sigma2)
xe = xesmd(apps, indices)
mi_256.append(2*(4 - (xe) / np.log(2)))
Plot
plt.rcParams[‘figure.figsize’] = [8, 6]
plt.plot(SNR_dBs, mi_16, label = ‘16QAM’)
plt.plot(SNR_dBs, mi_64, label = ‘64QAM’)
plt.plot(SNR_dBs, mi_mb_64, label = ‘64-MB’)
plt.plot(SNR_dBs, mi_256, label = ‘256QAM’)
plt.plot(SNR_dBs, np.log2(1+10**(SNR_dBs/10)), color=‘black’, label=‘Capacity’)
plt.legend()
plt.grid()
plt.title(‘AIR of QAM constellation’)
plt.ylabel(‘bits per channel use’)
plt.xlabel(‘SNR in dB’)
tikzplotlib.save(“QAM_MI.tex”)
SNR_dBs = np.arange(5,25)
plt.plot(SNR_dBs, np.log2(1+10**(SNR_dBs/10)), color=‘black’, label=‘
C
(
P
/
σ
2
)
C(P/σ
2
)’)
plt.plot(SNR_dBs, np.log2(1+10**(SNR_dBs/10)) - 0.5np.log2((np.pinp.e)/6) , linestyle=‘dashed’, color=‘C1’, label=‘
C
(
P
/
σ
2
)
−
1
2
log
2
π
e
6
C(P/σ
2
)−
2
1
log
2
6
πe
’)
plt.grid()
plt.ylabel(‘bits per channel use’)
plt.xlabel(‘SNR in dB’)
plt.xlim([5, 25])
plt.ylim([1, 8])
plt.title(‘AWGN channel capacity gap’)
plt.legend()
tikzplotlib.save(“SNR_GAP.tex”)
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))
解释代码，并进行注释，请注意，可以对代码进行规范化，但不要修改原本的结构。另外，将原本的信噪比SNR_dBs = np.arange(5,22)修改为仅包含5,10,15,20.增加星座点可视化，可参考下列代码：
星座图可视化
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()
本文详细介绍了NumPy库中生成随机数的方法，包括uniform函数的使用及其参数解析，以及如何通过ones和empty函数创建特定形状和类型的数组。同时，提供了多个代码示例帮助理解。
扫一扫
554
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
