简介一套 MATLAB 代码资源聚焦 GPS、GLONASS、COMPASS 等 GNSS 扩频信号环境下的窄带干扰抑制算法模拟面向通信与导航方向的研究人员、工程师及信号处理学习者。源码包共 23 个文件体积约 35KB以 19 个 .m 脚本为主体附带 3 个备份文件与 1 个 README 说明文档代码覆盖干扰生成、扩频调制、FFT 滤波、自适应处理、误码率统计等功能模块结构清晰便于直接运行和修改验证。已有 478 人学习下载说明相关方向关注度较高。通过这套程序读者可掌握窄带干扰对导航接收机的影响机制理解基于 FFT 的频域处理与自适应滤波两类抑制思路并结合误码率、频谱对比等结果对算法性能做量化评估。整体适合课程设计、毕业设计或科研入门时作为可复现的参考实现。1. 扩频信号里的窄带干扰为什么必须单独用一个算法去抑制GPS L1 C/A 信号有 40dB 量级的解扩增益但这份增益是预留给热噪声的。把一个干信比 30dB 的单音干扰注入接收机前端相关器输出依然能抓到卫星码相位鉴别器却会抖出十几米的偏差捕获搜索矩阵的主峰还会变钝。这就是标题里那套 MATLAB 代码要处理的问题在解扩之前先把落在带内的窄带干扰分量压下去把处理增益还给真正需要的噪声。它介于射频前端滤波和捕获跟踪之间是基带信号处理里一个独立的算法层。适合做 GPS、GLONASS、COMPASS 接收机仿真或者正在为抗干扰评估发愁的工程师和学生。2. 扩频信号模式下的窄带干扰建模与仿真信号构造2.1 处理增益只对“类白噪声”有效窄带干扰会让相关峰展宽扩频系统的处理增益定义为伪码速率与信息速率的比值GPS L1 C/A 用 1.023MHz 伪码承载 50bps 电文经典公式给出的解扩增益在 43dB 左右。但这个结论成立的前提是干扰经过解扩后仍然近似均匀地散布在带宽内。窄带干扰不满足这个前提它的能量集中在载频附近的窄范围内解扩时相当于用一个带宽极窄的频谱与本地伪码频谱做卷积结果是在码相位域产生一顶带裙边的“帽子”叠在相关主峰旁边。主峰不再锐利超前滞后鉴别器输出出现偏差对应于定位结果里那几米到几十米的 gps 误差。理解这一点之后仿真信号构造的重点就清楚了不需要把星历、电文、电离层全做进去只需要把扩频码、载波、多普勒和干扰这几样东西按真实参数摆出来让干扰抑制算法有一个可信的测试床。2.2 GPS C/A、GLONASS L1 与 COMPASS B1I 的伪码参数与 MATLAB 生成三种信号在算法层面最大的差异是伪码速率、码长和多址方式。GPS L1 C/A 是码分多址每颗星用不同 Gold 码抽头GLONASS L1 用频分多址中心频率在 1602MHz 基础上按 0.5625MHz 步进分配COMPASS 也就是北斗二号时期的 B1I 信号采用 2.046MHz 伪码速率的码分多址。信号伪码速率伪码长度载波/多址相对处理增益GPS L1 C/A1.023 MHz10231575.42 MHz / CDMA43 dBGLONASS L10.511 MHz51116020.5625k MHz / FDMA40 dBCOMPASS B1I2.046 MHz20461561.098 MHz / CDMA46 dB仿真时我一般先用 GPS 的 C/A 码把算法链路跑通再改码速率和中心频率验证通用性。GPS 的 C/A 码生成是一个 10 级移位寄存器结构G1 多项式为 1x^3x^10G2 多项式为 1x^2x^3x^6x^8x^9x^10不同 PRN 靠 G2 的不同抽头组合区分。function ca generate_ca_code(prn) % 生成GPS L1 C/A码返回1023个±1码片 g1 ones(1, 10); g2 ones(1, 10); % PRN1使用G2的第2和第6个抽头其它PRN查表替换这两个数 tap1 2; tap2 6; ca zeros(1, 1023); for i 1:1023 ca(i) xor(g1(10), xor(g2(tap1), g2(tap2))); g1 [xor(g1(3), g1(10)), g1(1:9)]; % G1反馈 g2 [mod(sum(g2([2 3 6 8 9 10])), 2), g2(1:9)]; % G2反馈 end ca 2*ca - 1; end这段代码里最关键的是 G2 反馈抽头的异或顺序写错一个索引整个码序列就变成另一族不相关的序列捕获阶段会搜不出主峰。生成后建议先用[~, lag] xcorr(ca, ca)看一眼自相关副峰确认主峰单值、副峰小于主峰 21dB 再往下走。2.3 接收信号构造SNR、JSR 与干扰类型的换算窄带干扰抑制算法的验证需要同时给出信噪比和干信比。信号功率按 ±1 码片归一化为 1噪声标准差用 SNR 换算干扰幅值用 JSR 换算三者叠加后就是接收机中频采样信号。fs 20e6; % 采样率20MHz N 20480; % 1ms采样点数 t (0:N-1)/fs; ca_code generate_ca_code(1); ca20 repelem(ca_code, 20); % 每码片20点 ca20 ca20(1:N); % 截断到整数个码片周期 SNR -20; % 信噪比-20dB JSR 40; % 干信比40dB sigma 10^(-SNR/20); % 噪声幅度 A_j 10^(JSR/20); % 干扰幅度 c1 A_j * cos(2*pi*4.2e6*t); % 单音 c2 A_j * cos(2*pi*4.8e6*t 0.5*sin(2*pi*25e3*t)); % 窄带FM rx ca20 sigma*randn(1, N) c1 c2;JSR 40dB 对应干扰幅值是信号的 100 倍在时域图上扩频信号完全被淹没只能靠频谱观察。c2 的调制指数 0.5、调制频率 25kHz产生的干扰瞬时频率在 4.8MHz 附近摆动3dB 带宽大致是几十 kHz 量级比信号带宽窄两个数量级符合窄带假设。用spectrogram(rx, 256, 128, 256, fs)能直接看到一条亮线加一条抖动亮带这就是后续算法要处理的两种典型窄带形态。3. 频域陷波与时域预测两类窄带干扰抑制算法的 MATLAB 实现3.1 频域陷波加窗、FFT、中位数门限、重叠保留频域处理的思路是直接把频谱上的凸起按下去。FFT 的每个频点相当于一个窄带滤波器单音干扰只占少数几个频率单元把超过门限的频点置零再做 IFFT 恢复时域信号窄带干扰就被剔除。工程实现不能一帧一帧独立做否则帧边界处会出现瞬态跳变我常用的做法是重叠保留加窗再重叠相加。function y notch_freq(x, frame_len, overlap, thr) % 频域窄带干扰抑制 % x: 输入信号frame_len: FFT长度overlap: 重叠比thr: 门限系数 n length(x); hop max(round(frame_len * (1 - overlap)), 1); w hann(frame_len, periodic); y zeros(n, 1); win_sum zeros(n, 1); idx 1; while idx frame_len - 1 n seg x(idx:idxframe_len-1) .* w; spec fft(seg); amp abs(spec); med median(amp); % 中位数估计噪声底 mask amp thr * med; % 超过门限判为干扰 spec(mask) 0; % 干扰频点置零 out real(ifft(spec)) .* w; y(idx:idxframe_len-1) y(idx:idxframe_len-1) out; win_sum(idx:idxframe_len-1) win_sum(idx:idxframe_len-1) w.^2; idx idx hop; end y y ./ (win_sum eps); % 归一化窗增益 end门限用中位数而不是均值是因为强干扰抬高了频谱均值后门限会跟着漂中位数对少数强谱线不敏感更能代表噪声底。thr 默认取 6意思是只陷那些功率是噪声底 6 倍以上的频点。frame_len 取 1024 到 4096频谱分辨率分别对应约 20kHz 到 5kHz单音用 1024 就够两个频率靠得很近的干扰需要更长的帧来分辨。overlap 取 0.5 能抑制大部分窗泄漏引起的边界振铃如果后续还要做码相位测量建议提高到 0.75。3.2 时域 LMS 自适应预测用窄带信号的可预测性做对消窄带信号在时间上有强相关性几个码片之后它依然能由过去的采样点预测出来扩频信号经过伪码调制后接近白噪声时移一个码片以上就不再有相关性。利用这一点可以用一个延迟后的参考信号去预测当前输入中的窄带分量再用输入减去预测值剩下的就是被保留的扩频信号和噪声。function [y, w] lms_notch(x, D, order, mu) % 归一化LMS窄带干扰预测对消 % D: 预测延迟(采样点)order: 滤波器阶数mu: 归一化步长 n length(x); w zeros(order, 1); y zeros(n, 1); for k order D 1 : n ref x(k-D:-1:k-D-order1); % 延迟D后的参考向量 y_hat ref. * w; err x(k) - y_hat; % 预测误差作为输出 w w mu * err * ref / (ref.*ref 1e-8); % NLMS更新 y(k) err; end endD 的取值决定了“可预测”和“不相关”之间的分界。D 取 1 时预测效果最强单音干扰对消深度最大但也会牺牲一部分扩频信号的低频分量D 取一个伪码码片对应的时间也就是本例里的 20 个采样点能在抑制窄带干扰的同时把对扩频信号的损伤降到最小。阶数 order 取 16 到 64覆盖窄带干扰的时域拖尾长度。步长 mu 在 NLMS 里取 0.05 到 0.5值越大收敛越快但权值噪声也越大干扰抑制后输出端的底噪会被抬高。3.3 两种算法在 GLONASS 频分多址下的行为差异GLONASS 每颗卫星用不同载波频率接收机前端通常按频点分别下变频。频域陷波对每个频点单独做时不需要互相协调门限处理天然隔离实现最简单。LMS 预测器则要注意延迟 D 与码速率的关系GLONASS 伪码速率只有 0.511MHz同样一个码片对应约 40 个采样点D 要相应增大。两类算法的取舍可以按下面这张表来判断。对比项频域陷波时域LMS预测单音干扰抑制深门限好调抑制深收敛后残余小窄带FM干扰需要门限覆盖整个频带裙边自适应跟随效果更好多音干扰同时陷多个峰简单可靠需要更高阶数权值互相牵扯计算复杂度FFT 帧处理固定开销每个采样点都要做滤波更新对扩频信号损伤陷波带宽过宽会削掉信号频谱D 和阶数选择不当会压低信号我一般在干扰个数已知、频谱形态固定的场景优先用频域陷波在干扰类型动态变化的场景改用 LMS。频域陷波最怕的是把信号频谱当成干扰削掉这一点在下一章讲参数时单独展开。4. 门限、步长与子空间法把抑制算法调到不伤信号4.1 三个必调参数的边界与后果频域陷波的门限系数 thr 是最容易调崩的参数。thr 太小会把正常噪声的峰值当干扰置零频谱上形成一个个坑解扩后相关峰出现额外副瓣太大则只压掉了干扰最强的中心频点窄带 FM 的裙边残留下来抑制增益达不到预期。经验法则是先用spectrogram看干扰的 3dB 带宽占多少个频点再让 thr 保证门限只覆盖这些频点而不是全部频谱的若干分之一。重叠比和 FFI 长度是第二组参数。frame_len 决定频谱分辨率两个中心频率相差不到一个分辨率带宽的干扰在频谱上会连成一片这时候加长帧比调门限更有效。overlap 不影响频谱分辨率只影响帧间连续性码环对相位跳变敏感时把 overlap 从 0.5 提到 0.75 能明显减少跟踪环路的抖动。参数典型范围调大的后果调小的后果frame_len1024 ~ 4096分辨率高帧延迟增加无法分离邻近干扰overlap0.25 ~ 0.75帧边界更连续算力增加出现振铃伪迹thr3 ~ 10残留弱干扰信号损伤小压得干净但可能削信号mu0.05 ~ 0.5收敛快权值噪声大收敛慢动态干扰跟不住D1 个码片附近对扩频信号损伤小单音抑制更深4.2 用捕获相关峰和载噪比判断抑制效果判断抑制算法有没有把信号一起削掉只看输出频谱是不够的。标准做法是把抑制后的信号直接接入捕获相关器用本地伪码做循环相关观察主峰形态。% 假设rx_notch是抑制后的20MHz采样信号 ca_ref ca20; % 本地伪码 corr_full ifft(fft(rx_notch) .* conj(fft(ca_ref))); corr_abs abs(corr_full) / N; [peak, peak_idx] max(corr_abs); % 取主峰两侧一个码片的位置做噪底估计 side_a corr_abs(mod(peak_idx-20, N) 1); side_b corr_abs(mod(peak_idx20, N) 1); floor_est (side_a side_b) / 2; peak_to_floor 10*log10(peak / floor_est); % 主峰与噪底比这个主峰与噪底的比值可以直接反映抑制后信号还能不能被捕获环路锁定。原始信号在无干扰下通常能到 25dB 以上加入窄带干扰被算法抑制后如果还大于 20dB说明信号频谱基本保住了。低于 15dB 时要么门限过小削了信号要么 LMS 阶数不够残余干扰仍然过大需要回到参数表逐项排查。载噪比估计更严格的做法是用宽带功率和窄带功率的比值来计算但仿真阶段用上述峰值噪底比做相对比较已经足够定位问题。4.3 干扰数多时改用子空间投影当带内同时存在三四个以上窄带干扰频域陷波需要逐个定位门限LMS 阶数需要抬得很高权值间互相牵扯导致收敛变慢。这时候剩下的方案是子空间法把观测数据按快拍排列求协方差矩阵特征值大的一组对应干扰子空间特征值小的对应扩频信号与噪声把数据投影到干扰子空间的正交补上即可完成抑制。M 32; % 子空间维数 X buffer(rx_notch, M, M-1); % 按帧排列快拍矩阵 X X(:, 1:floor(length(rx_notch)/M)); % 截齐 R X * X / size(X, 2); % 协方差矩阵 [V, d] eig(R); d diag(d); [~, idx] sort(d, descend); V V(:, idx); p 2; % 干扰个数估计值 I_sub V(:, 1:p); % 干扰子空间 P_proj eye(M) - I_sub * I_sub; % 正交投影矩阵 y_sub P_proj * rx_notch(1:M:end); % 投影后输出这里的 p 取干扰个数的估计值实践中可以用特征值谱观察大特征值的个数或者直接用 MDL 准则计算。p 设大会把信号子空间也投影掉一部分输出信噪比下降p 设小则干扰抑制不净。子空间法的优点是稳态性能好缺点是协方差矩阵估计需要足够长的数据动态干扰下性能会打折。5. 从仿真到工程验证级联结构与时频结合的调试技巧5.1 频域粗陷波加 LMS 细对消的级联顺序单一算法很难同时对付强单音和宽带稍大的窄带 FM。频域陷波对强单音的深度抑制很干净但门限固定时对动态窄带 FM 的裙边压不净LMS 能跟踪动态干扰但强干扰下权值噪声会把输出底噪抬起来。常见做法是级联两级第一级频域陷波把干扰总功率压到接近噪声底第二级 LMS 只负责把残余的动态窄带分量对消掉。顺序不要颠倒反过来让 LMS 先面对百倍强度的干扰收敛后权值已经扭曲再进频域陷波会额外削掉信号附近的有效频谱。级联后两个算法各让一步参数频域门限从 6 放宽到 8只压最强的谱峰LMS 阶数降到 16步长调小到 0.1防止它对已被压平的频谱过度反应。5.2 用二次陷波确认抑制是否收敛判断抑制效果有一个很直接的验证技巧把抑制后的输出再送入同一个陷波函数跑一遍对比两次输出的功率谱。第二遍频谱与第一遍几乎一致说明第一遍已经把窄带分量压到了门限以下第二遍又能看到明显凹陷说明门限设定偏高有干扰成分残留。这个办法不需要额外的参考信号在仿真阶段比读数更快。逐帧观察时用spectrogram看两个地方一是干扰频点附近是否还有残余亮线二是陷波带的宽度是否明显大于干扰本身的 3dB 带宽。正常范围内陷波带比干扰带宽宽 1 到 2 个频率分辨率单元宽太多就说明把信号频谱一并削掉了要回退 thr 或者改用带保护间隔的邻域置零也就是只陷峰值附近 3 个频点而不是整个连续超门限区域。把抑制后的信号再接入原有捕获跟踪环路对比无干扰情况下相关主峰位置偏移偏移控制在一个采样点以内就可以把抑制算法正式放出链路。本文还有配套的精品资源点击获取 基于MATLAB的GPS窄带干扰抑制算法实现
简介一套 MATLAB 代码资源聚焦 GPS、GLONASS、COMPASS 等 GNSS 扩频信号环境下的窄带干扰抑制算法模拟面向通信与导航方向的研究人员、工程师及信号处理学习者。源码包共 23 个文件体积约 35KB以 19 个 .m 脚本为主体附带 3 个备份文件与 1 个 README 说明文档代码覆盖干扰生成、扩频调制、FFT 滤波、自适应处理、误码率统计等功能模块结构清晰便于直接运行和修改验证。已有 478 人学习下载说明相关方向关注度较高。通过这套程序读者可掌握窄带干扰对导航接收机的影响机制理解基于 FFT 的频域处理与自适应滤波两类抑制思路并结合误码率、频谱对比等结果对算法性能做量化评估。整体适合课程设计、毕业设计或科研入门时作为可复现的参考实现。1. 扩频信号里的窄带干扰为什么必须单独用一个算法去抑制GPS L1 C/A 信号有 40dB 量级的解扩增益但这份增益是预留给热噪声的。把一个干信比 30dB 的单音干扰注入接收机前端相关器输出依然能抓到卫星码相位鉴别器却会抖出十几米的偏差捕获搜索矩阵的主峰还会变钝。这就是标题里那套 MATLAB 代码要处理的问题在解扩之前先把落在带内的窄带干扰分量压下去把处理增益还给真正需要的噪声。它介于射频前端滤波和捕获跟踪之间是基带信号处理里一个独立的算法层。适合做 GPS、GLONASS、COMPASS 接收机仿真或者正在为抗干扰评估发愁的工程师和学生。2. 扩频信号模式下的窄带干扰建模与仿真信号构造2.1 处理增益只对“类白噪声”有效窄带干扰会让相关峰展宽扩频系统的处理增益定义为伪码速率与信息速率的比值GPS L1 C/A 用 1.023MHz 伪码承载 50bps 电文经典公式给出的解扩增益在 43dB 左右。但这个结论成立的前提是干扰经过解扩后仍然近似均匀地散布在带宽内。窄带干扰不满足这个前提它的能量集中在载频附近的窄范围内解扩时相当于用一个带宽极窄的频谱与本地伪码频谱做卷积结果是在码相位域产生一顶带裙边的“帽子”叠在相关主峰旁边。主峰不再锐利超前滞后鉴别器输出出现偏差对应于定位结果里那几米到几十米的 gps 误差。理解这一点之后仿真信号构造的重点就清楚了不需要把星历、电文、电离层全做进去只需要把扩频码、载波、多普勒和干扰这几样东西按真实参数摆出来让干扰抑制算法有一个可信的测试床。2.2 GPS C/A、GLONASS L1 与 COMPASS B1I 的伪码参数与 MATLAB 生成三种信号在算法层面最大的差异是伪码速率、码长和多址方式。GPS L1 C/A 是码分多址每颗星用不同 Gold 码抽头GLONASS L1 用频分多址中心频率在 1602MHz 基础上按 0.5625MHz 步进分配COMPASS 也就是北斗二号时期的 B1I 信号采用 2.046MHz 伪码速率的码分多址。信号伪码速率伪码长度载波/多址相对处理增益GPS L1 C/A1.023 MHz10231575.42 MHz / CDMA43 dBGLONASS L10.511 MHz51116020.5625k MHz / FDMA40 dBCOMPASS B1I2.046 MHz20461561.098 MHz / CDMA46 dB仿真时我一般先用 GPS 的 C/A 码把算法链路跑通再改码速率和中心频率验证通用性。GPS 的 C/A 码生成是一个 10 级移位寄存器结构G1 多项式为 1x^3x^10G2 多项式为 1x^2x^3x^6x^8x^9x^10不同 PRN 靠 G2 的不同抽头组合区分。function ca generate_ca_code(prn) % 生成GPS L1 C/A码返回1023个±1码片 g1 ones(1, 10); g2 ones(1, 10); % PRN1使用G2的第2和第6个抽头其它PRN查表替换这两个数 tap1 2; tap2 6; ca zeros(1, 1023); for i 1:1023 ca(i) xor(g1(10), xor(g2(tap1), g2(tap2))); g1 [xor(g1(3), g1(10)), g1(1:9)]; % G1反馈 g2 [mod(sum(g2([2 3 6 8 9 10])), 2), g2(1:9)]; % G2反馈 end ca 2*ca - 1; end这段代码里最关键的是 G2 反馈抽头的异或顺序写错一个索引整个码序列就变成另一族不相关的序列捕获阶段会搜不出主峰。生成后建议先用[~, lag] xcorr(ca, ca)看一眼自相关副峰确认主峰单值、副峰小于主峰 21dB 再往下走。2.3 接收信号构造SNR、JSR 与干扰类型的换算窄带干扰抑制算法的验证需要同时给出信噪比和干信比。信号功率按 ±1 码片归一化为 1噪声标准差用 SNR 换算干扰幅值用 JSR 换算三者叠加后就是接收机中频采样信号。fs 20e6; % 采样率20MHz N 20480; % 1ms采样点数 t (0:N-1)/fs; ca_code generate_ca_code(1); ca20 repelem(ca_code, 20); % 每码片20点 ca20 ca20(1:N); % 截断到整数个码片周期 SNR -20; % 信噪比-20dB JSR 40; % 干信比40dB sigma 10^(-SNR/20); % 噪声幅度 A_j 10^(JSR/20); % 干扰幅度 c1 A_j * cos(2*pi*4.2e6*t); % 单音 c2 A_j * cos(2*pi*4.8e6*t 0.5*sin(2*pi*25e3*t)); % 窄带FM rx ca20 sigma*randn(1, N) c1 c2;JSR 40dB 对应干扰幅值是信号的 100 倍在时域图上扩频信号完全被淹没只能靠频谱观察。c2 的调制指数 0.5、调制频率 25kHz产生的干扰瞬时频率在 4.8MHz 附近摆动3dB 带宽大致是几十 kHz 量级比信号带宽窄两个数量级符合窄带假设。用spectrogram(rx, 256, 128, 256, fs)能直接看到一条亮线加一条抖动亮带这就是后续算法要处理的两种典型窄带形态。3. 频域陷波与时域预测两类窄带干扰抑制算法的 MATLAB 实现3.1 频域陷波加窗、FFT、中位数门限、重叠保留频域处理的思路是直接把频谱上的凸起按下去。FFT 的每个频点相当于一个窄带滤波器单音干扰只占少数几个频率单元把超过门限的频点置零再做 IFFT 恢复时域信号窄带干扰就被剔除。工程实现不能一帧一帧独立做否则帧边界处会出现瞬态跳变我常用的做法是重叠保留加窗再重叠相加。function y notch_freq(x, frame_len, overlap, thr) % 频域窄带干扰抑制 % x: 输入信号frame_len: FFT长度overlap: 重叠比thr: 门限系数 n length(x); hop max(round(frame_len * (1 - overlap)), 1); w hann(frame_len, periodic); y zeros(n, 1); win_sum zeros(n, 1); idx 1; while idx frame_len - 1 n seg x(idx:idxframe_len-1) .* w; spec fft(seg); amp abs(spec); med median(amp); % 中位数估计噪声底 mask amp thr * med; % 超过门限判为干扰 spec(mask) 0; % 干扰频点置零 out real(ifft(spec)) .* w; y(idx:idxframe_len-1) y(idx:idxframe_len-1) out; win_sum(idx:idxframe_len-1) win_sum(idx:idxframe_len-1) w.^2; idx idx hop; end y y ./ (win_sum eps); % 归一化窗增益 end门限用中位数而不是均值是因为强干扰抬高了频谱均值后门限会跟着漂中位数对少数强谱线不敏感更能代表噪声底。thr 默认取 6意思是只陷那些功率是噪声底 6 倍以上的频点。frame_len 取 1024 到 4096频谱分辨率分别对应约 20kHz 到 5kHz单音用 1024 就够两个频率靠得很近的干扰需要更长的帧来分辨。overlap 取 0.5 能抑制大部分窗泄漏引起的边界振铃如果后续还要做码相位测量建议提高到 0.75。3.2 时域 LMS 自适应预测用窄带信号的可预测性做对消窄带信号在时间上有强相关性几个码片之后它依然能由过去的采样点预测出来扩频信号经过伪码调制后接近白噪声时移一个码片以上就不再有相关性。利用这一点可以用一个延迟后的参考信号去预测当前输入中的窄带分量再用输入减去预测值剩下的就是被保留的扩频信号和噪声。function [y, w] lms_notch(x, D, order, mu) % 归一化LMS窄带干扰预测对消 % D: 预测延迟(采样点)order: 滤波器阶数mu: 归一化步长 n length(x); w zeros(order, 1); y zeros(n, 1); for k order D 1 : n ref x(k-D:-1:k-D-order1); % 延迟D后的参考向量 y_hat ref. * w; err x(k) - y_hat; % 预测误差作为输出 w w mu * err * ref / (ref.*ref 1e-8); % NLMS更新 y(k) err; end endD 的取值决定了“可预测”和“不相关”之间的分界。D 取 1 时预测效果最强单音干扰对消深度最大但也会牺牲一部分扩频信号的低频分量D 取一个伪码码片对应的时间也就是本例里的 20 个采样点能在抑制窄带干扰的同时把对扩频信号的损伤降到最小。阶数 order 取 16 到 64覆盖窄带干扰的时域拖尾长度。步长 mu 在 NLMS 里取 0.05 到 0.5值越大收敛越快但权值噪声也越大干扰抑制后输出端的底噪会被抬高。3.3 两种算法在 GLONASS 频分多址下的行为差异GLONASS 每颗卫星用不同载波频率接收机前端通常按频点分别下变频。频域陷波对每个频点单独做时不需要互相协调门限处理天然隔离实现最简单。LMS 预测器则要注意延迟 D 与码速率的关系GLONASS 伪码速率只有 0.511MHz同样一个码片对应约 40 个采样点D 要相应增大。两类算法的取舍可以按下面这张表来判断。对比项频域陷波时域LMS预测单音干扰抑制深门限好调抑制深收敛后残余小窄带FM干扰需要门限覆盖整个频带裙边自适应跟随效果更好多音干扰同时陷多个峰简单可靠需要更高阶数权值互相牵扯计算复杂度FFT 帧处理固定开销每个采样点都要做滤波更新对扩频信号损伤陷波带宽过宽会削掉信号频谱D 和阶数选择不当会压低信号我一般在干扰个数已知、频谱形态固定的场景优先用频域陷波在干扰类型动态变化的场景改用 LMS。频域陷波最怕的是把信号频谱当成干扰削掉这一点在下一章讲参数时单独展开。4. 门限、步长与子空间法把抑制算法调到不伤信号4.1 三个必调参数的边界与后果频域陷波的门限系数 thr 是最容易调崩的参数。thr 太小会把正常噪声的峰值当干扰置零频谱上形成一个个坑解扩后相关峰出现额外副瓣太大则只压掉了干扰最强的中心频点窄带 FM 的裙边残留下来抑制增益达不到预期。经验法则是先用spectrogram看干扰的 3dB 带宽占多少个频点再让 thr 保证门限只覆盖这些频点而不是全部频谱的若干分之一。重叠比和 FFI 长度是第二组参数。frame_len 决定频谱分辨率两个中心频率相差不到一个分辨率带宽的干扰在频谱上会连成一片这时候加长帧比调门限更有效。overlap 不影响频谱分辨率只影响帧间连续性码环对相位跳变敏感时把 overlap 从 0.5 提到 0.75 能明显减少跟踪环路的抖动。参数典型范围调大的后果调小的后果frame_len1024 ~ 4096分辨率高帧延迟增加无法分离邻近干扰overlap0.25 ~ 0.75帧边界更连续算力增加出现振铃伪迹thr3 ~ 10残留弱干扰信号损伤小压得干净但可能削信号mu0.05 ~ 0.5收敛快权值噪声大收敛慢动态干扰跟不住D1 个码片附近对扩频信号损伤小单音抑制更深4.2 用捕获相关峰和载噪比判断抑制效果判断抑制算法有没有把信号一起削掉只看输出频谱是不够的。标准做法是把抑制后的信号直接接入捕获相关器用本地伪码做循环相关观察主峰形态。% 假设rx_notch是抑制后的20MHz采样信号 ca_ref ca20; % 本地伪码 corr_full ifft(fft(rx_notch) .* conj(fft(ca_ref))); corr_abs abs(corr_full) / N; [peak, peak_idx] max(corr_abs); % 取主峰两侧一个码片的位置做噪底估计 side_a corr_abs(mod(peak_idx-20, N) 1); side_b corr_abs(mod(peak_idx20, N) 1); floor_est (side_a side_b) / 2; peak_to_floor 10*log10(peak / floor_est); % 主峰与噪底比这个主峰与噪底的比值可以直接反映抑制后信号还能不能被捕获环路锁定。原始信号在无干扰下通常能到 25dB 以上加入窄带干扰被算法抑制后如果还大于 20dB说明信号频谱基本保住了。低于 15dB 时要么门限过小削了信号要么 LMS 阶数不够残余干扰仍然过大需要回到参数表逐项排查。载噪比估计更严格的做法是用宽带功率和窄带功率的比值来计算但仿真阶段用上述峰值噪底比做相对比较已经足够定位问题。4.3 干扰数多时改用子空间投影当带内同时存在三四个以上窄带干扰频域陷波需要逐个定位门限LMS 阶数需要抬得很高权值间互相牵扯导致收敛变慢。这时候剩下的方案是子空间法把观测数据按快拍排列求协方差矩阵特征值大的一组对应干扰子空间特征值小的对应扩频信号与噪声把数据投影到干扰子空间的正交补上即可完成抑制。M 32; % 子空间维数 X buffer(rx_notch, M, M-1); % 按帧排列快拍矩阵 X X(:, 1:floor(length(rx_notch)/M)); % 截齐 R X * X / size(X, 2); % 协方差矩阵 [V, d] eig(R); d diag(d); [~, idx] sort(d, descend); V V(:, idx); p 2; % 干扰个数估计值 I_sub V(:, 1:p); % 干扰子空间 P_proj eye(M) - I_sub * I_sub; % 正交投影矩阵 y_sub P_proj * rx_notch(1:M:end); % 投影后输出这里的 p 取干扰个数的估计值实践中可以用特征值谱观察大特征值的个数或者直接用 MDL 准则计算。p 设大会把信号子空间也投影掉一部分输出信噪比下降p 设小则干扰抑制不净。子空间法的优点是稳态性能好缺点是协方差矩阵估计需要足够长的数据动态干扰下性能会打折。5. 从仿真到工程验证级联结构与时频结合的调试技巧5.1 频域粗陷波加 LMS 细对消的级联顺序单一算法很难同时对付强单音和宽带稍大的窄带 FM。频域陷波对强单音的深度抑制很干净但门限固定时对动态窄带 FM 的裙边压不净LMS 能跟踪动态干扰但强干扰下权值噪声会把输出底噪抬起来。常见做法是级联两级第一级频域陷波把干扰总功率压到接近噪声底第二级 LMS 只负责把残余的动态窄带分量对消掉。顺序不要颠倒反过来让 LMS 先面对百倍强度的干扰收敛后权值已经扭曲再进频域陷波会额外削掉信号附近的有效频谱。级联后两个算法各让一步参数频域门限从 6 放宽到 8只压最强的谱峰LMS 阶数降到 16步长调小到 0.1防止它对已被压平的频谱过度反应。5.2 用二次陷波确认抑制是否收敛判断抑制效果有一个很直接的验证技巧把抑制后的输出再送入同一个陷波函数跑一遍对比两次输出的功率谱。第二遍频谱与第一遍几乎一致说明第一遍已经把窄带分量压到了门限以下第二遍又能看到明显凹陷说明门限设定偏高有干扰成分残留。这个办法不需要额外的参考信号在仿真阶段比读数更快。逐帧观察时用spectrogram看两个地方一是干扰频点附近是否还有残余亮线二是陷波带的宽度是否明显大于干扰本身的 3dB 带宽。正常范围内陷波带比干扰带宽宽 1 到 2 个频率分辨率单元宽太多就说明把信号频谱一并削掉了要回退 thr 或者改用带保护间隔的邻域置零也就是只陷峰值附近 3 个频点而不是整个连续超门限区域。把抑制后的信号再接入原有捕获跟踪环路对比无干扰情况下相关主峰位置偏移偏移控制在一个采样点以内就可以把抑制算法正式放出链路。本文还有配套的精品资源点击获取