用Matlab实现DSSS盲估计：载频、码速率与码周期参数的工程估计方法

简介：这套Matlab实现的直接序列扩频信号参数盲估计系统，面向通信工程专业学生与信号处理研究人员，解决DSSS/BPSK信号在未知发射端参数下的载频、码速率与码周期自动估计问题，适用于非合作信号分析、认知无线电、频谱监测等典型场景。包内共35个文件，以24个.m源程序为核心，分别实现载频估计、码速率估计、码周期估计、BPSK解调、m/平衡Gold序列生成与信号同步等环节；另有6个.asv自动备份、4个.fda滤波器设计文件和1个.mexw32编译文件，方便对比原始设计和加速执行。压缩包整体仅444KB，轻量易用，便于快速部署与二次开发。当前已有651人学习。通过该资源可获得一套完整的盲估计处理流程，覆盖周期图与自相关载频检测、滑动窗码速率估计、互相关码周期搜索等经典方法，并配有相对误差计算与信号波形绘图模块，可直接运行或修改，帮助深入理解扩频通信与非合作信号处理。

1. Matlab 做直接序列扩频盲估计，先分清载频、码速率和码周期各自在哪

非合作场景里拿到一段信号，往往是只知道采样率和大致频段，不知道调制参数。对直接序列扩频（DSSS）信号做参数盲估计，目标就是在没有导频、不知道扩频码的条件下，把载频 f_c、码速率 R_c 和码周期 T_b 三个量估出来。后续解扩、解调、信息恢复都要以这三个参数为前提。

这里有个常被新手低估的结论：BPSK 调制反而降低了盲估计难度。信号平方之后，信息位和扩频码都变成 1，2 倍载频处会留下一根干净的谱线，载频因此成了三个参数里最容易估的；码速率藏在谱线展宽的一阶零点附近，码周期藏在延迟相乘输出的自相关副峰里，这两个才是重头戏。

这篇按“信号建模 → 生成截获数据 → 估计载频 → 估计码速率与码周期 → 闭环验证”的顺序，给出一套能在 Matlab 里直接跑的流程。做无线电监测、认知无线电和通信抗截获方向的工程师或研究生可以照搬参数，做课程设计的人也能一步步对照着改。

2. 建立 BPSK 直扩信号模型：采样率、码片速率和码周期如何设定

2.1 短码直扩信号的数学模型与参数换算关系

盲估计的第一步不是写估计算法，而是先把信号模型固定下来。常见的短码直扩 BPSK 信号写作：s(t)=A·d(t)·c(t)·cos(2πf_c t+φ)。d(t) 是 ±1 的信息序列，每个信息位持续一个符号周期 T_b；c(t) 是 ±1 的扩频码序列，由 N 个码片组成，周期恰好等于符号周期。也就是说一个信息位期间扩频码完整重复一遍，这是短码体制的定义，也是码周期能被盲估出来的前提。

四个参数的换算关系是：码速率 R_c=N/T_b，信息速率 R_b=1/T_b=R_c/N。盲估计端口的任务是得到 (f_c, R_c, T_b) 三元组，扩频码长度 N 不需要单独估，取 round(T_b·R_c) 即可。仿真时还要额外定两个量：采样率 fs 和每码片采样点数 sps=fs/R_c。fs 要同时满足 fs>2f_c（保证平方谱的 2f_c 谱线不混叠）和 sps≥8（保证延迟相乘的延迟量落在半个码片之内），两者冲突时先带通滤波再降采样。

2.2 用 Matlab 生成可盲估计的仿真信号（m 序列扩频）

代码把“发射端已知参数”和“接收端盲估计输入”分开写，便于后面做误差对比。扩频码用 5 级 m 序列，长度 31，自带周期特性，比随机序列更适合观察码周期副峰。

fs = 100e6; % 采样率 100 MHz fc = 20e6; % 载频 20 MHz Rc = 5e6; % 码片速率 5 Mchip/s N = 31; % 扩频码长度，短码：一符号一周期 Nb = 512; % 信息符号个数 SNRdB = 5; % 接收信噪比 % 5 级 LFSR 生成 31 位 m 序列，本原多项式 x^5+x^3+1 reg = ones(1,5); pn = zeros(N,1); for n = 1:N pn(n) = reg(5); fb = xor(reg(5), reg(3)); % 抽头 5、3 反馈 reg = [fb, reg(1:4)]; end code = 1 - 2*pn; % 双极性：0->+1, 1->-1 % 信息位扩频：每行一个符号，与扩频码逐码片相乘 data = randi([0 1], Nb, 1); dbpsk = 1 - 2*data; chips = reshape(dbpsk * code.', [], 1); % Nb*N 行码片级基带 % 矩形成型上采样：每码片 sps 个采样 sps = fs / Rc; bb = reshape(repmat(chips.', sps, 1), [], 1); tt = (0:length(bb)-1).'/fs; s = bb .* cos(2*pi*fc*tt); % BPSK 调制 % 加性高斯白噪声 Ps = mean(s.^2); Pn = Ps / 10^(SNRdB/10); r = s + sqrt(Pn)*randn(size(s)); % 盲估计输入

dbpsk * code.是矩阵外积，得到 Nb×N 矩阵，第 k 行正好是第 k 个符号的扩频结果，reshape 成列向量就是码片级基带序列，比循环快一个量级。repmat(chips.', sps, 1)把每个码片重复 sps 次，得到矩形脉冲成型；如果截获信号经过升余弦滤波，谱零点不会正好落在整数倍码速率上，盲估计前建议先用带通滤波器压掉旁瓣。

噪声按 Pn=Ps/10^(SNRdB/10) 换算，这里的 SNR 指整个信号带宽内的信噪比，不是解扩后的处理增益。整段信号长度 NbNsps=317440 个采样，约 3.17 ms，包含 16.5 个码周期。码周期估计要求观测覆盖至少 5 个码周期，这段留了余量。

2.3 估计前先自检：频段位置与观测时长

盲估计算法对输入有两条最低要求：2f_c 谱线落在 0~fs/2 内，观测时长至少 5 个码周期。拿到截获数据先画一次功率谱，确认谱峰大致位置再进盲估。

3. 载频盲估计：平方谱谱线定位在 Matlab 里的实现与精度修正

3.1 平方谱为什么在 2 倍载频处形成谱线

对接收信号直接平方：r²(t)=A²b²(t)cos²(2πf_c t+φ)，其中 b(t)=d(t)c(t)∈{±1}，b²(t)≡1。展开后 r²(t)=(A²/2)[1+cos(4πf_c t+2φ)] 再加噪声项。数据位和扩频码在一平方之后全部消失，只剩频率 2f_c 的余弦，FFT 后谱峰频率除以 2 就是载频。

这个性质只对 BPSK 成立。QPSK 信号符号取 ±1、±j，平方后符号变 ±1 仍带调制，要四次方才能消掉；DSSS 的码片和信息位都是实 ±1，平方后自然消调。所以“平方谱线”是 BPSK 直扩盲估计里性价比最高的一步。另一个前提是 fs>2f_c，载频接近 fs/2 时谱线折叠到低频段，估计值出错且不易察觉。

3.2 平方谱找峰、抛物线插值与谱线突出度判别的 Matlab 代码

% 载频盲估计：平方谱谱峰定位 Nfft = 2^nextpow2(length(r)); S = abs(fft(r.^2, Nfft)).^2; % 平方谱 f = (0:Nfft-1)/Nfft*fs; % 单边频率轴 % 搜索区间：2fc 粗先验 ±10 MHz，实践中由接收机扫描带宽给出 flo = 2*fc - 10e6; fhi = 2*fc + 10e6; ig = find(f >= flo & f <= fhi); [~, kk] = max(S(ig)); k = ig(1) + kk - 1; % 谱峰序号 % 抛物线插值修正峰位置 if k > 2 && k < Nfft-1 a = S(k-1); b = S(k); c = S(k+1); delta = 0.5*(a - c) / (a - 2*b + c); if abs(delta) < 0.5 f2c = (k + delta) * fs / Nfft; else f2c = k * fs / Nfft; end else f2c = k * fs / Nfft; end fc_est = f2c / 2; % 谱线突出度：峰功率 / 两侧中值功率，工程上要求 >10 dB ref = median([S(ig(1):k-3), S(k+3:ig(end))]); sharp = 10*log10(S(k) / ref); fprintf('fc=%.5f MHz，突出度 %.1f dB\n', fc_est/1e6, sharp);

搜索区间不能用全频带最大值。平方后的噪声谱低频分量偏高，DC 附近还有残差，直接全谱取最大值容易锁到低频噪声上。flo/fhi 应由接收机扫描或信道先验给出，没有先验就分频段取局部峰再做突出度判决。

抛物线插值用峰和左右两点拟合抛物线，顶点对应的频率就是修正后的谱线位置。delta 限制在 ±0.5 个 bin 内，超出说明存在旁瓣干扰，插值不可信，退回整数 bin 结果。突出度是“是不是真载频线”的判据，单根谱线在 SNR≥0 dB 时通常超过 20 dB；低于 10 dB 时信噪比过低或区间里有强窄带干扰，应把结果标记为低置信度。

误差随观测时长和 Nfft 变化，上表量级说明插值在 5 dB 以上收益明显，0 dB 时谱线本身展宽，插值收益有限。

3.3 载频估计的三个易错点

1. 混叠没检查。fs 只比 fc 大一点时，2f_c 超过 fs/2，谱线折叠到 fs-2f_c 处，除以 2 得到的不是 fc。开工前先确认 2f_c 的名义频率小于 fs/2，拿不准就先把载频搬到低频段。
2. 把平方后的有色噪声当白噪声处理。白噪声平方后低频分量明显抬高，搜索区间尽量避开 DC，突出度阈值用中值而非均值，避免被少数高旁瓣带偏。
3. 插值公式符号写错。抛物线顶点公式分母为负是正常的，delta 可正可负，符号丢掉会得到对称的偏差。写完用频偏 0.3 bin 的测试信号自检一次，确认插值方向正确。

4. 码速率和码周期盲估计：延迟相乘谱零点与自相关副峰

4.1 延迟相乘为什么同时保留码速率与码周期信息

平方谱只有载频一根线，码速率和码周期的信息在平方过程中被当作符号抹掉了。把平方换成延迟相乘：y(t)=r(t)r(t-τ)，τ 取 1~4 个采样点，远小于一个码片。展开后 y(t) 分成三块：

y(t)= 0.5A²cos(2πf_cτ)·w(t) + 0.5A²·w(t)·cos(4πf_c t-2πf_cτ) + 噪声块

其中 w(t)=c(t)c(t-τ)，信息位 d(t) 在同一个符号内 d²(t)≡1，被自然消掉。w(t) 在码片内部为 +1，码片翻转位置附近为 -1，所以第二项相当于以 2f_c 为载频、w(t) 为基带调制的 BPSK 信号，频谱主瓣的两个一阶零点落在 2f_c±R_c，码速率从零点距离读出。w(t) 又由周期为 T_b 的扩频码决定，边界模式每 T_b 重复一次，自相关在延迟 T_b 处出现副峰，码周期从副峰位置读出。一个乘积运算同时保留两组信息，这是直扩盲估计里最常用的算子。

τ 的选取是折中：τ 越小，w(t) 越接近常数 1，2f_c 谱线越尖锐，但码片翻转的凹槽越窄，零点越浅；sps=20 时取 1~2 个采样点通常最稳。τ 接近或超过一个码片时，w(t) 变成普通随机 ±1 序列，周期结构和零点都会退化。

4.2 从延迟相乘谱的一阶零点估计码速率

tau = 2; % 延迟 2 个采样点，约 1/10 码片 y = r(1:end-tau) .* r(1+tau:end); % 延迟相乘 Nf2 = 2^nextpow2(length(y)); Y = abs(fft(y - mean(y), Nf2)).^2; % 去 DC 后取功率谱 fy = (0:Nf2-1)/Nf2*fs; Ysm = movmean(Y, 64); % 滑动平均压噪声 % 在 2fc_est 邻域定位谱线 bw = 3e6; ig = find(fy >= 2*fc_est-bw & fy <= 2*fc_est+bw); [~, p] = max(Ysm(ig)); pk = ig(1) + p - 1; % 右侧零点：峰后找第一个由降转升的局部极小 jR = pk + 3; boundR = pk + round(8e6 / (fs/Nf2)); while jR < boundR && Ysm(jR+1) <= Ysm(jR) jR = jR + 1; end fR = fy(jR); % 左侧零点：对称处理 jL = pk - 3; boundL = pk - round(8e6 / (fs/Nf2)); while jL > boundL && Ysm(jL-1) <= Ysm(jL) jL = jL - 1; end fL = fy(jL); Rc_est = ((fR - 2*fc_est) + (2*fc_est - fL)) / 2; fprintf('码速率估计 = %.4f Mchip/s（真值 %.1f）\n', Rc_est/1e6, Rc/1e6);

先减 mean(y) 再 FFT：延迟相乘输出在 DC 附近有大分量，不除掉会泄漏到 2f_c 邻域。movmean 的窗口长度直接决定零点可靠性，64 点是 sps=20、SNR=5 dB 时的经验值；SNR 到 0 dB 时加大到 128，并相应放宽搜索上界。零点搜索用“第一个由降转升的点”，不用区间内全局最小，因为第二旁瓣的极小值可能更低，全局最小会把码速率估大。左右零点距离取平均，能抵消谱线定位偏差对单侧的影响。

4.3 用自相关第一副峰定位码周期

% 码周期：延迟相乘输出的低频分量做 FFT 自相关 BW = Rc_est / 2; % 低通带宽取码速率一半 ylp = lowpass(y, BW, fs, ImpulseResponse="fir", Steepness=0.9); n = length(ylp); Yl = fft(ylp, 2^nextpow2(2*n-1)); Rw = fftshift(ifft(Yl .* conj(Yl))); Rw = Rw / max(Rw); % 零延迟处归一化 lg = (-(n-1):(n-1)) / fs; % 延迟时间轴 % 跳过主瓣（约±1.5 个码片），找第一副峰 guard = round(1.5*sps); rg = find(abs(lg) > guard & abs(lg) < 0.5*n/fs); [~, m] = max(Rw(rg)); Tb_est = abs(lg(rg(m))); Nchip_est = round(Tb_est * Rc_est); fprintf('码周期=%.3f us，折算码长=%d（真值 %d）\n', ... Tb_est*1e6, Nchip_est, N);

低通带宽取 Rc_est/2，把 2f_c 邻域彻底隔开。若直接用 y 做自相关，2f_c 邻域分量会以高频振荡形式混入，副峰位置不变但幅度被干扰。lowpass 需要 Signal Processing Toolbox（R2018a 起）；不想依赖工具箱就用 fft 把 DC 邻域谱抽出来加窗再 ifft。guard=1.5 个码片，因为自相关主瓣是宽度约 2 个码片的三角，跳过主瓣后第一个显著副峰就是码周期，m 序列的周期自相关特性保证这个峰接近满幅。折算码长 round(Tb_est*Rc_est) 是最终输出，与真实 N 一致说明该次估计成功；差 1 说明副峰定位偏了，调低通带宽或 guard 再试。

4.4 SNR 与观测长度的工程边界

边界背后是两个物理限制：码速率估计的极限分辨率由观测时长决定，频率分辨率 1/T_obs 之内零点移动可忽略，观测 3 ms 对应 0.3 kHz 量级；码周期估计要求自相关副峰的统计波动小于峰高，观测包含的码周期数越多越稳。仿真时先固定 SNR 扫观测长度，再固定长度扫 SNR，两个方向各给一组曲线，就能看到这套方法的失效面。

5. 盲估计结果验证技巧：谱线标记、周期折叠与 Monte Carlo 误差核对

5.1 把估计值标回选定谱图，用周期折叠验证参数一致性

盲估计最容易犯的错是“自洽但不正确”：算法给了三个数，画图看也像那么回事，其实零点或副峰认错了对象。先做一次廉价检查：把 2f_c、2f_c±R_c_est 和 T_b_est 三条竖线画在延迟相乘谱与自相关曲线上，确认竖线是否压住谱峰、零点和副峰，这一步能挡掉大半错误。

% 按估计码周期分段折叠，码片边界应出现规则过渡 Lper = round(Tb_est * fs); Mfold = floor(length(y) / Lper); fold = mean(reshape(y(1:Mfold*Lper), Lper, Mfold), 2); plot((0:Lper-1)/fs, abs(fold - mean(fold)), '.-');

按估计码周期分段求平均的折叠操作等效于 Mfold 次相干积累。参数正确时折叠波形在码片翻转位置出现固定凹槽；参数差一个码片时，折叠波形被随机数据抹平。这是比数值误差更直观的自验证手段。Lper 必须取整到采样点，亚采样偏差会在折叠边缘产生渐变，看到渐进式衰减说明码周期估计精度不够。

5.2 Monte Carlo 误差表与通过标准

% 20 次独立实验，仅更换噪声种子 for iter = 1:20 rng(iter); % 重新生成含噪信号 r（复用 2.2 节代码段） % 依次执行第 3、4 章的估计，记录 fc_est/Rc_est/Tb_est err_fc(iter) = abs(fc_est - fc)/fc; err_Rc(iter) = abs(Rc_est - Rc)/Rc; err_Tb(iter) = abs(Tb_est - (N/Rc)) / (N/Rc); end

典型 20 次统计结果（sps=20，观测 16 个码周期）：

通过标准按用途分层：能用于后续解扩的参数，要求码周期检出率 100%、码速率误差小于 1%；只做信号识别和粗分选时，码速率误差小于 5% 已可接受。载频误差用相对量衡量，1e-5 量级在多数应用里优于解调端需求。两个典型失败案例值得记住：码速率零点认到第二旁瓣时误差跳到 10% 以上，码周期副峰认到整数倍延迟时码长翻倍，这两类错误在谱线图上都有明显特征，画图代码保留在工程目录里常驻备用。

提示：所有估计步骤都假设目标频段内没有强窄带干扰。输入信号先过一遍陷波滤波器，把明显强于信号线谱的干扰剔掉，三个参数的估计成功率会显著回升。

本文还有配套的精品资源，点击获取