简介在无线通信与信号处理研究中这套MATLAB代码包聚焦BPSK循环谱估计仿真面向通信工程专业学生、算法工程师以及需要深入分析信号循环平稳特性的研究人员。压缩包共3个MATLAB脚本整体仅2KB代码量精简且步骤清晰分别为噪声生成、BPSK/MSK基带信号构建、以及自相关计算与循环谱三维绘图便于观察信号循环谱图这种划分也方便按模块学习自相关、傅里叶变换与三维可视化的关联。目前已有207人学习/下载。借助该脚本用户可以快速掌握循环谱估计的基本流程理解BPSK信号在多径与时变信道下的循环特征为后续调制方式识别、盲信道估计和频谱感知提供可复现的仿真基础并可在不同信噪比条件下对比谱特征变化。由于代码紧凑也适合二次修改用于教学演示或科研验证。1. BPSK循环谱估计到底是干什么的比功率谱多出一维特征把一段 BPSK 信号丢进普通的 FFT 功率谱你只能看到载频附近一根凸起的谱线码速率、符号周期这类信息几乎全部藏在相位跳变的细节里。换成 BPSK 循环谱估计事情就完全不一样了二维谱图上会在循环频率 α ±RbRb 为码速率的位置出现明显的峰值码速率直接“照”出来。这就是这个 code 方案的价值——用可运行的实现把 BPSK 信号从一维频谱扩展到循环频率与谱频率构成的二维平面为调制识别、盲参数估计、频谱感知这类任务提供一张信息量高得多的“指纹图”。这篇文章面向需要实际跑通循环谱估计的从业者你可能在做信号检测、干扰识别或是软件无线电里的调制分类。我从原理讲到代码再到参数设置和几个高频踩坑点争取你照着复现完就能在自己的数据集上用起来。2. 循环谱怎么把 BPSK 的码速率“照”出来原理与第一版代码2.1 功率谱看不出码速率BPSK 的循环平稳性从哪来BPSK 信号有一个被功率谱“隐藏”掉的性质它是循环平稳的。所谓循环平稳指的是信号的统计特性尤其是自相关函数随时间周期性变化。对 BPSK 而言这个周期性来自符号序列的重复结构——每个码元持续时间 Tb 内相位保持不变到了下一个 Tb 可能跳变这种以 Tb 为周期的统计规律会让自相关函数在时延域之外还多出一个“循环频率”维度。循环自相关函数的定义是R_x^α(τ) ⟨x(t τ/2) · x*(t − τ/2) · e^{−j2παt}⟩_t其中 α 就是循环频率。对 R_x^α(τ) 再做一次关于 τ 的傅里叶变换就得到循环谱密度函数 S_x^α(f)。普通功率谱只是 α 0 时的一个切片而完整的循环谱把每个循环频率 α 下的谱密度都算出来BPSK 的很多特征参数就从这个二维分布里暴露出来了。对 BPSK 信号理论上最显著的特征峰出现在三处α 0 处对应常规功率谱α ±Rb 处对应码速率峰值出现在 f 0 附近α ±2fc 附近也有载频相关的循环特征。实际做中频实信号仿真时由于双边带镜像的存在谱图会比复基带信号多出一些对称副本但 α ±Rb 的峰始终是最稳定、最容易提取的调制指纹。2.2 先造一段能复现的 BPSK 信号仿真信号生成器没有公共数据集的情况下我一般先自制一段 BPSK 中频信号来做算法验证。生成器的核心是把随机 bit 映射成 ±1 的矩形脉冲序列再乘以载波 cos(2πfc·t)。这里的常见做法是直接对符号序列做零阶保持上采样前提是采样率与码速率成整数倍关系import numpy as np def generate_bpsk(fc, Rb, fs, duration, snr_dbNone): 生成中频 BPSK 实信号。 fc: 载波频率(Hz) Rb: 码速率(bit/s)要求 fs/Rb 为整数 fs: 采样率(Hz) duration: 信号时长(s) snr_db: 信噪比(dB)None 表示不加噪声 t np.arange(0, duration, 1.0 / fs) n_symbols int(duration * Rb) bits np.random.randint(0, 2, n_symbols) symbols 2 * bits - 1 # 0 - -1, 1 - 1 sps int(fs / Rb) # 每符号采样点数 baseband np.repeat(symbols, sps) # 零阶保持上采样 baseband baseband[:len(t)] # 截断到与 t 等长 x baseband * np.cos(2 * np.pi * fc * t) if snr_db is not None: signal_power np.mean(x ** 2) noise_power signal_power / (10 ** (snr_db / 10.0)) noise np.sqrt(noise_power) * np.random.randn(len(x)) x x noise return x, t这段代码里最容易踩坑的是sps int(fs / Rb)如果 fs 不能被 Rb 整除取整后实际符号长度会漂移仿真出来的信号码速率和设定值不一致后面循环谱上峰的位置也会对不上。我在第 5 章专门讲这个问题。另一个需要注意的点是snr_db的核算方式这里以信号功率为基准折算噪声功率和直接给噪声幅度是两个概念。验证生成器是否正常可以先看时域波形是否在相位跳变点出现幅度过零或者直接拿功率谱看载频处是否有一根谱线。2.3 把循环谱定义落成第一版代码循环自相关直接估计先给一个“教科书直译版”的循环谱估计实现思路完全从定义出发对每个循环频率 α按公式计算循环自相关再对 τ 做傅里叶变换得到 S_x^α(f)。这个版本慢但逻辑透明适合验证你的信号模型是否正确def cyclic_spectrum_direct(x, fs, alpha_res, tau_max): 基于循环自相关定义的直接估计。 alpha_res: 循环频率扫描步长(Hz) tau_max: 最大时延(采样点)决定频率分辨率 n len(x) alpha_axis np.arange(-fs/2, fs/2, alpha_res) tau_axis np.arange(-tau_max, tau_max 1) S np.zeros((len(alpha_axis), len(tau_axis)), dtypecomplex) for ai, alpha in enumerate(alpha_axis): for ti, tau in enumerate(tau_axis): t_idx np.arange(n - abs(tau)) # 移位的两个采样序列注意负时延的方向 if tau 0: prod x[t_idx tau] * np.conj(x[t_idx]) else: prod x[t_idx] * np.conj(x[t_idx - tau]) R np.mean(prod * np.exp(-2j * np.pi * alpha * t_idx / fs)) S[ai, ti] R # 对时延维做 FFT 得到该 alpha 下的谱切片 S[ai, :] np.fft.fftshift(np.fft.fft(S[ai, :])) return S, alpha_axis, tau_axis注意这里为了演示做了两层循环实际上数据量稍大就跑不动了alpha_axis 如果有 200 个点、tau 有 256 个点就要执行 51200 次内层计算每次还涉及数组乘法。这个版本的价值在于对照公式查错——如果你的信号生成器有 bug拿这个直译版先算一个 α 0 的切片和普通功率谱对比能立刻定位问题。等确认模型无误后再切换到第 3 章的频域相关实现那个版本快得多适合做参数扫描和批量估计。3. 从定义到可运行代码基于 STFT 频域相关法的循环谱估计实现3.1 为什么从循环自相关切换到频域相关法直接估计版代码能跑通但没法在实际项目里用。比如 fs 10 MHz、信号时长 5 ms 的数据时域自相关法算一次要几分钟而频域相关法只需要把信号分帧做 STFT再在帧维度上做统计平均计算量低了一个量级。频域相关法的理论依据是循环谱密度的频域形式S_x^α(f) 等于频谱分量 X(f α/2) 与 X*(f − α/2) 的时间平均。换句话说把信号切成一帧一帧做 FFT然后把相隔 α 的两个频点相乘、对所有帧求平均就得到了该循环频率下的谱值。这种实现方式对 BPSK 特别友好。因为 BPSK 的码速率峰出现在 α ±Rb、f ≈ 0 的位置频域相关法天然把频率对 (f α/2, f − α/2) 组织成了一个可索引的二维矩阵峰值定位非常直接。相比时域平滑法如 FAM 算法频域相关法的循环频率轴分辨率可以直接由帧长控制不需要额外的低通滤波器设计代码实现的出错率低很多——这正是我把它作为主推方案的原因。3.2 基于 STFT 的频域相关实现可直接复现的完整代码下面这段代码就是标题里的 code 核心。它把上一步的循环自相关直译版替换为基于 STFT 的频域相关估计输入一段 BPSK 信号输出循环谱矩阵 S、频率轴和循环频率轴import numpy as np def bpsk_cyclic_spectrum(x, fs, nfft1024, overlap0.5, max_alpha_hzNone): 基于 STFT 频域相关的循环谱估计。 x: 输入信号实信号或复基带信号均可 fs: 采样率 nfft: STFT 的 FFT 点数决定频率分辨率 df fs/nfft overlap: 相邻帧重叠率0~1 max_alpha_hz: 循环频率扫描范围默认 fs/2 返回: S: 循环谱矩阵形状 (len(alpha_axis), nfft) f_axis: 频率轴长度 nfft alpha_axis: 循环频率轴 if max_alpha_hz is None: max_alpha_hz fs / 2.0 # 分帧、加窗、FFT win np.hanning(nfft) step int(nfft * (1 - overlap)) n_frames max(1, 1 (len(x) - nfft) // step) X np.empty((n_frames, nfft), dtypecomplex) for i in range(n_frames): seg x[i * step : i * step nfft] X[i, :] np.fft.fft(seg * win, nfft) # 频率轴搬到零中频方便观察 f0 附近的码速率峰 X np.fft.fftshift(X, axes1) df fs / nfft max_m int(np.floor(max_alpha_hz / (2.0 * df))) alpha_axis 2.0 * df * np.arange(-max_m, max_m 1) n_alpha len(alpha_axis) S np.zeros((n_alpha, nfft), dtypecomplex) for mi, m in enumerate(range(-max_m, max_m 1)): # 循环频率 alpha 2*m*df对应频谱分量间距为 2m 个频点 k_start abs(m) k_stop nfft - abs(m) k np.arange(k_start, k_stop) # 对每帧做 X[km] * conj(X[k-m])再对帧维平均 S[mi, k_start:k_stop] np.mean( X[:, k m] * np.conj(X[:, k - m]), axis0 ) f_axis np.linspace(-fs / 2, fs / 2, nfft) return S, f_axis, alpha_axis逻辑说明分帧后的 X 矩阵形状是帧数, nfft每一行是某一帧的频谱。外层循环的 m 控制循环频率间隔当 m 0 时 S 退化为常规功率谱当 m ≠ 0 时代码取频点 km 和 k−m 的共轭乘积在帧维度上做平均。这里的平均是循环谱估计的核心——它把噪声项在统计上压下去而 BPSK 的循环平稳特征项由于各帧相位一致会随帧数增加而线性增长。帧数越多谱图越干净代价是计算时间变长。参数说明nfft 1024 是默认值对应 fs 10 MHz 时频率分辨率约 9.77 kHzmax_alpha_hz 默认扫到 ±5 MHz但对码速率 1 Mbit/s 的信号来说扫到 ±2 MHz 就足够了扫太宽浪费时间。overlap 0.5 表示相邻帧重叠一半增加重叠率可以提升帧数、改善统计平均效果但计算量也随之上升。如果你的信号是复基带形式直接把 x 换成复数数组即可频率轴会相应变为 −fs/2 到 fs/2 的单边镜像结构代码不需要任何改动。3.3 从二维矩阵里读出 BPSK 的 α±Rb 特征峰跑完上面的函数后S 是一个二维复矩阵。直接看复数不方便通常取幅度谱做可视化或峰值搜索。常见做法是取S_mag np.abs(S)再用imshow画出来注意把 α 轴设为纵轴、f 轴设为横轴动态范围用对数刻度压缩一下否则 α0 处的强峰会淹没 α±Rb 的弱峰。下面是粗略的谱峰定位代码def find_bpsk_peaks(S, f_axis, alpha_axis, f_target0.0, exclude_alpha_zeroTrue): 在指定频率切片上搜索循环频率峰值。 f_idx np.argmin(np.abs(f_axis - f_target)) profile np.abs(S[:, f_idx]) ap alpha_axis.copy() if exclude_alpha_zero: # 排除 α0 附近的常规功率谱成分 guard np.max(np.abs(np.diff(alpha_axis))) * 2 mask np.abs(alpha_axis) guard profile profile[mask] ap ap[mask] peak_order np.argsort(profile)[::-1][:4] # 取最大的 4 个峰 return ap[peak_order], profile[peak_order]对第 2 节生成的 BPSK 信号fc2 MHz、Rb1 Mbit/s、fs10 MHz、时长 5 ms、信噪比 10 dB运行后返回的前两个峰应当非常接近 α ±1 MHz。如果峰值偏离超过 10%优先检查信号生成器里的 sps 是否精确等于 fs/Rb。我踩过的坑是用这个函数直接找峰时不排除 α0 的直流分量结果前四个峰全落在 α0 附近误导排查方向。4. 决定谱图质量的三个参数nfft、循环频率扫描范围与帧重叠率4.1 nfft 怎么定频率分辨率和循环频率分辨率同时受它控制nfft 是循环谱估计里牵一发动全身的参数。它同时决定两件事频率分辨率 Δf fs/nfft以及循环频率轴的最小间隔 Δα 2Δf。看第 3 节的代码就知道循环频率轴是用2.0 * df * np.arange(...)构造的也就是说 α 轴的网格密度直接由 nfft 决定。nfft 越大Δf 越小α 轴越密码速率峰可以被定位得更准但分帧后每帧数据更长在总数据量不变的前提下帧数会变少统计平均质量下降。这是一个需要权衡的矛盾。我在实际项目中一般按这个流程定 nfft先根据你对码速率的先验估计让 α 分辨率至少小于预期码速率的 1/20。比如 Rb 约 1 Mbit/s那么 Δα 至少要 ≤ 50 kHz即 2Δf ≤ 50 kHz推得 nfft ≥ 2fs/50000fs10 MHz 时 nfft ≥ 400取 512 或 1024 比较合理。如果你是完全盲估计没有先验信息就先用 512 快速扫一遍确认峰的大致位置后再加大 nfft 精估。4.2 max_alpha_hz扫太宽看不到峰扫太窄丢特征max_alpha_hz 参数控制循环频率轴的扫描范围。它的上限是 fs/2但绝大多数情况下不需要扫那么宽。BPSK 的特征峰集中在 α0 和 α±Rb 附近如果你知道 Rb 的范围把 max_alpha_hz 设为 2~3 倍的 Rb 上限就足够。扫太宽有两个坏处一是计算量大增循环次数和 α 轴长度线性相关二是谱图动态范围被拉大弱峰更容易被噪声淹没。扫太窄的坏处更隐蔽——你可能会把 α±2fc 附近的镜像峰误判成码速率峰尤其在载频比较低的时候。一个实用的检查方法先用 max_alpha_hz fs/4 跑一遍打印出 α 轴上峰值出现的位置如果最高峰落在 α0 之外的区域再把范围收窄到该峰的 2 倍附近精扫。这样做虽然多跑一次但在盲估计场景下比盲目猜测参数靠谱得多。另外扫宽不变时增大 nfft 会自动让 α 轴变密谱图会更细腻但扫描范围不变。4.3 重叠率、窗函数与数据长度计算量和谱图干净度的取舍overlap 参数控制相邻帧的重叠程度。默认 0.5 是个不错的起点如果你发现谱图噪声太大、峰不够锐利可以提高到 0.75 甚至 0.875。重叠率每提高一倍帧数约增加一倍统计平均的帧数更多谱峰更锐利但计算时间也近似翻倍。超过 0.875 后收益迅速递减我一般不会设更高因为窗函数的加窗效应会让连续帧之间的信息高度冗余。窗函数方面代码里默认用了 Hanning 窗。它只影响每一帧 FFT 的频谱泄漏对循环谱峰位置没有系统性影响。如果你处理的是单音干扰较强的场景可以换 Blackman 窗获得更高的旁瓣抑制但主瓣会变宽频率分辨率略降。在需要精确估计峰值的场景我反而推荐矩形窗即不加窗因为它的主瓣最窄峰定位误差最小代价是旁瓣较高。噪声环境恶劣时用 Hanning信噪比足够高时用矩形窗这是我的默认规则。数据长度也是一个容易被忽略的维度。总时长 T 决定了帧数上限帧数 ≈ (T·fs − nfft)/(nfft·(1−overlap))。帧数少于 30 时谱峰统计起伏会非常明显看起来像是随机噪声里凸起几个毛刺。建议在参数扫描前先确认帧数如果帧数太少优先增大信号时长而不是减小 nfft。下表给出一组经验参考值fs10 MHz、nfft1024、overlap0.5 的条件下信号时长帧数谱峰质量适用场景1 ms约 19峰可见但毛刺多快速侦察、粗扫5 ms约 96峰明显可用于估计常规参数估计20 ms约 389峰锐利低信噪比可用弱信号检测、精估5. BPSK 循环谱估计避坑指南5 个翻车现场5.1 谱图一片噪点找不到 α±Rb 的峰现象跑完代码后画出来的谱图全是细碎噪点完全没有明显的峰状结构。原因通常是信噪比太低或者帧数不够导致统计平均没有收敛。BPSK 循环谱在低信噪比下的噪声基底比较高α±Rb 处的信号特征峰如果没有足够的帧数积累会被起伏淹没。解决先增大信号时长到 10 ms 以上确认帧数超过 100再检查是否加了噪声如果加了把信噪比先调到 20 dB 做算法连通性测试。算法本身没问题再逐步降信噪比观察极限。我一般会在调试时加一个打印语句输出最大帧数和平均帧能量方便快速定位是数据短了还是真没信号。5.2 α 轴糊成一团码速率峰被抹平现象谱图能隐约看到能量分布但沿 α 轴方向峰与峰之间没有明显间隔整个图“糊”在一起。原因nfft 太小导致 Δα2Δf 粗于码速率αRb 和 α−Rb 的峰之间只有一两个采样点根本分不开。解决把 nfft 从 256 加到 1024 或 2048Δα 会成比例变细。同时确认 max_alpha_hz 没有把 α0 附近的能量展得太宽适当收窄扫描范围也能让谱图对比度更高。这个坑在 fs 较高但 Rb 较低时特别常见比如 fs100 MHz、Rb1 Mbit/snfft 用 256 时 Δα≈781 kHz和 2Rb 差不多分辨率完全不够。5.3 α0 处直流量巨大把细节全淹了现象谱图上 α0 这一条线异常明亮其他区域即使有峰也看不见。原因α0 对应常规功率谱载波能量和信号能量全部集中在这里动态范围比 α±Rb 处的循环谱特征高出几十 dB。画图时如果线性刻度小峰全部被压平。解决一是绘图时用10*np.log10(S_mag eps)取对数二是对显示范围做截断比如把色标上限设在最大值的 60%三是在找峰逻辑里显式排除 α0 附近的区域参考第 3 节 find_bpsk_peaks 里的 mask 逻辑——严谨的处理是用常规功率谱的旁瓣宽度作为保护间隔而不是随意取一个固定值。5.4 估出的码速率总是实际值的一半或两倍现象峰值搜索找到了很锐利的峰但峰位置对应的 α 值和设定的 Rb 差一倍。原因循环谱的天然对称性会在某些条件下引入镜像峰实信号 BPSK 的循环谱在 α±Rb 处有峰同时在 α±2fc±Rb 处也会出现镜像峰当 fc 较低或频谱重叠时这些镜像峰可能比真实峰更突出。另一个常见原因是 nfft 太小时 α 网格太粗峰值落在两个网格点之间抛物线插值缺失时直接取了邻近整数点误差刚好可以到半个网格。解决先用理论位置校验峰数量BPSK 应当看到对称的一对峰再用抛物线插值精估最后检查 fc 选择是否过低——通常让 fc 大于 3 倍 Rb 可以有效分离镜像峰。5.5 仿真信号本身的坑sps 非整数导致仿真失真现象生成信号脚本跑得很顺利循环谱也画出来了但峰位置和预期 Rb 系统性偏差 3%~10%。原因第 2 节代码里sps int(fs / Rb)做了取整当 fs 不能被 Rb 整除时实际的符号采样点数偏离理想值等效码速率就变了。比如 fs10 MHz、Rb1.3 Mbit/ssps7 取整实际 Rb≈1.429 Mbit/s偏差接近 10%。解决要么选择使 fs/Rb 为整数的参数组合最省事要么修改生成器对符号序列做时分插值让每个符号的采样点数可以是小数。作为快速修复我在工程里常用前者仿真验证算法足够做半实物联调时再用任意码速率的生成方案。另外要提醒这个误差在功率谱里几乎看不出来只在循环谱这类对时序敏感的参数估计里暴露属于最容易漏掉的隐藏 bug。6. 用循环谱自动估码速率和频偏峰值搜索与抛物线插值拿到干净的循环谱之后下一步就是把“用眼睛看峰”变成“用代码读数”。我常用的技巧是在 f0 附近取一个切片沿 α 轴搜索峰值再用抛物线插值把估计精度推到亚网格级别。频偏的估计则是看 α0 处 f 轴的峰值偏移量两者结合就能同时输出码速率和载波频偏两个关键参数。下面给出一个可直接套用的峰值提取片段def refine_alpha_peak(profile, alpha_axis, peak_idx): 用三点抛物线插值细化峰值位置。 d_alpha alpha_axis[1] - alpha_axis[0] if 0 peak_idx len(profile) - 1: denom (profile[peak_idx-1] - 2*profile[peak_idx] profile[peak_idx1]) if abs(denom) 1e-12: delta 0.5 * (profile[peak_idx-1] - profile[peak_idx1]) / denom else: delta 0.0 return alpha_axis[peak_idx] delta * d_alpha return alpha_axis[peak_idx] # 实际用法取 f0 切片排除 α0 后找最大峰 f_idx np.argmin(np.abs(f_axis - 0.0)) profile np.abs(S[:, f_idx]) candidate np.argmax(profile[np.abs(alpha_axis) 2e4]) # 排除 α≈0 alpha_est refine_alpha_peak(profile, alpha_axis, candidate)这里的抛物线插值用的是峰值点及其左右两点的二次拟合delta 表示偏离中心格点的比例。三点抛物线只在主峰形状接近二次曲线时精度高如果谱峰旁瓣大可以先对 profile 做一次滑动平均再搜索能明显改善稳定性。一个经验是信噪比 10 dB 以上、帧数 100 左右时插值后的码速率估计误差通常能控制在 1% 以内低信噪比下不建议强行精估先做帧累积或跳频积累更实际。频偏估计是类似的逻辑在 α0 切片上找 f 轴的最大峰峰位对应载波频率估计如果信号发射时知道标称载频两者相减就是频偏。需要提醒的是频偏估计精度受 nfft 限制Δf 是多少估计误差就近似是多少想提高就得加大 nfft。我在做干扰识别的时候会把码速率、频偏、峰宽三个量合成一个特征向量喂给分类器这比单纯用能量检测可靠得多。我的习惯是在每次实验前先固定跑一遍参考信号。把同一段干净的 BPSK 信号存成 .npy 文件作为回归测试基准任何参数改动后先拿它验证峰位不变再上真实数据。这个小习惯帮我避免了不少“参数调了半天结果是生成器坏了”的尴尬。整套流程从生成信号到循环谱计算再到参数提取核心代码不到一百行却覆盖了调制识别和盲估计里最常用的一环。希望帮到你。本文还有配套的精品资源点击获取