动手仿真：用Python从零实现BPSK/QPSK/MSK调制与解调（附完整代码）-程序员充电站

用Python从零实现BPSK/QPSK/MSK调制与解调：通信工程师的代码实践指南

从理论到实践：数字调制技术的Python实现

作为一名通信工程师，我经常需要在理论知识和实际实现之间架起桥梁。数字调制技术是无线通信系统的核心，但教科书上的公式往往让人难以直观理解其工作原理。本文将带您用Python从零实现BPSK、QPSK和MSK调制解调系统，通过代码让抽象概念变得触手可及。

在开始编码前，我们需要明确几个关键概念。数字调制本质上是用数字信号控制载波的某些特性（幅度、频率或相位）的过程。BPSK（二进制相移键控）通过改变载波相位来传输信息，QPSK（四相相移键控）是其扩展，每个符号携带2比特信息。MSK（最小频移键控）则是一种特殊的连续相位频移键控，具有频谱效率高和包络恒定的优点。

1. 环境准备与基础构建

1.1 Python环境配置

首先确保您已安装以下Python库：

import numpy as np import matplotlib.pyplot as plt from scipy import signal from scipy.stats import norm

这些库将帮助我们完成信号生成、处理和可视化工作。我推荐使用Jupyter Notebook进行实验，因为它支持交互式开发和即时可视化。

1.2 基础信号生成

我们先创建一个简单的函数来生成随机比特流：

def generate_bits(n_bits): """生成随机比特序列""" return np.random.randint(0, 2, n_bits)

接下来，定义一些基本参数：

# 系统参数 bit_rate = 1000 # 比特率 (bps) sample_rate = 10000 # 采样率 (Hz) fc = 2000 # 载波频率 (Hz) n_bits = 100 # 传输的比特数 Tb = 1/bit_rate # 比特周期 (s)

2. BPSK调制与解调实现

2.1 BPSK调制器设计

BPSK是最简单的数字调制方式，用0°相位表示"0"，180°相位表示"1"。实现代码如下：

def bpsk_modulate(bits, fc, sample_rate, bit_rate): """BPSK调制""" samples_per_bit = int(sample_rate / bit_rate) t = np.arange(0, len(bits)*Tb, 1/sample_rate) carrier = np.cos(2*np.pi*fc*t) # 将比特映射到相位：0->1, 1->-1 symbols = 1 - 2*bits modulated = np.repeat(symbols, samples_per_bit) * carrier return modulated, t, carrier

2.2 BPSK解调技术

BPSK解调通常采用相干解调（同步检测）方式：

def bpsk_demodulate(signal, fc, sample_rate, bit_rate): """BPSK解调""" samples_per_bit = int(sample_rate / bit_rate) t = np.arange(0, len(signal)/sample_rate, 1/sample_rate) carrier = np.cos(2*np.pi*fc*t[:len(signal)]) # 相干解调 demodulated = signal * carrier # 积分判决 bits_recovered = np.zeros(len(signal)//samples_per_bit) for i in range(len(bits_recovered)): start = i * samples_per_bit end = (i+1) * samples_per_bit bits_recovered[i] = np.sum(demodulated[start:end]) < 0 return bits_recovered

2.3 BPSK性能测试与可视化

让我们测试一下我们的BPSK实现：

bits = generate_bits(n_bits) bpsk_signal, t, carrier = bpsk_modulate(bits, fc, sample_rate, bit_rate) bits_recovered = bpsk_demodulate(bpsk_signal, fc, sample_rate, bit_rate) # 绘制前10个比特的波形 plt.figure(figsize=(12, 6)) plt.subplot(3,1,1) plt.plot(t[:10*sample_rate//bit_rate], bits[:10].repeat(samples_per_bit)) plt.title('原始比特流') plt.subplot(3,1,2) plt.plot(t[:10*sample_rate//bit_rate], bpsk_signal[:10*sample_rate//bit_rate]) plt.title('BPSK调制信号') plt.subplot(3,1,3) plt.plot(t[:10*sample_rate//bit_rate], (carrier * (1-2*bits[:10].repeat(samples_per_bit)))[:10*sample_rate//bit_rate]) plt.title('载波相位反转') plt.tight_layout() plt.show()

注意：在实际系统中，载波同步是解调的关键。我们这里假设接收端已经完美同步，实际应用中需要使用锁相环等技术实现载波恢复。

3. QPSK调制系统实现

3.1 QPSK调制原理

QPSK通过将两个比特组合为一个符号，使频谱效率提高一倍。四个相位分别对应00、01、10、11：

def qpsk_modulate(bits, fc, sample_rate, bit_rate): """QPSK调制""" # 确保比特数为偶数 if len(bits) % 2 != 0: bits = np.append(bits, 0) samples_per_symbol = int(2 * sample_rate / bit_rate) t = np.arange(0, len(bits)//2 * 2/bit_rate, 1/sample_rate) # 比特到符号映射 (Gray编码) symbols = np.zeros(len(bits)//2, dtype=complex) for i in range(0, len(bits), 2): if bits[i] == 0 and bits[i+1] == 0: # 00 -> 1+j symbols[i//2] = 1 + 1j elif bits[i] == 0 and bits[i+1] == 1: # 01 -> -1+j symbols[i//2] = -1 + 1j elif bits[i] == 1 and bits[i+1] == 0: # 10 -> 1-j symbols[i//2] = 1 - 1j else: # 11 -> -1-j symbols[i//2] = -1 - 1j # 生成同相和正交分量 I = np.real(symbols).repeat(samples_per_symbol) Q = np.imag(symbols).repeat(samples_per_symbol) # 调制 carrier_I = np.cos(2*np.pi*fc*t[:len(I)]) carrier_Q = np.sin(2*np.pi*fc*t[:len(Q)]) modulated = I * carrier_I - Q * carrier_Q return modulated, t[:len(modulated)], symbols

3.2 QPSK解调实现

QPSK解调需要同时处理同相和正交分量：

def qpsk_demodulate(signal, fc, sample_rate, bit_rate): """QPSK解调""" samples_per_symbol = int(2 * sample_rate / bit_rate) t = np.arange(0, len(signal)/sample_rate, 1/sample_rate) # 相干解调 I = signal * np.cos(2*np.pi*fc*t[:len(signal)]) Q = -signal * np.sin(2*np.pi*fc*t[:len(signal)]) # 低通滤波 b, a = signal.butter(5, 0.2, 'low') I_filtered = signal.lfilter(b, a, I) Q_filtered = signal.lfilter(b, a, Q) # 符号判决 bits_recovered = np.zeros(len(signal)//samples_per_symbol * 2) for i in range(len(signal)//samples_per_symbol): start = i * samples_per_symbol end = (i+1) * samples_per_symbol avg_I = np.mean(I_filtered[start:end]) avg_Q = np.mean(Q_filtered[start:end]) # Gray解码 if avg_I > 0 and avg_Q > 0: # 1+j -> 00 bits_recovered[2*i] = 0 bits_recovered[2*i+1] = 0 elif avg_I < 0 and avg_Q > 0: # -1+j -> 01 bits_recovered[2*i] = 0 bits_recovered[2*i+1] = 1 elif avg_I > 0 and avg_Q < 0: # 1-j -> 10 bits_recovered[2*i] = 1 bits_recovered[2*i+1] = 0 else: # -1-j -> 11 bits_recovered[2*i] = 1 bits_recovered[2*i+1] = 1 return bits_recovered

3.3 QPSK星座图绘制

星座图是评估调制质量的重要工具：

def plot_constellation(symbols, title): """绘制星座图""" plt.figure(figsize=(6,6)) plt.scatter(np.real(symbols), np.imag(symbols)) plt.axhline(0, color='gray', linestyle='--') plt.axvline(0, color='gray', linestyle='--') plt.xlabel('In-phase') plt.ylabel('Quadrature') plt.title(title) plt.grid(True) plt.axis('equal') plt.show() # 测试QPSK并绘制星座图 bits = generate_bits(n_bits) qpsk_signal, t, symbols = qpsk_modulate(bits, fc, sample_rate, bit_rate) plot_constellation(symbols, 'QPSK星座图')

4. MSK调制技术实现

4.1 MSK调制原理

MSK是一种特殊的连续相位频移键控(CPFSK)，具有恒定包络特性：

def msk_modulate(bits, fc, sample_rate, bit_rate): """MSK调制""" samples_per_bit = int(sample_rate / bit_rate) t = np.arange(0, len(bits)*Tb, 1/sample_rate) # 差分编码 diff_bits = np.zeros_like(bits) diff_bits[0] = bits[0] for i in range(1, len(bits)): diff_bits[i] = bits[i] ^ diff_bits[i-1] # 生成相位轨迹 phase = np.zeros(len(t)) for i in range(len(bits)): start = i * samples_per_bit end = (i+1) * samples_per_bit if diff_bits[i] == 1: phase[start:end] = np.pi/2 * (t[start:end] - i*Tb)/Tb else: phase[start:end] = -np.pi/2 * (t[start:end] - i*Tb)/Tb # 累积相位 accumulated_phase = np.cumsum(phase) * (1/sample_rate) # 生成调制信号 modulated = np.cos(2*np.pi*fc*t + accumulated_phase) return modulated, t, phase

4.2 MSK解调实现

MSK解调可以采用相干或非相干方式，这里实现相干解调：

def msk_demodulate(signal, fc, sample_rate, bit_rate): """MSK相干解调""" samples_per_bit = int(sample_rate / bit_rate) t = np.arange(0, len(signal)/sample_rate, 1/sample_rate) # 同相和正交分量 I = signal * np.cos(2*np.pi*fc*t[:len(signal)]) Q = signal * np.sin(2*np.pi*fc*t[:len(signal)]) # 低通滤波 b, a = signal.butter(5, 0.2, 'low') I_filtered = signal.lfilter(b, a, I) Q_filtered = signal.lfilter(b, a, Q) # 比特判决 bits_recovered = np.zeros(len(signal)//samples_per_bit) for i in range(len(bits_recovered)): start = i * samples_per_bit end = (i+1) * samples_per_bit # 比较I和Q通道的能量 energy_I = np.sum(I_filtered[start:end]**2) energy_Q = np.sum(Q_filtered[start:end]**2) bits_recovered[i] = energy_Q > energy_I # 差分解码 decoded_bits = np.zeros_like(bits_recovered) decoded_bits[0] = bits_recovered[0] for i in range(1, len(bits_recovered)): decoded_bits[i] = bits_recovered[i] ^ bits_recovered[i-1] return decoded_bits

4.3 MSK信号特性分析

MSK的相位连续性和频谱特性是其最大优势：

bits = generate_bits(20) # 少量比特便于观察 msk_signal, t, phase = msk_modulate(bits, fc, sample_rate, bit_rate) plt.figure(figsize=(12, 8)) plt.subplot(3,1,1) plt.plot(t[:len(msk_signal)], msk_signal) plt.title('MSK调制信号') plt.subplot(3,1,2) plt.plot(t[:len(msk_signal)], phase) plt.title('MSK相位轨迹') plt.subplot(3,1,3) plt.specgram(msk_signal, Fs=sample_rate, NFFT=1024) plt.title('MSK频谱') plt.tight_layout() plt.show()

5. 性能比较与误码率分析

5.1 误码率测试框架

为了比较不同调制技术的性能，我们实现一个误码率测试框架：

def calculate_ber(mod_func, demod_func, snr_db, n_bits=1000, fc=2000, sample_rate=10000, bit_rate=1000): """计算给定SNR下的误码率""" bits = generate_bits(n_bits) modulated = mod_func(bits, fc, sample_rate, bit_rate)[0] # 添加高斯白噪声 signal_power = np.mean(modulated**2) noise_power = signal_power / (10**(snr_db/10)) noise = np.random.normal(0, np.sqrt(noise_power), len(modulated)) noisy_signal = modulated + noise # 解调 bits_recovered = demod_func(noisy_signal, fc, sample_rate, bit_rate) # 确保长度一致 min_len = min(len(bits), len(bits_recovered)) errors = np.sum(bits[:min_len] != bits_recovered[:min_len]) return errors / min_len

5.2 三种调制技术性能比较

snr_range = np.arange(0, 11, 1) ber_bpsk = [] ber_qpsk = [] ber_msk = [] for snr in snr_range: ber_bpsk.append(calculate_ber(bpsk_modulate, bpsk_demodulate, snr)) ber_qpsk.append(calculate_ber(qpsk_modulate, qpsk_demodulate, snr)) ber_msk.append(calculate_ber(msk_modulate, msk_demodulate, snr)) # 绘制BER曲线 plt.figure(figsize=(10,6)) plt.semilogy(snr_range, ber_bpsk, 'o-', label='BPSK') plt.semilogy(snr_range, ber_qpsk, 's-', label='QPSK') plt.semilogy(snr_range, ber_msk, 'd-', label='MSK') plt.xlabel('SNR (dB)') plt.ylabel('Bit Error Rate') plt.title('不同调制技术的误码率比较') plt.grid(True) plt.legend() plt.show()

5.3 频谱效率分析

def plot_psd(signal, sample_rate, label): """绘制功率谱密度""" f, Pxx = signal.welch(signal, fs=sample_rate, nperseg=1024) plt.semilogy(f, Pxx, label=label) plt.figure(figsize=(10,6)) bits = generate_bits(1000) plot_psd(bpsk_modulate(bits, fc, sample_rate, bit_rate)[0], sample_rate, 'BPSK') plot_psd(qpsk_modulate(bits, fc, sample_rate, bit_rate)[0], sample_rate, 'QPSK') plot_psd(msk_modulate(bits, fc, sample_rate, bit_rate)[0], sample_rate, 'MSK') plt.xlim(0, 5000) plt.xlabel('Frequency (Hz)') plt.ylabel('Power Spectral Density') plt.title('不同调制技术的频谱特性') plt.grid(True) plt.legend() plt.show()

6. 实际应用中的考量与优化

6.1 载波同步问题

在实际系统中，接收端需要从接收信号中恢复载波。一个简单的方法是使用平方环：

def carrier_recovery(signal, fc, sample_rate): """平方环载波恢复""" # 平方处理 squared = signal**2 # 带通滤波提取2fc分量 b, a = signal.butter(5, [1.9*fc, 2.1*fc], 'bandpass', fs=sample_rate) filtered = signal.lfilter(b, a, squared) # 分频 recovered_carrier = np.cos(np.angle(signal.hilbert(filtered))/2) return recovered_carrier

6.2 定时同步实现

比特定时同步可以使用早迟门算法：

def timing_recovery(signal, samples_per_bit): """早迟门定时同步""" # 计算能量 energy = signal**2 # 早采样点 early_samples = energy[:-2*samples_per_bit:samples_per_bit//2] # 迟采样点 late_samples = energy[samples_per_bit//2:-samples_per_bit:samples_per_bit//2] # 计算误差 error = np.sum(early_samples - late_samples) return error

6.3 信道均衡技术

多径信道会导致码间干扰，可以使用简单的线性均衡器：

def linear_equalizer(signal, taps=5): """线性均衡器""" # LMS自适应算法 mu = 0.01 # 步长 w = np.zeros(taps) for i in range(taps, len(signal)): x = signal[i-taps:i] # 这里假设已知训练序列，实际中需要前导码 d = signal[i] # 理想情况下 y = np.dot(w, x) e = d - y w = w + mu * e * x return w

7. 完整系统集成与测试

7.1 端到端通信系统模拟

将各个模块组合成完整系统：

def simulate_communication_system(mod_type='bpsk', snr_db=10, n_bits=1000): """模拟完整通信系统""" # 生成随机数据 bits = generate_bits(n_bits) # 调制 if mod_type == 'bpsk': modulated, t, _ = bpsk_modulate(bits, fc, sample_rate, bit_rate) elif mod_type == 'qpsk': modulated, t, _ = qpsk_modulate(bits, fc, sample_rate, bit_rate) else: # msk modulated, t, _ = msk_modulate(bits, fc, sample_rate, bit_rate) # 添加噪声 signal_power = np.mean(modulated**2) noise_power = signal_power / (10**(snr_db/10)) noise = np.random.normal(0, np.sqrt(noise_power), len(modulated)) received = modulated + noise # 解调 if mod_type == 'bpsk': recovered = bpsk_demodulate(received, fc, sample_rate, bit_rate) elif mod_type == 'qpsk': recovered = qpsk_demodulate(received, fc, sample_rate, bit_rate) else: # msk recovered = msk_demodulate(received, fc, sample_rate, bit_rate) # 计算误码率 min_len = min(len(bits), len(recovered)) errors = np.sum(bits[:min_len] != recovered[:min_len]) ber = errors / min_len return ber, bits, recovered # 测试BPSK系统 ber, tx_bits, rx_bits = simulate_communication_system('bpsk', 8) print(f"BPSK系统误码率: {ber:.2e}")

7.2 眼图分析

眼图是评估数字通信系统性能的重要工具：

def plot_eye_diagram(signal, samples_per_bit, title, offset=0, n_eyes=5): """绘制眼图""" span = 2 * samples_per_bit eye = np.zeros((n_eyes, span)) for i in range(n_eyes): start = offset + i * samples_per_bit end = start + span eye[i,:] = signal[start:end] plt.figure(figsize=(10,6)) for i in range(n_eyes): plt.plot(np.arange(span)/samples_per_bit, eye[i,:]) plt.title(title) plt.xlabel('Time (bit periods)') plt.ylabel('Amplitude') plt.grid(True) plt.show() # 生成BPSK眼图 bits = generate_bits(100) bpsk_signal, t, _ = bpsk_modulate(bits, fc, sample_rate, bit_rate) plot_eye_diagram(bpsk_signal, int(sample_rate/bit_rate), 'BPSK眼图')

7.3 实际应用建议

根据我们的实现和分析，在实际系统设计中应考虑以下因素：

频谱效率：QPSK比BPSK高一倍，MSK具有更好的频谱滚降特性
功率效率：BPSK在低SNR下表现最好，MSK次之
实现复杂度：BPSK最简单，MSK需要处理连续相位
同步要求：QPSK和MSK对载波和定时同步更敏感

在带宽受限的应用中（如卫星通信），QPSK通常是更好的选择。而在功率受限且需要恒定包络的场景（如移动通信），MSK可能更合适。

动手仿真：用Python从零实现BPSK/QPSK/MSK调制与解调（附完整代码）