[学习]M-QAM的数学原理与调制解调原理详解（仿真示例）

M-QAM的数学原理与调制解调原理详解

QAM（正交幅度调制）作为现代数字通信的核心技术，其数学原理和实现方法值得深入探讨。本文将分为数学原理、调制解调原理和实现要点三个部分进行系统阐述。

一、数学原理

1. 信号空间理论
   M-QAM信号可表示为二维信号空间中的点集：
   $$s_{mn}(t)=A_{m}g(t)\cos (2\pi f_{c}t)-B_{n}g(t)\sin (2\pi f_{c}t),m,n=1,2,...,\sqrt{M}$$
   其中：

• $A_m$, $B_n$：离散幅度值
• $g(t)$：脉冲成形函数
• $f_c$：载波频率

2. 星座点能量计算
   对于矩形星座（如16QAM），平均符号能量：
   $$E_{avg}=\frac{2d^2}{M}\sum _{i=1}^{\sqrt{M/2}}(2i-1{)}^2$$
   其中$d$为最小距离，4-QAM时$E_{avg}=2d^2$，16-QAM时$E_{avg}=10d^2$

3. 正交性证明
   载波正交性验证：
   $${\int }_0^T\cos (2\pi f_{c}t)\sin (2\pi f_{c}t)dt=0$$
   满足正交条件，实现频谱重叠而不干扰

二、调制原理

1. 调制过程数学描述
   $$s(t)=\sum _k[I_{k}g(t-k{T}_s)\cos (2\pi f_{c}t)-Q_{k}g(t-k{T}_s)\sin (2\pi f_{c}t)]$$
   其中$I_k$, $Q_k$为第k个符号的同相和正交分量

调制核心方程

$$s(t)=\text{Re}\left[\sum _k(a_k+j{b}_k)g(t-k{T}_s)e^{j2\pi f_{c}t}\right]$$

2. 实现步骤

   • 比特分组：每$k={\log }_{2}M$比特一组
   • 星座映射：查表法或计算法实现

     # 64QAM映射示例
     def map_64qam(bits):
         bits = bits.reshape(-1,6)
         I = 2*(32*(bits[:,0]^bits[:,1]^bits[:,2]) + \
             4*(bits[:,0]^bits[:,1]) + 2*bits[:,0] - 7
         Q = 2*(32*(bits[:,3]^bits[:,4]^bits[:,5]) + \
             4*(bits[:,3]^bits[:,4]) + 2*bits[:,3] - 7
         return (I + 1j*Q)/np.sqrt(42)
     

   • 脉冲成形：常用升余弦滤波器
     $$g(t)=\text{sinc}(\frac{t}{T_s})\frac{\cos (\pi \alpha t/T_s)}{1-(2\alpha t/T_s{)}^2}$$
   • 正交调制：I/Q两路分别调制

3. 调制器结构

三、解调原理

1. 解调过程数学描述
   相干解调：
   $$\{\begin{array}{l}{I}_k=\int r(t)g(t-k{T}_s)\cos (2\pi f_{c}t)dt\\ Q_k=-\int r(t)g(t-k{T}_s)\sin (2\pi f_{c}t)dt\end{array}$$

解调核心方程：
$$r_k={\int }_{k{T}_s}^{(k+1)T_s}r(t)g(t-k{T}_s)e^{-j2\pi f_{c}t}dt$$

2. 关键处理步骤

   • 载波同步：Costas环实现

     def costas_loop(signal, fs, fc):
         phase = 0
         for n in range(len(signal)):
             # 相位检测
             error = np.real(signal[n]) * np.imag(signal[n])
             # 环路滤波
             phase -= 0.01 * error
             # 相位补偿
             signal[n] *= np.exp(-1j*phase)
         return signal
     

   • 定时恢复：Gardner算法
   • 均衡处理：LMS自适应均衡
     $${\mathbf{w}}_{n+1}={\mathbf{w}}_n+\mu e_n{\mathbf{x}}_n^{\ast }$$

3. 解调器结构

四、实现要点

1. 星座图设计

   • Gray编码：相邻星座点仅1比特差异
   • 非矩形星座：如十字形星座（用于32QAM等非2^m情况）

2. 性能优化技术

   • 软判决解码：利用幅度信息提高解码性能

     def soft_demod(symbol, constellation):
         distances = np.abs(symbol - constellation)
         llr = np.log((np.min(distances[::2]) + 1e-10) - \
               np.log((np.min(distances[1::2]) + 1e-10)
         return llr
     

   • 自适应调制：根据信道质量动态调整M值

     def adapt_modulation(snr_db):
         if snr_db > 30: return 256
         elif snr_db > 20: return 64
         elif snr_db > 15: return 16
         else: return 4
     

3. 性能指标计算

   • 理论BER：
     $$P_b\approx \frac{4}{{\log }_{2}M}\left(1-\frac{1}{\sqrt{M}}\right)Q\left(\sqrt{\frac{3{\log }_{2}M}{M-1}\frac{E_b}{N_0}}\right)$$
   • EVM计算：

     def calc_evm(tx_symbols, rx_symbols):
         return 100*np.sqrt(np.mean(np.abs(tx_symbols-rx_symbols)**2)/\
                np.mean(np.abs(tx_symbols)**2))
     

五、16QAM的Python仿真实现

5.1 完整仿真代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import upfirdn

plt.close('all')
# 设置全局字体为支持中文的字体
plt.rcParams['font.sans-serif'] = ['WenQuanYi Micro Hei']  # 黑体
# 解决负号显示问题
plt.rcParams['axes.unicode_minus'] = False

# ================== 参数设置 ==================
M = 16                      # QAM阶数
bits_per_symbol = int(np.log2(M))
num_symbols = 10000         # 发送符号数
sps = 8                     # 每符号采样数
rolloff = 0.4               # 升余弦滚降因子
fs = 1e6                    # 采样率
SNR_dB = 20                 # 信噪比

# ================== 16QAM调制 ==================
def qam16_modulate(bits):
    # 将比特流重塑为每行4比特
    bits_reshaped = bits.reshape(-1, bits_per_symbol)
    
    # 生成16QAM星座图 (Gray编码)
    constellation = np.array([
        -3+3j, -3+1j, -3-3j, -3-1j,
        -1+3j, -1+1j, -1-3j, -1-1j,
         3+3j,  3+1j,  3-3j,  3-1j,
         1+3j,  1+1j,  1-3j,  1-1j
    ], dtype=np.complex128)
    
    # 归一化因子 (确保平均功率为1)
    normalization = np.sqrt(np.mean(np.abs(constellation)**2))
    constellation /= normalization
    
    # 将比特映射到符号
    symbols = np.zeros(len(bits_reshaped), dtype=np.complex128)
    for i in range(len(bits_reshaped)):
        # 将4个比特转换为0-15的整数
        symbol_index = (bits_reshaped[i, 0] << 3) | (bits_reshaped[i, 1] << 2) | (bits_reshaped[i, 2] << 1) | bits_reshaped[i, 3]
        symbols[i] = constellation[symbol_index]
    
    return symbols, constellation

# ================== 脉冲成形 ==================
def rrc_filter(sps, rolloff, span=8):
    """生成根升余弦滤波器"""
    t = np.arange(-span*sps, span*sps+1) / sps
    h = np.zeros_like(t, dtype=np.float64)
    
    # 避免除以零
    for i, t_val in enumerate(t):
        if abs(t_val) < 1e-8:  # t=0
            h[i] = 1.0 - rolloff + (4*rolloff/np.pi)
        elif abs(abs(4*rolloff*t_val) - 1.0) < 1e-8:  # 奇异点
            h[i] = (rolloff/np.sqrt(2)) * ((1+2/np.pi)*np.sin(np.pi/(4*rolloff)) + (1-2/np.pi)*np.cos(np.pi/(4*rolloff)))
        else:
            num = np.sin(np.pi*t_val*(1-rolloff)) + 4*rolloff*t_val*np.cos(np.pi*t_val*(1+rolloff))
            den = np.pi*t_val*(1 - (4*rolloff*t_val)**2)
            h[i] = num / den
    
    # 归一化滤波器能量
    h /= np.sqrt(np.sum(h**2))
    return h

# ================== 主程序 - 基带仿真 ==================
if __name__ == "__main__":
    np.random.seed(42)  # 固定随机种子
    
    # 生成随机比特流
    tx_bits = np.random.randint(0, 2, num_symbols * bits_per_symbol)
    
    # =============== 调制过程 ===============
    # 星座映射
    tx_symbols, constellation = qam16_modulate(tx_bits)
    
    # 创建根升余弦滤波器
    h_rrc = rrc_filter(sps, rolloff)
    filter_delay = (len(h_rrc) - 1) // 2  # 计算滤波器延迟
    
    # 脉冲成形
    tx_shaped = upfirdn(h_rrc, tx_symbols, up=sps)
    
    # =============== 信道传输 ===============
    # 添加高斯白噪声 (复噪声)
    signal_power = np.mean(np.abs(tx_shaped)**2)
    snr = 10**(SNR_dB/10)
    noise_power = signal_power / snr
    noise_real = np.sqrt(noise_power/2) * np.random.randn(len(tx_shaped))
    noise_imag = np.sqrt(noise_power/2) * np.random.randn(len(tx_shaped))
    rx_baseband = tx_shaped + (noise_real + 1j*noise_imag)
    
    # =============== 解调过程 ===============
    # 匹配滤波
    rx_filtered = upfirdn(h_rrc, rx_baseband)
    
    # 修正：计算总延迟（发送滤波+接收滤波）
    total_delay = 2 * filter_delay
    
    # 寻找最佳采样点 (考虑总延迟)
    start_idx = total_delay
    # 确保不会超出范围
    if start_idx + (num_symbols-1)*sps < len(rx_filtered):
        sampled_points = rx_filtered[start_idx:start_idx+num_symbols*sps:sps]
    else:
        # 如果超出范围，只取有效部分
        sampled_points = rx_filtered[start_idx::sps]
    
    # 修正：截断到与发射符号相同长度
    min_len = min(len(tx_symbols), len(sampled_points))
    sampled_points = sampled_points[:min_len]
    
    # 修正：相位补偿
    if min_len > 0:
        # 使用已知符号估计相位误差
        phase_error = np.angle(np.sum(sampled_points * np.conj(tx_symbols[:min_len])))
        compensated_points = sampled_points * np.exp(-1j * phase_error)
    else:
        compensated_points = sampled_points
    
    # 最小距离判决（使用相位补偿后的信号）
    symbol_indices = np.zeros(len(compensated_points), dtype=int)
    for i, s in enumerate(compensated_points):
        # 计算到所有星座点的距离
        distances = np.abs(s - constellation)
        # 找到最小距离的星座点索引
        symbol_indices[i] = np.argmin(distances)
    
    # 符号转比特
    rx_bits = np.zeros(len(symbol_indices) * bits_per_symbol, dtype=int)
    for i, idx in enumerate(symbol_indices):
        # 将符号索引转换为比特
        for j in range(bits_per_symbol):
            # 从最高位开始提取比特
            shift = bits_per_symbol - 1 - j
            rx_bits[i*bits_per_symbol + j] = (idx >> shift) & 1
    
    # =============== 性能分析 ===============
    # 截断到相同长度
    min_bit_len = min(len(tx_bits), len(rx_bits))
    tx_bits_trunc = tx_bits[:min_bit_len]
    rx_bits_trunc = rx_bits[:min_bit_len]
    
    # 误码率计算
    bit_errors = np.sum(tx_bits_trunc != rx_bits_trunc)
    ber = bit_errors / min_bit_len if min_bit_len > 0 else 1.0
    print(f"传输比特数: {min_bit_len}")
    print(f"误码数: {bit_errors}")
    print(f"误码率(BER): {ber:.2e}")
    
    # EVM计算（使用补偿后信号）
    if min_len > 0:
        ref_symbols = tx_symbols[:min_len]
        rx_symbols_subset = compensated_points
        
        # 计算误差矢量幅度
        evm_numerator = np.sqrt(np.mean(np.abs(rx_symbols_subset - ref_symbols)**2))
        evm_denominator = np.sqrt(np.mean(np.abs(ref_symbols)**2))
        evm = evm_numerator / evm_denominator
        print(f"EVM: {evm*100:.2f}%")
    else:
        print("EVM: N/A (无接收符号)")
    
    # 星座图绘制
    plt.figure(figsize=(14, 10))
    
    # 发射星座图
    plt.subplot(221)
    plt.scatter(np.real(tx_symbols), np.imag(tx_symbols), alpha=0.6)
    plt.title(f"发射星座图 ({len(tx_symbols)} 符号)")
    plt.grid(True)
    plt.axis('equal')
    plt.xlim(-2, 2)
    plt.ylim(-2, 2)
    
    # 接收星座图（使用补偿后信号）
    plt.subplot(222)
    if min_len > 0:
        plt.scatter(np.real(compensated_points), np.imag(compensated_points), alpha=0.6, c='r')
        plt.title(f"接收星座图 (SNR={SNR_dB}dB, BER={ber:.2e})")
        plt.grid(True)
        plt.axis('equal')
        plt.xlim(-2, 2)
        plt.ylim(-2, 2)
    else:
        plt.text(0.5, 0.5, "无接收符号", ha='center', va='center')
        plt.title("接收星座图 - 无数据")
    
    # 眼图绘制
    if len(rx_filtered) > 3 * sps:
        # 跳过初始过渡区
        stable_start = total_delay * 2
        stable_end = len(rx_filtered) - total_delay
        stable_signal = rx_filtered[stable_start:stable_end]
        
        # I路眼图
        plt.subplot(223)
        for i in range(100, min(300, len(stable_signal)//sps - 1)):
            start = i * sps
            end = start + 2 * sps
            if end < len(stable_signal):
                segment = np.real(stable_signal[start:end])
                plt.plot(np.linspace(0, 2, len(segment)), segment, 'b', alpha=0.3)
        plt.title("I路眼图")
        plt.grid(True)
        
        # Q路眼图
        plt.subplot(224)
        for i in range(100, min(300, len(stable_signal)//sps - 1)):
            start = i * sps
            end = start + 2 * sps
            if end < len(stable_signal):
                segment = np.imag(stable_signal[start:end])
                plt.plot(np.linspace(0, 2, len(segment)), segment, 'r', alpha=0.3)
        plt.title("Q路眼图")
        plt.grid(True)
    
    plt.tight_layout()
    plt.savefig('qam16_baseband_simulation_corrected.png', dpi=300)
    plt.show()

5.2 关键代码解析

1. 星座映射：

constellation = {
    0: (-3+3j), 1: (-3+1j), 2: (-3-3j), 3: (-3-1j),
    4: (-1+3j), 5: (-1+1j), 6: (-1-3j), 7: (-1-1j),
    8: (3+3j),  9: (3+1j),  10:(3-3j), 11:(3-1j),
    12:(1+3j), 13:(1+1j),  14:(1-3j), 15:(1-1j)
}

采用Gray编码排列，最小化相邻星座点的比特差异

2. 脉冲成形：

h = np.sinc(t) * np.cos(np.pi*rolloff*t) / (1 - (2*rolloff*t)**2)

升余弦滤波器设计，控制符号间干扰(ISI)

3. 最佳接收：

filtered = upfirdn(h, signal, down=sps)

匹配滤波最大化信噪比

5.3 仿真结果分析

结果输出：

传输比特数: 40000
误码数: 0
误码率(BER): 0.00e+00
EVM: 3.55%

结果可视化：

更多仿真
调制SNR_dB 参数，可以观察不同信噪比下的仿真结果。

六、性能优化方向

1. 信道均衡：

# LMS自适应均衡示例
def lms_equalizer(rx_signal, training_seq, num_taps=5, mu=0.01):
    w = np.zeros(num_taps, dtype=complex)
    for n in range(len(training_seq)-num_taps):
        x = rx_signal[n:n+num_taps]
        e = training_seq[n] - np.dot(w, x)
        w += mu * e * np.conj(x)
    return w

2. 载波同步：

   • Costas环实现
   • 频偏估计算法

3. 自适应调制：

   • 根据信道质量动态调整QAM阶数
   • 算法实现：

     def adapt_modulation(snr_db):
         if snr_db > 25: return 64
         elif snr_db > 15: return 16
         else: return 4
     

七、MATLAB验证示例

提供一个简单的matlab示例，有兴趣的可以根据这个调整一下，用于和前文Python仿真的对比。

% 16QAM调制解调示例
M = 16; 
data = randi([0 M-1],1000,1);
txSig = qammod(data,M,'UnitAveragePower',true);

% 加噪
SNR = 20;
rxSig = awgn(txSig,SNR,'measured');

% 解调
rxData = qamdemod(rxSig,M,'UnitAveragePower',true);

% 计算BER
ber = sum(data ~= rxData)/length(data);
constellation = qammod(0:M-1,M,'UnitAveragePower',true);
scatterplot(rxSig);
hold on; plot(constellation,'r*'); hold off;

八、工程实现挑战

1. 载波同步：频偏补偿精度需达$0.01\%$以内
2. 定时同步：采样时钟偏差需小于$100$ppm
3. 非线性失真：功率放大器AM-AM/AM-PM效应补偿

   def digital_predistorter(signal, pa_model):
       # 基于查找表的预失真
       return np.interp(np.abs(signal), pa_model['amp_in'], pa_model['amp_out']) * \
              np.exp(1j*np.interp(np.abs(signal), pa_model['phase_in'], pa_model['phase_out']))
   

---

代码调试是真耗费时间精力，对您有帮助请点赞关注支持~~

---

研究学习不易，点赞易。
工作生活不易，收藏易，点收藏不迷茫 ：）

---