简介：完整OFDM仿真程序面向通信工程相关专业学生、课程设计者及科研人员，用于系统理解OFDM系统建模流程、调制解调原理及关键参数设置。压缩包共5个m文件，大小约7KB，涵盖主程序、调制模块与解调模块，支持BPSK、QPSK、16QAM、64QAM四种映射方式，并集成循环前缀添加、加窗处理、频谱图绘制、PAPR计算及误码率统计等核心环节，覆盖从发射端到接收端的完整链路。已有2537人学习使用，资料可直接在MATLAB中运行，便于对比不同调制方式下的频谱特性与误码性能，也可观察加窗前后的频谱差异。代码结构清晰，模块划分明确，每个功能独立成m文件，方便读者按需调用与修改；适合课程设计、毕业设计或OFDM原理验证，也可作为二次开发的基础框架，帮助读者快速掌握OFDM仿真实现要点，提升无线通信系统建模能力。
1. 完整的OFDM仿真程序：从一条看似简单却极易翻车的链路说起
一套完整的OFDM仿真程序，包含QPSK和16QAM调制，听起来是通信仿真里最基础的练习。但真正动手写过的人都知道，这个“基础练习”从比特流走到星座图再回到误码率曲线，中间每一步都有让你怀疑人生的细节：FFT的缩放因子、循环前缀的长度、导频插值的边界、EbN0和SNR的换算，任何一处没对齐，整条链路的表现就会莫名其妙变差。这套仿真程序解决的核心问题，就是帮你在一套可运行的框架里，同时看清OFDM的收发结构和两种调制方式的行为差异。它适合正在做课程设计、项目预研，或者想验证一个想法又不想从零搭框架的开发者。我下面讲的是我用过的、能跑通也经得起分析的一套做法，顺着链路一步步拆给你看。
2. 先把OFDM仿真链路拆开：从比特流到星座点，每一步在做什么
2.1 发射端链路：编码、映射、串并变换、IFFT、加CP
OFDM发射端做的事情，本质上就是把一长串比特流按块切成多个并行的窄带子载波来发送。先经过调制映射（QPSK或16QAM），把比特变成复数符号；然后做串并变换，把这些符号分配到N个并行的子载波上；对这个频域向量做IFFT，变成时域波形；最后在时域符号前面加上循环前缀，再串行发送出去。
很多刚开始做仿真的人，会把注意力放在调制映射上，觉得星座图画出来就完成了。但实际上OFDM仿真的核心在串并变换和IFFT这一步：它决定了每个子载波在频谱上的位置，也决定了接收端该怎么把信号解出来。理解了发射端的五个环节，后续接收方向的所有操作都是它的逆过程。
2.2 最小可运行的OFDM发射机：参数表与第一版代码
动手之前先定一组参数。这组参数不是随便拍的，它决定了仿真能覆盖什么样的信道场景、频谱利用率是多少、接收端同步的容错空间有多大。我常用的参数如下表，算是一个比较标准的基带OFDM配置。
| 参数 | 取值 | 说明 |
|---|---|---|
| N_FFT | 64 | IFFT/FFT变换点数，即子载波总数 |
| N_CP | 16 | 循环前缀长度，覆盖的信道时延扩展上限 |
| N_data | 48 | 数据子载波数 |
| N_pilot | 16 | 导频子载波数，每4个子载波插入一个 |
| 子载波间隔 | 15 kHz | 常见于LTE的配置 |
| 采样率 | 960 kHz | N_FFT × 子载波间隔 |
| 调制方式 | QPSK / 16QAM | 均使用Gray映射 |
发射机的第一个可用版本，我建议用QPSK先跑通链路。注意下面代码里我用了norm='ortho'参数，这保证了时域信号功率和频域符号功率一致，后面的EbN0换算才不用反复补偿FFT的缩放因子。
import numpy as np N_FFT = 64 N_CP = 16 PILOT_STEP = 4 pilot_idx = np.arange(0, N_FFT, PILOT_STEP) data_idx = np.array([i for i in range(N_FFT) if i not in pilot_idx]) N_DATA = len(data_idx) PILOT_VALUE = (1 + 1j) / np.sqrt(2) def generate_qpsk_symbols(bits): # 每2比特映射一个QPSK符号，比特0对应+1，比特1对应-1 bits = bits[:2 * (len(bits) // 2)] i_bits = bits[0::2] q_bits = bits[1::2] symbols = ((2 * i_bits - 1) + 1j * (2 * q_bits - 1)) / np.sqrt(2) return symbols def ofdm_tx(bits): symbols = generate_qpsk_symbols(bits) n_symbols = len(symbols) // N_DATA symbol_frame = symbols[:n_symbols * N_DATA].reshape(n_symbols, N_DATA) X = np.zeros((n_symbols, N_FFT), dtype=complex) X[:, data_idx] = symbol_frame X[:, pilot_idx] = PILOT_VALUE # ortho归一化让IFFT前后能量一致，避免FFT缩放因子干扰SNR换算 x_time = np.fft.ifft(X, axis=1, norm='ortho') # 加循环前缀：把每个符号的最后N_CP个采样复制到符号头部 x_cp = np.concatenate([x_time[:, -N_CP:], x_time], axis=1) return x_cp这段代码的逻辑很直白：先把比特切成QPSK符号，然后按N_DATA个一组排成帧，填入数据子载波位置，导频位置填入固定值，其余子载波补零。IFFT之后，每个OFDM符号的时域长度是N_FFT个采样点，加上CP后变成N_FFT+N_CP个采样点。
参数说明里最需要注意的是norm='ortho'。NumPy默认的ifft带1/N缩放，fft不带缩放，两者混用会导致频域幅度差N倍。很多仿真程序跑出来星座图幅度不对，都是栽在这一步。用norm='ortho'之后，正反变换互为共轭转置关系，能量守恒，信噪比定义清爽得多。另外PILOT_VALUE的功率要设为1，和数据符号的平均功率对齐，否则信道估计出来的幅度会整体偏大或偏小。
2.3 接收端链路：同步、去CP、FFT、均衡、解映射
发射端的每个环节，接收端几乎都有一个逆操作。先做符号同步找到OFDM符号的起始位置，去掉循环前缀，对每个符号做FFT变换回频域，然后做信道估计与均衡，最后把数据子载波上的符号取出来解映射成比特。
在这个基础版本里，我们先假设同步完全正确、信道是理想AWGN（即信道响应H=1），那么接收端只需要去CP、FFT、取数据子载波、解映射四个动作。对应的代码如下：
def ofdm_rx_awgn(x_cp): n_symbols = x_cp.shape[0] x_time = x_cp[:, N_CP:] # 去循环前缀 Y = np.fft.fft(x_time, axis=1, norm='ortho') # 理想信道：H=1，不需要估计和均衡 Y_data = Y[:, data_idx] return Y_data去CP之后数组形状从(n_symbols, N_FFT+N_CP)变回(n_symbols, N_FFT)，FFT之后Y的每一个列索引对应一个子载波位置，X[0]对应直流子载波。然后用data_idx把数据子载波捞出来。这个版本的接收端没有做导频提取，因为AWGN信道下不需要信道估计。
这里有一个隐蔽的坑：如果发射端用默认ifft（带1/N）而接收端用默认fft（不带缩放），那么Y会等于X乘以N倍，星座图整体放大，看起来还能解调，但噪声也被放大，误码率会完全对不上理论值。所以再次强调，正反变换要么都用ortho，要么手动补偿1/N或N因子。
2.4 链路参数的联动关系：子载波数、CP长度、导频间隔怎么定
OFDM的参数不是独立决定的，它们之间互相约束。子载波间隔决定了一个OFDM符号的持续时间，N_FFT越大，符号越长，频谱效率越高，但对相位噪声和频偏越敏感。循环前缀长度必须大于信道的最大时延扩展，否则产生符号间干扰；但CP越长，信噪比开销越大，因为CP本身不携带新信息。导频间隔则由信道的相干带宽决定，信道频率选择性越强，导频就要越密。
我的习惯是先用AWGN信道把整条链路跑通，验证调制解调和OFDM收发没问题，再引入多径信道和信道估计。不要一上来就全信道全流程，否则出了问题你分不清是调制写错，还是IFFT缩放错了，还是信道估计错了。分阶段验证，是通信仿真里最务实的排错思路。
3. QPSK与16QAM调制映射：星座图、功率归一化与解调实现
3.1 QPSK映射与解调：映射表与最小命令
QPSK每2比特映射一个符号，四个星座点分别对应四个相位。最常见的Gray映射方式是：I路比特决定实部，Q路比特决定虚部，比特0映射为+1，比特1映射为-1，最后除以√2让符号平均功率为1。这样星座点落在(±1/√2, ±1/√2)上，相邻星座点之间只差1个比特，误码性能最优。
def qpsk_demod(symbols): # 硬判决：实部大于0判为0，否则判为1；虚部同理 i_bits = (symbols.real > 0).astype(int) q_bits = (symbols.imag > 0).astype(int) bits = np.zeros(2 * len(symbols), dtype=int) bits[0::2] = i_bits bits[1::2] = q_bits return bits这个解调器在AWGN信道下就是最优的，因为QPSK的四个星座点等概率，判决区域就是I-Q平面的四个象限。注意解映射得到的比特顺序必须和发射端一致。
3.2 16QAM映射与功率归一化：为什么必须除以sqrt(10)
16QAM每4比特映射一个符号，I路和Q路各用2比特，从{-3,-1,+1,+3}四个电平里选一个。平均符号功率等于(I电平均方 + Q电平均方)。I路的四个电平平方为9,1,1,9，均方为5；Q路同样为5，所以星座图总的平均功率为10。要让平均符号功率归一化为1，必须除以√10。这是16QAM仿真里最经典的错误来源——忘记归一化，星座图看起来没问题，但误码率曲线会整体偏移10log10(10)=10dB，肉眼可见地偏离理论值。
QAM16_LEVELS = np.array([-3, -1, 1, 3]) / np.sqrt(10) QAM16_GRAY = {0: -3, 1: -1, 3: 1, 2: 3} # Gray码到电平的映射 def qam16_map(bits): bits = bits[:4 * (len(bits) // 4)] symbols = np.zeros(len(bits) // 4, dtype=complex) for k in range(len(symbols)): b = bits[4 * k: 4 * k + 4] i_level = QAM16_GRAY[(b[0] << 1) | b[1]] q_level = QAM16_GRAY[(b[2] << 1) | b[3]] symbols[k] = i_level + 1j * q_level return symbols def qam16_demod(symbols): constellation = np.array( [i + 1j * q for i in QAM16_LEVELS for q in QAM16_LEVELS] ) # 最小欧氏距离硬判决 dist = np.abs(symbols[:, None] - constellation[None, :]) idx = np.argmin(dist, axis=1) bits = np.zeros(4 * len(symbols), dtype=int) rev_gray = {v: k for k, v in QAM16_GRAY.items()} for k, j in enumerate(idx): i_idx = (j >> 2) & 3 # 取高2位对应I路电平 q_idx = j & 3 # 取低2位对应Q路电平 i_bits = [rev_gray[QAM16_LEVELS[i_idx] * np.sqrt(10)], ] # 用整数电平索引反查Gray码更方便 level_int = int(np.round(QAM16_LEVELS[i_idx] * np.sqrt(10))) g_code = rev_gray[level_int] bits[4*k:4*k+2] = [(g_code >> 1) & 1, g_code & 1] level_int_q = int(np.round(QAM16_LEVELS[q_idx] * np.sqrt(10))) g_code_q = rev_gray[level_int_q] bits[4*k+2:4*k+4] = [(g_code_q >> 1) & 1, g_code_q & 1] return bits这段代码里QAM16_GRAY的键是2比特的整数值，0对应比特00，1对应01，3对应11，2对应10。这样Gray码相邻电平只差一比特。解调时先做全星座的最小距离判决，再反查Gray映射得到比特。因为星座点功率已经归一化，所以QAM16_LEVELS里存的就是归一化后的电平值，反查时乘回√10得到原始整数电平，方便查表。
实际使用中，如果你追求简洁，可以直接把16个星座点和对应比特做成一张表，解调时用np.argmin找最近星座点，然后再从表里取出比特。上面这个逐点循环的写法只是为了逻辑透明，性能不是重点——在几千个符号的仿真量级下，Python循环完全可接受。
3.3 软判决与硬判决：LLR计算的近似做法
硬判决只看接收符号落在哪个判决区域，直接输出比特；软判决则输出每个比特的似然信息，供后级的信道译码使用。在OFDM仿真里，如果后面接LDPC或Turbo译码，就需要软比特。一个常用近似是最大对数MAP，对第i个比特计算接收符号到两类星座点的距离差。
def qpsk_soft_llr(symbols, noise_var): # 简化的LLR：实部和虚部各自除以噪声方差的两倍 llr_i = 2 * np.sqrt(2) * symbols.real / noise_var llr_q = 2 * np.sqrt(2) * symbols.imag / noise_var return llr_i, llr_q这个近似对QPSK来说特别干净，因为QPSK的I路和Q路各自是独立BPSK。符号映射做了÷√2归一化后，幅度因子就是2×√2。如果noise_var估计不准，LLR的整体缩放会偏差，但这在硬判决仿真里无所谓，做软译码时才需要认真估计噪声方差。
3.4 星座图检查法：发射端最容易出错的三个表现
跑仿真时我几乎每次都会在发射端IFFT之前和接收端均衡之后各画一次星座图，这是定位问题最快的手段。星座图出错通常有三种表现：整体旋转、整体缩放、点集错位。整体旋转一般是映射表符号方向反了，比如Q路取反；整体缩放几乎都是归一化因子漏了或者FFT缩放不一致；点集错位则多半是Gray映射表写错或者数据子载波索引没对齐。
一个检查小技巧：发射端IFFT之前的星座图应该严格落在标准星座点上，任何偏差都是映射代码的问题；接收端均衡之后的星座图应该围绕标准星座点聚成簇，发散程度和SNR相关。如果接收端星座图是旋转但簇很紧，基本可以断定是信道估计里导频相位没校准。
4. OFDM核心收发：IFFT/FFT、循环前缀与频域均衡
4.1 为什么用IFFT/FFT而不是DCT：正交性来源
OFDM用IFFT做调制，用FFT做解调，本质上是因为复指数基函数在周期内正交。在接收端对时域采样做FFT，每个子载波上的信号能量被积分到对应的频点，子载波之间的干扰在理想同步下等于零。DCT虽然也有正交基，但它对应实数基函数，频谱利用率只有复指数的一半，而且处理复信道不自然，所以实际OFDM系统全用FFT。
我在第2章给的代码里用norm='ortho'，某种意义上就是要把IFFT当成一个正交变换来用。工程上很多实现为了省计算量用非规整FFT，结果接收端必须手动乘系数，这对仿真来说是纯粹的负担，没有必要。
4.2 循环前缀的作用与长度选择：从时延扩展到ISI
循环前缀解决两个问题：一是作为保护间隔，让前一符号的多径拖尾落在CP区间内，不污染当前符号的有用部分；二是把线性卷积变成循环卷积，使信道在频域呈现为逐子载波的乘积关系。第二点才是OFDM频域均衡的理论基础，也是CP被称为“循环”前缀的原因——复制符号尾部到头部，人为构造循环结构。
CP长度怎么选？如果信道最大时延扩展是L个采样间隔，那么CP必须大于L，否则去CP后仍有符号间干扰残留。CP也不是越长越好，它不携带信息，每加长一个采样点，信噪比开销就增加10log10((N_FFT+CP)/(N_FFT))。在第2章的参数里，CP取16，对应时延16/960kHz≈16.7微秒，能覆盖相当多室内场景。做仿真时你可以把CP分别设为4和16对比一下误码率，会看到CP不足时误码率曲线出现明显地板。
4.3 频域均衡：LS信道估计与一抽头均衡器代码
频域均衡的思想极其简单：信道在第k个子载波上的响应是H[k]，接收到的频域符号Y[k]=H[k]X[k]+W[k]，那么X的估计就是Y[k]/H[k]。问题变成怎么得到H[k]。最常见的是LS估计：在已知导频位置，H_hat = Y_pilot / X_pilot。
下面的代码实现了基于梳状导频的LS信道估计和线性插值均衡。导频每隔PILOT_STEP个子载波放一个，先估计导频位置的信道，再用np.interp插值得到所有子载波位置的信道响应。
def ofdm_rx_channel_estimation(x_cp): n_symbols = x_cp.shape[0] x_time = x_cp[:, N_CP:] Y = np.fft.fft(x_time, axis=1, norm='ortho') # LS信道估计：导频位置直接用接收值除以已知导频值 H_pilot = Y[:, pilot_idx] / PILOT_VALUE # 对每个OFDM符号，用线性插值补全所有子载波的信道响应 full_idx = np.arange(N_FFT) H_full = np.zeros_like(Y) for n in range(n_symbols): H_real = np.interp(full_idx, pilot_idx, H_pilot[n].real) H_imag = np.interp(full_idx, pilot_idx, H_pilot[n].imag) H_full[n] = H_real + 1j * H_imag # 一抽头均衡：逐子载波除以信道响应 Y_eq = Y / H_full return Y_eq[:, data_idx], H_full这个均衡器在AWGN信道下等于没做，因为H_full全接近1；在多径信道下，它能把频率选择性衰落压平。逻辑说明：H_pilot是一维复数数组，长度等于导频个数；插值时实部和虚部分开做，避免复数线性插值的不连续问题。参数说明：如果导频间隔太大，插值跟不上信道在频域的快速变化，均衡后剩余干扰会增大；信道时延扩展越大，导频就要越密。这正是2.4节参数联动关系在实际均衡中的体现。
4.4 参数设定：子载波间隔、采样率与带宽的关系
子载波间隔Δf、N_FFT和采样率fs三者的关系是fs = N_FFT × Δf。OFDM信号的总带宽约等于N_FFT×Δf，但因为边缘子载波通常留空作为保护带，实际占用带宽会小于这个值。在一个OFDM符号内，Δf越小，符号周期越长，对时延扩展的容忍度越高，但对频偏越敏感。
如果你只是做基带仿真，不打算把信号搬频到射频，那么记住一个原则：所有频率都用归一化角频率表示，采样率只影响时延扩展的采样点数和CP长度的换算。仿真时最省心的方式是直接以采样点为单位定义时延和CP长度，完全不涉及物理时间单位，代码更简洁，问题也更少。
5. 仿真避坑：OFDM+QPSK/16QAM最常见的5个翻车点
5.1 星座图整体旋转或幅度不对：忘了归一化或忘了能量调整
现象是星座图在接收端画出来，整体绕原点转了45°，或者幅度比标准星座点大了一圈。原因通常是QPSK映射表里I/Q顺序反了，或者发送端映射时忘了除以√2，又或者是IFFT用了默认缩放而接收端FFT没有对应调整。解决方法是先检查发射端IFFT之前的星座图：如果这里就不对，说明映射代码有问题；如果这里对而接收端不对，检查FFT缩放因子和信道估计的相位。用norm='ortho'之后，这个坑会少掉一半。
5.2 误码率曲线在高信噪比处掉头：同步误差或CP不够长
误码率随着EbN0增加先降后平，出现地板效应，是一个典型症状。先怀疑CP长度：如果信道时延扩展大于CP，那么高信噪比下残留的符号间干扰成为主导，误码率不可能再下降。再看同步：如果接收端ОФDM符号起始点偏了N_CP个采样甚至偏了任意采样，FFT窗口内混入相邻符号的数据，同样会出现地板。解决方式是先把信道设成纯AWGN，如果地板消失，说明问题在信道配置或CP；如果还在，说明收发帧结构有问题。
5.3 EbN0和SNR换算错误：归一化因子少算一次
现象是仿真误码率曲线和理论曲线平行但整体偏移几个dB。原因几乎都是EbN0到噪声功率的换算没有考虑CP开销和导频开销。我在做这套仿真时，噪声方差按下面的公式设置，它显式包含了N_FFT、N_CP和N_DATA三个因子：
def awgn_noise_var(ebno_db, mod_order): ebno_lin = 10 ** (ebno_db / 10) bits_per_sym = int(np.log2(mod_order)) # 频谱效率：数据比特数 / 实际发送采样数 spectral_eff = N_DATA * bits_per_sym / (N_FFT + N_CP) noise_var = 1.0 / (ebno_lin * spectral_eff) return noise_var如果漏掉CP开销系数，16QAM的曲线会偏移10log10(80/64)≈0.97dB；如果漏掉导频开销，QPSK曲线会偏移10log10(96/48)≈3dB。这个偏移和高斯白噪声下理论曲线的形状完全一样，所以特别容易漏判。建议画曲线时把带开销换算和不带开销换算的两条仿真曲线放在同一张图里，对比一目了然。
5.4 导频位置和解映射索引错位：打点符号与取点符号不一致
现象是误码率在小信噪比时正常，高信噪比时突然变差，或者星座图上有少数点严重偏离。这通常是导频索引和数据索引在发射端与接收端定义不一致。比如发射端用pilot_idx=np.arange(0,64,4)得到0,4,…,60，接收端却手写了一个数组，漏掉最后一个导频位置，导致数据子载波索引整体偏移。解决方式是让收发两端共享同一份pilot_idx和data_idx，不要在两处各写一遍。我一般在代码顶层只定义一次，后面所有函数都引用同一份变量。
5.5 16QAM在高信噪比出现地板效应：信道估计精度不足
跟5.2不同的是，这个地板出现在16QAM下而QPSK正常。原因大概率是导频密度不够或者插值方式太粗糙，导致信道估计误差在高阶调制下成为主导干扰。16QAM的星座点间距比QPSK密，对残余误差更敏感。解决方法是加密导频，把PILOT_STEP从4改成2，或者改用二阶插值代替线性插值。这个坑也是我在多径信道演示里踩过最多次的，一个简单的两径信道，导频间隔8时16QAM地板明显，间隔4时基本消失。
6. 用误码率曲线验证整套仿真：从“能跑”到“可信”
6.1 理论BER曲线与仿真曲线的对照方法
验证仿真程序是否正确，唯一的可靠手段是让仿真误码率曲线和理论曲线对得上。QPSK在AWGN信道下的误比特率公式是BER = 0.5×erfc(√EbN0)。16QAM没有精确闭式解，但用Gray映射时可以近似为BER ≈ 0.75×0.5×erfc(√(0.8×EbN0))，在误码率低于1e-2时误差已经很小。这两条理论曲线可以直接用scipy.special.erfc画出来。
6.2 完整仿真循环与耗时优化
完整仿真循环分三层：外层遍历EbN0点，中层生成多个OFDM符号并叠加噪声，内层做收发和比特统计。每个EbN0点跑的符号数决定曲线的稳定程度，一般每个点累计至少2000个数据符号，误码率才能画到1e-3量级。下面给出核心循环代码：
def run_ber_simulation(ebno_db_list, mod_order, n_symbols_per_point=2000): ber_list = [] for ebno_db in ebno_db_list: noise_var = awgn_noise_var(ebno_db, mod_order) total_bit_errors = 0 total_bits = 0 n_blocks = 20 block_symbols = n_symbols_per_point // n_blocks for _ in range(n_blocks): if mod_order == 2: bits = np.random.randint(0, 2, block_symbols * N_DATA * 2) tx = ofdm_tx(bits) noise = np.sqrt(noise_var / 2) * ( np.random.randn(*tx.shape) + 1j * np.random.randn(*tx.shape) ) rx_symbols = ofdm_rx_awgn(tx + noise) rx_bits = qpsk_demod(rx_symbols.flatten()) else: bits = np.random.randint(0, 2, block_symbols * N_DATA * 4) symbols = qam16_map(bits) # 组帧后按ofdm_tx流程发送，这里省略重复的组帧代码 # 噪声和接收流程与QPSK分支一致 rx_bits = qam16_demod(rx_symbols.flatten()) total_bit_errors += np.sum(rx_bits != bits[:len(rx_bits)]) total_bits += len(rx_bits) ber_list.append(total_bit_errors / max(total_bits, 1)) return ber_list这个循环把每个EbN0点拆成20个数据块来跑，一方面避免一次生成超大规模数组占用内存，另一方面方便以后做置信区间统计。参数说明：noise_var直接来自5.3节的换算函数，噪声复高斯生成时实部虚部分别用sqrt(noise_var/2)，总方差才等于noise_var。mod_order=2时每符号2比特，mod_order=4时每符号4比特，比特总数按块内符号数和数据子载波数计算。
一段值得记住的个人教训：我最初做这套仿真时，偷懒只画了QPSK的BER曲线，16QAM仅仅看了星座图觉得“差不多”，后来给别人演示多径场景才发现16QAM的信道估计问题。从那以后，我养成了一个习惯——任何调制方式改动，都强制跑一次0到12dB的BER对比曲线，仿真里多花十分钟，能省掉事后排错的半天。希望帮到你。
本文还有配套的精品资源，点击获取
