用Python实战解析数字通信：从零绘制BPSK/QPSK/16-QAM星座图

数字通信系统的核心在于如何高效可靠地将比特流转换为电磁波信号。在这个过程中，调制技术扮演着关键角色。而星座图，则是理解各种调制方式最直观的窗口。本文将带你用Python和Matplotlib，从基本原理出发，一步步实现BPSK、QPSK和16-QAM星座图的绘制与可视化分析。

1. 理解星座图的数学基础

在开始编码前，我们需要明确几个核心概念。任何载波信号都可以表示为：

s(t) = I·cos(2πft) - Q·sin(2πft)

其中：

• I（In-phase）：同相分量
• Q（Quadrature）：正交分量
• f：载波频率

这个表达式揭示了数字调制的本质——我们实际上是通过控制I和Q这两个参数来传递信息。将不同的(I,Q)组合映射到二维平面上，就形成了所谓的星座图。

1.1 直角坐标与极坐标转换

星座图中的每个点既可以用直角坐标(I,Q)表示，也可以用极坐标(A,φ)表示，两者之间的转换关系为：

直角坐标 → 极坐标	极坐标 → 直角坐标
A = √(I² + Q²)	I = A·cos(φ)
φ = atan2(Q, I)	Q = A·sin(φ)

提示：Python中的math.atan2(y, x)函数可以正确处理所有象限的角度计算，比单纯的atan(y/x)更可靠。

1.2 常见调制方式的星座点分布

不同调制方式对应不同的星座点布局：

• BPSK：2个点，沿I轴对称分布
• QPSK：4个点，均匀分布在单位圆上
• 16-QAM：16个点，排列成4×4的方阵

这些星座点的坐标可以通过以下Python代码预先计算：

import numpy as np

# BPSK星座点
bpsk_points = np.array([[-1, 0], [1, 0]])

# QPSK星座点
qpsk_points = np.array([
    [1/np.sqrt(2), 1/np.sqrt(2)],
    [-1/np.sqrt(2), 1/np.sqrt(2)],
    [-1/np.sqrt(2), -1/np.sqrt(2)],
    [1/np.sqrt(2), -1/np.sqrt(2)]
])

# 16-QAM星座点（归一化到单位能量）
qam16_points = np.array([
    [-3, -3], [-3, -1], [-3, 1], [-3, 3],
    [-1, -3], [-1, -1], [-1, 1], [-1, 3],
    [1, -3], [1, -1], [1, 1], [1, 3],
    [3, -3], [3, -1], [3, 1], [3, 3]
]) / np.sqrt(10)

2. 搭建Python绘图环境

在开始绘制星座图前，我们需要配置合适的Python环境。推荐使用Anaconda创建虚拟环境：

conda create -n digital_com python=3.8
conda activate digital_com
pip install numpy matplotlib

2.1 基础绘图函数设计

我们将创建一个通用的星座图绘制函数，支持多种调制方式：

import matplotlib.pyplot as plt

def plot_constellation(points, title="", noise_std=0):
    """绘制星座图
    
    参数：
        points: 星座点坐标，形状为(N,2)的数组
        title: 图表标题
        noise_std: 高斯噪声标准差
    """
    plt.figure(figsize=(8, 8))
    
    # 添加噪声（如果指定）
    if noise_std > 0:
        noisy_points = points + np.random.randn(*points.shape) * noise_std
        plt.scatter(noisy_points[:,0], noisy_points[:,1], 
                   c='red', alpha=0.3, label='含噪声信号')
    
    # 绘制理想星座点
    plt.scatter(points[:,0], points[:,1], 
               c='blue', marker='x', s=100, linewidth=2, label='理想位置')
    
    # 设置坐标轴
    plt.axhline(0, color='gray', linestyle='--', linewidth=0.5)
    plt.axvline(0, color='gray', linestyle='--', linewidth=0.5)
    plt.grid(True, linestyle=':', alpha=0.5)
    plt.axis('equal')
    
    # 添加标签和标题
    plt.xlabel('同相分量 (I)')
    plt.ylabel('正交分量 (Q)')
    plt.title(title)
    plt.legend()
    
    return plt

3. 绘制基础星座图

3.1 BPSK星座图实现

BPSK是最简单的数字调制方式，每个符号只携带1比特信息。让我们绘制其星座图：

# 绘制BPSK星座图
plot_constellation(bpsk_points, "BPSK星座图")
plt.show()

执行结果将显示I轴上的两个对称点。我们可以添加噪声观察效果：

# 带噪声的BPSK星座图
plot_constellation(bpsk_points, "BPSK星座图（含噪声）", noise_std=0.2)
plt.show()

3.2 QPSK星座图实现

QPSK通过四个相位状态传输2比特信息，频谱效率是BPSK的两倍：

# 绘制QPSK星座图
plot_constellation(qpsk_points, "QPSK星座图")
plt.show()

# 带噪声的QPSK星座图
plot_constellation(qpsk_points, "QPSK星座图（含噪声）", noise_std=0.15)
plt.show()

观察噪声影响时，注意四个象限的点如何形成"云团"。

3.3 16-QAM星座图实现

16-QAM通过幅度和相位的组合实现更高频谱效率：

# 绘制16-QAM星座图
plot_constellation(qam16_points, "16-QAM星座图")
plt.show()

# 带噪声的16-QAM星座图
plot_constellation(qam16_points, "16-QAM星座图（含噪声）", noise_std=0.1)
plt.show()

注意观察噪声如何使密集的星座点相互重叠，这是高阶调制对信道质量要求更高的直观体现。

4. 星座图的高级可视化

4.1 动态噪声演示

为了更直观地理解噪声影响，我们可以创建动态演示：

from matplotlib.animation import FuncAnimation

def animate_noise(points, title=""):
    fig, ax = plt.subplots(figsize=(8, 8))
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.grid(True)
    ax.set_title(title)
    ax.set_xlabel('I')
    ax.set_ylabel('Q')
    
    # 绘制理想点
    ideal = ax.scatter(points[:,0], points[:,1], 
                      c='blue', marker='x', s=100, label='理想位置')
    
    # 初始化噪声点
    noisy = ax.scatter([], [], c='red', alpha=0.5, label='含噪声信号')
    
    def update(frame):
        # 更新噪声点
        noise_std = 0.05 + frame * 0.005
        noisy_points = points + np.random.randn(*points.shape) * noise_std
        noisy.set_offsets(noisy_points)
        ax.set_title(f"{title} (噪声标准差: {noise_std:.3f})")
        return noisy,
    
    ani = FuncAnimation(fig, update, frames=100, interval=100, blit=True)
    plt.legend()
    plt.close()
    return ani

# 生成QPSK动态噪声演示
ani = animate_noise(qpsk_points, "QPSK星座图噪声演示")
from IPython.display import HTML
HTML(ani.to_jshtml())

4.2 星座图与误码率关系分析

星座点之间的最小欧氏距离直接影响系统抗噪声能力。我们可以计算并可视化这一关键参数：

def plot_min_distance(points, title=""):
    # 计算所有点对之间的距离
    dists = []
    n = len(points)
    for i in range(n):
        for j in range(i+1, n):
            dist = np.linalg.norm(points[i] - points[j])
            dists.append(dist)
    
    min_dist = min(dists)
    
    # 绘制距离分布直方图
    plt.figure(figsize=(10, 5))
    plt.hist(dists, bins=20, alpha=0.7)
    plt.axvline(min_dist, color='red', linestyle='--', 
               label=f'最小距离 = {min_dist:.3f}')
    plt.xlabel('星座点间欧氏距离')
    plt.ylabel('频次')
    plt.title(f'{title} - 星座点距离分布')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

# 分析不同调制方式的距离分布
plot_min_distance(bpsk_points, "BPSK")
plot_min_distance(qpsk_points, "QPSK")
plot_min_distance(qam16_points, "16-QAM")

从结果中可以明显看出，随着调制阶数提高，最小距离减小，系统对噪声更加敏感。

5. 实际应用与扩展

5.1 星座图在信号质量评估中的应用

星座图不仅能显示理想信号状态，还能反映各种信号损伤：

• 相位噪声：星座点呈现环形扩散
• 幅度压缩：外圈星座点向内收缩
• IQ不平衡：星座图整体变形

我们可以模拟这些效应：

def plot_impairments(points, title=""):
    plt.figure(figsize=(15, 5))
    
    # 相位噪声
    plt.subplot(131)
    phase_noise = np.random.randn(len(points)) * 0.3
    rotated = points @ np.array([[np.cos(phase_noise), -np.sin(phase_noise)],
                               [np.sin(phase_noise), np.cos(phase_noise)]]).T
    plt.scatter(rotated[:,0], rotated[:,1], c='red', alpha=0.5)
    plt.scatter(points[:,0], points[:,1], c='blue', marker='x', s=50)
    plt.title("相位噪声效应")
    plt.grid(True)
    
    # 幅度压缩
    plt.subplot(132)
    compressed = points * (0.8 + 0.1*np.abs(points))
    plt.scatter(compressed[:,0], compressed[:,1], c='red', alpha=0.5)
    plt.scatter(points[:,0], points[:,1], c='blue', marker='x', s=50)
    plt.title("幅度压缩效应")
    plt.grid(True)
    
    # IQ不平衡
    plt.subplot(133)
    imbalanced = points @ np.array([[1, 0.2], [0.1, 0.9]])
    plt.scatter(imbalanced[:,0], imbalanced[:,1], c='red', alpha=0.5)
    plt.scatter(points[:,0], points[:,1], c='blue', marker='x', s=50)
    plt.title("IQ不平衡效应")
    plt.grid(True)
    
    plt.suptitle(f"{title} - 常见信号损伤模拟")
    plt.tight_layout()
    plt.show()

plot_impairments(qpsk_points, "QPSK")
plot_impairments(qam16_points, "16-QAM")

5.2 扩展到其他调制方式

基于相同的框架，我们可以轻松实现更高阶的调制方式：

# 64-QAM星座点生成
qam64_points = np.array([
    [-7, -7], [-7, -5], [-7, -3], [-7, -1], [-7, 1], [-7, 3], [-7, 5], [-7, 7],
    [-5, -7], [-5, -5], [-5, -3], [-5, -1], [-5, 1], [-5, 3], [-5, 5], [-5, 7],
    [-3, -7], [-3, -5], [-3, -3], [-3, -1], [-3, 1], [-3, 3], [-3, 5], [-3, 7],
    [-1, -7], [-1, -5], [-1, -3], [-1, -1], [-1, 1], [-1, 3], [-1, 5], [-1, 7],
    [1, -7], [1, -5], [1, -3], [1, -1], [1, 1], [1, 3], [1, 5], [1, 7],
    [3, -7], [3, -5], [3, -3], [3, -1], [3, 1], [3, 3], [3, 5], [3, 7],
    [5, -7], [5, -5], [5, -3], [5, -1], [5, 1], [5, 3], [5, 5], [5, 7],
    [7, -7], [7, -5], [7, -3], [7, -1], [7, 1], [7, 3], [7, 5], [7, 7]
]) / np.sqrt(42)

plot_constellation(qam64_points, "64-QAM星座图")
plot_min_distance(qam64_points, "64-QAM")

在实际项目中，这些可视化工具可以帮助工程师快速诊断系统问题。例如，当观察到64-QAM星座图外圈点明显内缩时，可能提示功率放大器进入了非线性区。