[学习]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] - 7Q = 2*(32*(bits[:,3]^bits[:,4]^bits[:,5]) + \4*(bits[:,3]^bits[:,4]) + 2*bits[:,3] - 7return (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 = 0for 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 256elif snr_db > 20: return 64elif snr_db > 15: return 16else: 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 upfirdnplt.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) / spsh = np.zeros_like(t, dtype=np.float64)# 避免除以零for i, t_val in enumerate(t):if abs(t_val) < 1e-8:  # t=0h[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 / snrnoise_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 - jrx_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.0print(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_denominatorprint(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 * 2stable_end = len(rx_filtered) - total_delaystable_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 * spsend = start + 2 * spsif 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 * spsend = start + 2 * spsif 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 64elif snr_db > 15: return 16else: 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']))
   

---

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

---

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

---