DSSS抗窄带干扰MATLAB仿真：从扩频增益到频域置零与半解析BER

简介这是一份面向无线通信学习者与MATLAB仿真实践者的DSSS扩频通信抗窄带干扰仿真代码包聚焦直接序列扩频系统在窄带干扰环境下的建模与性能分析适合通信工程、电子信息类专业学生及需要理解扩频抗干扰机制的开发者参考。包内共3个文件以2个m脚本和1个md说明文档为主脚本承担扩频调制、干扰引入与误码率计算等核心仿真流程说明文档用于交代项目结构与运行方式压缩包整体约3KB体量轻便、便于快速上手。资源围绕伪随机码生成、信息信号扩频调制、窄带干扰建模、解扩与匹配滤波、误码率与信噪比评估等环节展开读者可据此复现完整的DSSS收发链路观察窄带干扰对系统性能的影响并在此基础上尝试自适应均衡、干扰抵消等抗干扰策略的改进。目前已有467人学习下载适合作为通信原理课程实验或扩频通信入门仿真的参考素材。1. DSSS 抗窄带干扰从一条被压扁的频谱说起做扩频通信仿真的人大概率都遇到过这样的场景接收端解扩之后星座图糊成一团误码率曲线在某个信噪比之后突然翘起来怎么调滤波器都没用。把频谱仪打开一看某个窄带信号像一根针一样扎在扩频带宽里功率还不低。这就是窄带干扰NBI对 DSSS 系统最典型的杀伤方式——它不需要覆盖整个带宽只要在扩频码的一个频点上持续注入能量解扩时的相关运算就会把这根针的能量搬回基带直接抬高判决前的噪声底。DSSS_matlab 这类工程的核心就是把扩频—加窄带干扰—抗干扰—解扩—误码统计这条链路在 MATLAB 里跑通让你能定量看到不同干扰策略下系统到底还剩多少余量。它适合三类人正在做扩频通信课程设计或毕设的学生、需要快速验证抗干扰算法如自适应陷波、频域置零的算法工程师、以及想把 MATLAB 仿真结果和硬件实测对齐的射频调试人员。整条链路不复杂但参数之间的耦合关系很容易让人翻车下面按先立原理、再搭链路、最后抠参数的顺序拆开讲。2. 扩频增益与窄带干扰的对抗关系先算清楚再动手2.1 处理增益到底能扛多少干扰DSSS 抗窄带干扰的底气来自处理增益。设扩频码速率为 ( R_c )信息速率为 ( R_b )则处理增益 ( G_p 10\lg(R_c/R_b) ) dB。以常见的 ( R_c 10.23 ) Mcps、( R_b 50 ) kbps 为例( G_p \approx 23 ) dB。这意味着一个功率比信号高 20 dB 的窄带干扰解扩后会被压低到信号之下约 3 dB。但这里有个前提干扰必须是窄带的也就是它的带宽远小于扩频带宽。如果干扰带宽接近甚至超过 ( R_c )处理增益就失效了因为相关器无法把干扰能量均匀摊开。所以仿真里第一步不是写代码而是先确定干扰带宽 ( B_j ) 和扩频带宽 ( B_s ) 的比值。我一般会把这个比值控制在 0.01 到 0.1 之间这样处理增益的模型才成立。另一个容易被忽略的点是干扰频率位置。窄带干扰落在扩频信号频谱的中心和落在边缘对解扩的影响不一样。中心频点的干扰经过相关器后产生的直流分量最大对判决的破坏最直接。仿真时至少要扫三个频点中心、( \pm B_s/4 )、( \pm B_s/2 )才能看出系统的最坏情况。2.2 在 MATLAB 里搭一条最小可跑链路下面这段代码是整条链路的最小骨架包含扩频、加窄带干扰、解扩和误码统计。先跑通它再往里加抗干扰模块。% DSSS 最小链路扩频 窄带干扰 解扩 误码统计 clear; clc; rng(42); % 固定随机种子保证结果可复现 % ---- 参数区 ---- Rb 50e3; % 信息速率 50 kbps Rc 10.23e6; % 扩频码速率 10.23 Mcps sf Rc / Rb; % 扩频因子约 204.6 Nbit 2000; % 仿真比特数 fc 2e6; % 载波频率仅用于干扰频点定位 fj 2e6; % 窄带干扰中心频率先对准载波 Bj 200e3; % 窄带干扰带宽 200 kHz JSR 20; % 干信比 dB % ---- 生成信息比特与扩频码 ---- bits randi([0 1], Nbit, 1); pn randi([0 1], sf, 1) * 2 - 1; % 双极性 PN 码 pn repmat(pn, ceil(Nbit*sf/length(pn)), 1); pn pn(1:Nbit*sf); % ---- 扩频 ---- tx zeros(Nbit*sf, 1); for k 1:Nbit tx((k-1)*sf1 : k*sf) (2*bits(k)-1) * pn((k-1)*sf1 : k*sf); end % ---- 生成窄带干扰带通滤波后的高斯噪声 ---- fs Rc; % 采样率取码速率 n (0:length(tx)-1).; noise randn(size(tx)); % 用简单 IIR 带通近似窄带干扰中心 fj带宽 Bj [b, a] butter(4, [(fj-Bj/2)/(fs/2), (fjBj/2)/(fs/2)], bandpass); nbi filter(b, a, noise); nbi nbi / std(nbi); % 归一化 sig_pow mean(tx.^2); nbi nbi * sqrt(sig_pow * 10^(JSR/10)); % ---- 叠加干扰与噪声 ---- SNR 10; % 信噪比 dB awgn randn(size(tx)) * sqrt(sig_pow / 10^(SNR/10)); rx tx nbi awgn; % ---- 解扩 ---- rx_des rx .* pn; % 按扩频因子积分判决 rx_int reshape(rx_des, sf, Nbit); decision sum(rx_int, 1).; bits_hat decision 0; % ---- 误码统计 ---- ber sum(bits_hat ~ bits) / Nbit; fprintf(JSR %d dB, SNR %d dB, BER %.4e\n, JSR, SNR, ber);这段代码的逻辑很直白先按扩频因子把每个比特铺成 ( sf ) 个码片再叠加一个带通滤波后的高斯噪声作为窄带干扰最后用同一组 PN 码做相关解扩。参数区里JSR是干信比SNR是信噪比Bj是干扰带宽。跑一遍你会看到JSR 到 20 dB 时 BER 已经明显恶化但还没到完全不可用的程度这正是处理增益在起作用。需要提醒的是butter的阶数不要设太高。4 阶已经能给出比较干净的窄带干扰阶数再高会引入明显的相位失真反而让干扰模型偏离实际。另外fs取Rc是为了让每个码片一个采样点如果你后续要做频域抗干扰采样率至少要取 ( 2R_c ) 以上否则频谱会混叠。2.3 干扰带宽和干信比怎么扫才有意义单跑一个点看不出趋势实际分析时我会固定 SNR然后扫 JSR 从 0 到 30 dB每个点跑 10 次取平均 BER。下面这段是扫描脚本的核心部分JSR_list 0:2:30; ber_avg zeros(size(JSR_list)); for idx 1:length(JSR_list) ber_tmp zeros(10,1); for trial 1:10 % 这里调用上面的链路函数传入 JSR_list(idx) ber_tmp(trial) run_dsss_link(JSR_list(idx), 10, 200e3); end ber_avg(idx) mean(ber_tmp); end semilogy(JSR_list, ber_avg, -o); xlabel(JSR (dB)); ylabel(BER); grid on;扫描时要注意两点一是每个 JSR 点至少跑 10 次独立噪声实现否则 BER 曲线会抖得没法看二是当 BER 低于 ( 10^{-4} ) 时2000 个比特已经不够统计了需要把Nbit加到 ( 10^5 ) 以上或者改用半解析方法。我一般会在 BER 降到 ( 10^{-3} ) 以下时切换成更大的比特数避免曲线尾部全是零。3. 频域抗窄带干扰门限置零与自适应陷波的 MATLAB 实现3.1 为什么频域处理比时域更适合窄带干扰窄带干扰在时域上表现为一个缓慢变化的周期性分量直接做时域自适应滤波需要很高的阶数才能跟踪。但换到频域它就是一个或几个尖锐的谱峰用门限检测加置零就能干掉。这就是频域抗干扰的基本思路对接收信号做 FFT找到超过门限的谱线把它们置零或衰减再 IFFT 回时域解扩。门限的设定是核心。设 FFT 点数为 ( N )噪声功率谱密度均匀时谱线幅度服从瑞利分布。门限 ( T ) 一般取噪声均值的 3 到 5 倍或者用中位数加倍数的方式自适应估计。中位数法更稳健因为窄带干扰本身会拉高均值用均值算门限容易把门限抬得过高导致该置零的没置零。3.2 频域置零的完整代码与参数说明function rx_out freq_notch(rx, Nfft, threshold_factor) % 频域窄带干扰抑制重叠保留法 门限置零 % rx: 输入时域信号 % Nfft: FFT 点数 % threshold_factor: 门限倍数典型 3~5 L length(rx); step Nfft / 2; % 50% 重叠 win hann(Nfft); % 加窗减少泄漏 rx_out zeros(L, 1); count zeros(L, 1); for start 1:step:(L - Nfft 1) seg rx(start : startNfft-1) .* win; X fft(seg); mag abs(X); % 用中位数估计噪声底避免干扰拉高均值 noise_floor median(mag); threshold threshold_factor * noise_floor; % 超过门限的谱线置零 mask mag threshold; X(mask) 0; seg_out real(ifft(X)); rx_out(start : startNfft-1) rx_out(start : startNfft-1) seg_out; count(start : startNfft-1) count(start : startNfft-1) win; end % 归一化重叠部分 count(count eps) 1; rx_out rx_out ./ count; end这段代码用的是重叠保留法Nfft一般取 256 或 512threshold_factor从 3 开始试。加 Hann 窗是为了减少频谱泄漏否则强干扰的旁瓣会被误判成多个干扰峰。置零之后一定要做重叠归一化不然窗函数的形状会调制信号包络解扩时反而引入额外误码。参数上最需要调的是threshold_factor。设太小噪声谱线的随机起伏会被误置零相当于在信号上开了很多小口子设太大弱一点的窄带干扰就漏过去了。我的经验是先在无干扰条件下跑一遍看纯噪声时误置零的比例控制在 1% 以下再把这个倍数用到有干扰的场景。3.3 自适应陷波器的替代方案与适用边界频域置零的缺点是它会同时干掉干扰频点上的有用信号分量。如果窄带干扰正好落在信号能量集中的频点置零的代价就很大。这时候可以用自适应陷波器ANF它只对干扰频点做窄带衰减保留其余频谱。LMS 自适应陷波器的核心是一个二阶 IIR 陷波中心频率由 LMS 自适应调整。MATLAB 里可以用adaptfilt或者手写更新公式。手写的话权值更新步长 ( \mu ) 要满足 ( 0 \mu 1/(N \cdot P_x) )其中 ( P_x ) 是输入功率。步长太大陷波器会震荡太小则跟踪不上干扰频率的漂移。实际选型时如果干扰频率固定频域置零更简单直接如果干扰频率有慢漂移或者同时存在多个窄带干扰自适应陷波更合适。但自适应方案的收敛时间通常在几百到几千个采样点对突发干扰的响应不如频域置零快。4. 避坑与排查DSSS 抗窄带干扰仿真里最容易翻车的五件事4.1 解扩后 BER 不降反升现象加了频域置零之后BER 比不加还高。原因重叠保留法的归一化没做对或者 FFT 点数选得太小导致置零时把信号主瓣也削掉了。解决先检查count归一化那一步确保每个采样点的窗函数累加和正确。然后把Nfft从 256 加到 1024 再试如果 BER 改善说明是频域分辨率不够。最后确认干扰带宽和 FFT bin 宽度的关系干扰带宽至少要覆盖 3 个 bin 以上否则置零会误伤信号。4.2 处理增益算出来和仿真对不上现象理论处理增益 23 dB但仿真里 JSR 到 15 dB 就崩了。原因PN 码的周期太短或者扩频因子不是整数导致码片和比特没对齐。解决检查sf Rc/Rb是否为整数不是整数就调整Rc或Rb。另外 PN 码周期要大于一帧数据长度否则重复的码型会产生周期性相关峰等效于降低了处理增益。我一般会让 PN 码周期至少是Nbit * sf的两倍。4.3 窄带干扰模型不像窄带现象频谱上看干扰占了好几个 MHz根本不是一根针。原因butter带通滤波器的过渡带太宽或者Bj设得太大。解决把Bj降到Rc的 1% 到 5%同时提高滤波器阶数到 6 阶但要注意相位失真。更好的做法是用正弦波加窄带噪声的混合模型正弦波代表单频干扰窄带噪声代表部分带宽干扰两者分开测试。4.4 误码统计在低 BER 时全是零现象BER 曲线在 ( 10^{-3} ) 以下直接掉到零看不出趋势。原因仿真比特数不够2000 个比特在 BER 为 ( 10^{-4} ) 时期望错误数只有 0.2 个。解决用半解析方法只对噪声做蒙特卡洛干扰用解析计算。或者把Nbit加到 ( 10^6 )但要注意内存和运行时间。折中方案是分段统计先跑 ( 10^4 ) 个比特看趋势再在关键点加密。4.5 固定随机种子导致结论不可信现象换一台机器跑BER 曲线完全不一样。原因rng种子固定但不同 MATLAB 版本的随机数生成器实现有差异。解决不要依赖固定种子来复现结果而是每个参数点跑多次取平均。种子只用来保证单次调试时可复现最终结论必须基于统计平均。我一般会在脚本里把rng(shuffle)和固定种子做成可切换的选项调试时固定出图时打乱。5. 把仿真推到边界用半解析法快速定位抗干扰门限跑完整蒙特卡洛虽然直观但要在 ( 10^{-5} ) 的 BER 上扫参数计算量会大到让人想砸键盘。我后来改用半解析法把解扩后的判决变量写成信号 干扰残余 噪声的形式其中干扰残余是确定性的噪声是高斯性的这样 BER 可以直接用 Q 函数算出来不用真的跑误码统计。具体做法是对每个 JSR 点只跑一次链路得到解扩后的干扰残余波形然后把它和理论噪声方差一起代入 Q 函数。这样每个点只需要一次 FFT 和一次解扩速度比蒙特卡洛快两个数量级。下面是对应的核心代码function ber semi_analytic_ber(JSR, SNR, Bj) % 半解析 BER干扰残余确定性计算噪声用 Q 函数 % 跑一次链路拿到解扩后的干扰残余 [~, intf_residual, sig_pow] run_dsss_once(JSR, Bj); % 噪声方差解扩后 noise_var sig_pow / 10^(SNR/10) / sf; % 判决变量 信号幅度 ± 干扰残余 噪声 % 对每个比特干扰残余不同逐比特算 BER 再平均 Nbit length(intf_residual); ber_per_bit zeros(Nbit, 1); for k 1:Nbit mu 1 intf_residual(k); % 归一化信号幅度为 1 ber_per_bit(k) qfunc(mu / sqrt(noise_var)); end ber mean(ber_per_bit); end这里qfunc是 MATLAB 通信工具箱里的 Q 函数没有工具箱的话用0.5*erfc(x/sqrt(2))代替。intf_residual是解扩后每个比特上的干扰残余它已经包含了扩频码和干扰波形的相关结果。这个方法的精度在干扰为高斯或类高斯时很好但如果干扰是单频正弦残余会有周期性需要多跑几个相位取平均。用这个方法我把 JSR 扫描从原来的半小时压缩到两分钟才有力气去扫干扰带宽和频点位置的二维网格。后来发现最坏情况不是干扰在中心频点而是在 ( \pm 0.3B_s ) 附近因为那里信号频谱的滚降还没完全起来干扰残余和信号主瓣的耦合最强。这个结论如果靠全蒙特卡洛跑我可能到现在还没扫完。最后说个习惯每次改完抗干扰模块我都会先把 JSR 设成 0 跑一遍确认 BER 和纯噪声理论值对得上再逐步加干扰。这一步花不了两分钟但能挡掉八成以上的算法无效误判。希望帮到你。本文还有配套的精品资源点击获取