简介资源主要面向雷达信号处理与FMCW雷达回波杂波抑制应用提供基于CLEAN算法的频域杂波抑制仿真实现。作者针对真实雷达回波数据设计了相减函数cos形式并预留了改为指数形式的扩展思路适用于对频域杂波抑制算法进行仿真验证或参数调整的工程师与科研人员。压缩包共1个文件包含一个.m源文件大小仅1KB代码精简便于快速阅读和修改参数后直接运行。该资源已有895人学习下载虽然体量小巧但聚焦CLEAN算法在雷达杂波抑制中的核心实现逻辑能够帮助读者理解频域对消处理的流程并可迁移至其他相似频域杂波抑制场景作为算法初探或功能验证的参考脚本。1. 雷达回波数据里的杂波CLEAN 算法为什么值得一试处理真实FMCW雷达回波数据时最头痛的不是目标太弱而是墙体、地面和大型静止物体形成的杂波太强。一次实测中一个静止金属板的回波把后方的行人完全压在噪声里经过MTD之后依然残留一条横向拖尾。频域CLEAN算法就是用来对付这类强杂波的它把信号频谱看成“若干强谱峰加连续弱分量”的组合逐次找出最大峰值、估计幅度和相位再从时域信号里减去一个对应的余弦分量迭代几轮后杂波被剥离目标才能露出来。这个思路最早来自射电天文图像处理但用在雷达杂波仿真和真实回波预处理上同样顺手。注意这套仿真基于相减函数为 cos() 的实信号假设换成复指数就能直接处理复数回波。适合雷达信号处理工程师和做检测前数据清洗的人。2. 频域CLEAN算法原理与FMCW杂波建模2.1 从射电天文“脏图”到雷达频域去杂波CLEAN算法最初解决的是射电天文里的“脏图”问题望远镜阵列得到的图像中真实源周围会叠加由稀疏基线产生的旁瓣整体反卷积会放大噪声。CLEAN不做整体逆滤波而是把图像分解成“源 点扩散函数”从脏图里一次次找到最大峰值再减去该源对应的点扩散函数贡献最后把所有源重新摆放得到接近无旁瓣的干净图。把CLEAN搬到雷达频域对象从二维图像变成一维频谱或者距离-慢时间矩阵。静止杂波在慢时间维上是强窄带谱峰相当于“脏谱”里的强源微弱目标多普勒展宽较宽幅度低相当于背景弱源。如果一次用陷波器把所有强频点全部挖掉目标分量也会跟着丢失。CLEAN的迭代消去更谨慎每轮只消掉当前最强峰而且消掉的是由该峰估计出的时域余弦分量不会动其他频段这样能保留目标较宽的频谱细节不像零频陷波那样一刀切。从工程参数上看CLEAN还有一个经典循环增益概念并不是每轮必须减掉整个峰而是只减掉峰值幅度的一部分。常见做法是设增益 0.5~1增益越小收敛越慢但对强杂波旁瓣的抑制越稳定。源码里没有把这个增益单独暴露成参数我习惯通过 threshold 和 max_iter 来控制等效剥离量后面第四章会具体说明。2.2 FMCW回波谱结构与杂波可稀疏性假设FMCW雷达发射线性调频连续波每个chirp的回波经过混频后得到差拍信号差拍频率正比于目标距离所以对单个chirp做FFT得到距离维。一次测量包含M个chirp取同一个距离门的M个采样点做慢时间FFT就得到该距离门的多普勒谱。地面、护栏这类静止目标的回波在慢时间频谱上集中在零频附近摆动的树木、旋转机械等形成特定窄带谱线。杂波源的数量远小于距离-多普勒单元总数因此具有明显的稀疏性这正是CLEAN能用的前提。下面这张表可以帮助判断什么场景值得用CLEAN什么场景应该换方法杂波类型慢时间谱特征CLEAN可处理性静止地物、墙体零频附近强谱峰很好谱峰少且集中旋转风扇、电机单一或多条窄带好只要带宽远小于频谱分辨率移动多人编队多普勒展宽谱峰较宽一般需要提高频率分辨率大面积植被、海浪宽带有色噪声不建议迭代会把噪声抽成假峰需要强调CLEAN并不区分强峰是目标还是杂波它只保证“把能量最强的一批窄带分量剥离”。如果目标本身也是窄带强回波就必须靠 threshold 控制迭代深度让目标峰值不进入剥离集合。2.3 默认工程里 cos() 相减函数到底在减什么摘要里提到“相减函数是 cos()如果需要复数形式改成指数就行了”这是理解这套源码的关键。假设输入是一段实信号采样一个窄带杂波在时域表现为A*cos(2*pi*f0*t phi)。CLEAN在频域找到它的正频率峰值后先还原幅度、频率和相位然后生成同样的余弦波从原信号中减去。下面这段示意代码正是这个还原过程% 频域峰值剥离在正频域找峰值反推时域余弦并减去 N numel(sig); spec fft(sig); pos_range 2:floor(N/2)1; [~, idx_local] max(abs(spec(pos_range))); idx pos_range(idx_local); Amp 2 * abs(spec(idx)) / N; % 实信号时域幅度换算 Phase angle(spec(idx)); % 峰值相位 f0 (idx - 1) / N * fs; % 峰值频率 sig_clean sig - Amp * cos(2*pi*f0*t Phase);代码里的要点有三个一是idx从2开始跳过直流分量直流通常由高通滤波器或后续杂波减法一并处理二是Amp必须乘以2再除以N因为一个实余弦的FFT能量被平均分到了正负两个频率三是Phase直接取自频谱峰值处的相位不要额外加pi/2。实数信号的余弦减法会同时消掉正负频率两条谱线因此不会产生单边泄露。如果换成复数基带数据比如相干接收机输出的I/Q信号杂波不再是实余弦而是复指数A*exp(1i*(2*pi*f0*t phi))频谱只有一个单边峰幅度换算要去掉2倍因子相减代码也要改成指数形式。这不只是把cos换成exp的机械替换还涉及对负频率是否存在判断。动手前先用isreal(sig)确认数据类型避免把实信号误判成复数。3. main.m 工程拆解读回波、做距离维 FFT、循环CLEAN3.1 解压后的文件结构与数据流源码包名称是main_雷达回波数据_雷达杂波_雷达CLEAN算法_雷达杂波仿真_雷达_源码.zip解压后最关键的文件是main.m。第一次使用这套代码的人容易按C语言习惯去找main函数的参数类型实际上MATLAB里的 main.m 就是一个普通脚本可以直接运行。工程的数据流可以按四条线理解先读入真实雷达回波数据通常是一个矩阵行是chirp内的ADC采样点数列是chirp编号然后对每一列做距离维FFT得到距离-慢时间矩阵接着对矩阵的每一行也就是每个距离单元做CLEAN迭代消去慢时间维上的强杂波谱峰最后把清理后的数据重新排列输出给后续CFAR检测或参数提取。参数集中在 main.m 顶部的结构体param中常用参数如下参数名含义推荐初始值param.fs差拍信号采样率Hz2e6param.prfchirp重复频率Hz1000param.num_chirp慢时间chirp数量256param.max_iterCLEAN最大迭代次数150param.threshold峰值幅度阈值8×噪声中位数param.window慢时间窗函数hamming(N)真实回波数据不一定是矩阵形式。有些录波文件把I/Q分成两个通道读入时要x qi(1:2:end) 1i*qi(2:2:end)也有的设备已经做过距离FFT直接输出[距离门×chirp]复数矩阵。这时 main.m 里应该跳过距离维FFT直接把数据送入CLEAN。我一般先用size()和isreal()打印工程结构确认数据格式后再决定是否走距离FFT避免把原始差拍信号当成距离数据。3.2 距离-慢时间数据组织与预处理代码对于标准FMCW差拍数据第一步是沿chirp维度做距离FFT。该步骤把时域差拍频率映射成距离窗函数和采样率决定距离分辨率。下面是一个通用距离维预处理函数function slow_time range_processing(adc_data, fs, range_win) % adc_data: [ADC采样点数, chirp数]每一列是一个chirp的差拍信号 % fs: 差拍采样率 % range_win: 距离窗向量长度与adc_data行数一致 n_samples size(adc_data, 1); range_fft fft(adc_data .* range_win, n_samples, 1); slow_time range_fft.; % 转置成 [range_bins, chirps] end在 main.m 中调用load(real_echo.mat, adc_data); param.fs 2e6; range_win hamming(size(adc_data, 1)); slow_time range_processing(adc_data, param.fs, range_win); % 对慢时间维加窗抑制多普勒谱旁瓣 slow_time slow_time .* hamming(size(slow_time, 2));第二行代码里hamming生成列向量长度与adc_data行数一致用.*按列广播。fft第三个参数1表示沿第一维做FFT也就是每个chirp内做距离变换。如果不写1MATLAB默认按列处理在矩阵不是方阵时容易出错。转置成[range_bins, chirps]是刻意的因为CLEAN循环希望每个距离单元的一维序列放在行上方便用for遍历。3.3 CLEAN主迭代循环实现CLEAN主循环是工程的核心。下面的clean_1d函数对一维慢时间序列做杂波剥离输入是加窗后的信号输出仍在加窗域function sig_out clean_1d(sig_in, param) N numel(sig_in); t (0:N-1) / param.prf; spec fft(sig_in); sig_out sig_in; for iter 1:param.max_iter pos_range 2:floor(N/2)1; [peak, idx_local] max(abs(spec(pos_range))); if peak param.threshold break; end idx pos_range(idx_local); amp 2 * abs(spec(idx)) / N; % 实信号余弦幅度 phase angle(spec(idx)); f0 (idx-1) / N * param.prf; comp amp * cos(2*pi*f0*t phase); sig_out sig_out - comp; spec fft(sig_out); end end这个函数每轮做五件事在当前频谱的正频率区域找最大峰值判断峰值是否低于阈值若低于阈值直接停止否则计算频率、幅度和相位生成同频同幅同相的余弦波从信号中减去该余弦并重新FFT。由于spec每次都由sig_out重算后续迭代能真实反映残留峰而不是依赖上一轮的残余输出。pos_range从2开始到floor(N/2)1确保负频率不参与峰值索引也避免直流偏置抬高门限。调用时按距离单元循环x_clean zeros(size(slow_time)); for rb 1:size(slow_time, 1) x_clean(rb, :) clean_1d(slow_time(rb, :), param); end这段调用代码要注意slow_time(rb,:)是复数还是实数。如果数据是实中频cos相减没问题如果已经是IQ复数要改成指数形式详见第四章。另一个问题是因为输入已经加窗CLEAN输出幅度也被窗函数调制过。如果后续要做多普勒积累或能量估计需要在外层把窗函数补偿回来常见处理是除以窗函数均值或者在做CFAR前用x_clean ./ hamming(size(slow_time,2))恢复幅度。4. 调阈值、控迭代、换相减模型把仿真调到真实数据可用4.1 峰值幅度与噪声底threshold 参数怎么给threshold 是CLEAN里最敏感的参数控制迭代何时停止。定得太低CLEAN会把噪声底上的随机尖峰也当成杂波剥离削弱真实目标的宽带能量定得太高强杂波的剩余旁瓣又会进入后续CFAR造成虚警。我一般不用绝对电压值而是根据当前信号的噪声中位数动态计算spec_abs abs(fft(slow_time(1, :))); noise_floor median(spec_abs); param.threshold 8 * noise_floor; param.threshold min(param.threshold, max(spec_abs) / 20);第一行取中位数而不是均值是为了避免少数强杂波峰把噪声基底估计拉高第二行做上限约束确保threshold不会超过最强峰的二十分之一。这里的8和20在空气清洁但地面反射严重的停车场场景比较稳定到了植被茂密场景需要重新标定。可选做法是先跑一轮不设threshold的CLEAN把每轮峰值幅度打印出来观察“峰值急剧下降”的拐点再取拐点对应的幅度作为threshold。前几轮峰值下降快后面下降变慢继续迭代只是在抽取噪声threshold应卡在拐点附近。调试时还可以把每一轮估计出的[f0, amp]记录成数组绘制一张“杂波谱线散点图”。如果发现少数频率被反复估计说明窗函数旁瓣泄漏正在干扰迭代应该先回到窗函数参数上。4.2 迭代次数与窗函数谱泄漏才是隐藏的坑max_iter 不能随便设置。迭代次数过少强窄带杂波剥离不干净过多则会进入“挖旁瓣”状态。当一个强峰不在FFT整数频率上时能量会分散到相邻几个binCLEAN每轮减掉一个cos后剩余旁瓣在下轮形成残留峰频率落在附近看起来就像迭代曲线周期跳动这对应射电天文里的“火车轮效应”。缓解办法首先是加窗。矩形窗旁瓣只有-13dBHamming窗可以压到-43dB左右对CLEAN稳定性提升非常明显。下面是我常用的一组对照表窗函数主瓣宽度最高旁瓣对CLEAN迭代的影响矩形2 bin-13 dB容易产生残峰需更细致阈值Hamming4 bin-43 dB推荐主瓣稍宽但旁瓣抑制好Hann4 bin-31 dB和Hamming接近频率估计偏置稍大Blackman6 bin-58 dB旁瓣最小但幅度动态损失明显在 main.m 中加窗很简单调用CLEAN前执行slow_time slow_time .* hamming(size(slow_time,2));。加了Hamming窗后param.max_iter设到150~200通常足够。由于 threshold 已经提供提前退出机制max_iter更像保护上限不是精确循环次数。clean_1d的输出仍在加窗域我通常在之后用x_clean ./ hamming(size(x_clean,2))恢复幅度再计算信杂比。4.3 从 cos() 到 exp(iθ)复数雷达数据的正确切法如果雷达接收机输出I/Q复基带信号slow_time本身是复数。此时CLEAN的相减函数要从cos()改成exp(1i*...)幅度换算系数也要变。复数指数形式如下amp abs(spec(idx)) / N; % 复指数不需要乘2 comp amp * exp(1i * (2*pi*f0*t phase)); sig_out sig_out - comp;这里的phase仍然是angle(spec(idx))不要因为复信号就额外加pi/2。对一个复指数A*exp(1i*(2*pi*f0*tphi))做FFT峰值bin上的频谱值是A*N所以除以N就能还原幅度实余弦乘2是因为负频占掉一半能量。工程中改成复信号后还要多验证一步连续迭代后各峰相位是否稳定。复信号的杂波通常只有单边谱CLEAN减去分量不会影响负频率但如果数据里有直流偏置或I/Q不平衡负频率也会出现伪峰。这时要先做直流消除和I/Q校准再进CLEAN。简单检查方法是mean(sig)复信号直流分量应接近0明显不为0就先在时域减去均值否则CLEAN会在零频附近反复估计白白消耗迭代次数。反过来如果数据是实信号却直接换成指数降噪正负频能量只从单边估计幅度反而不准所以isreal那一句检查不能省。5. 验证技巧把 CLEAN 输出送进 CFAR 才算闭环5.1 用合成回波验证CLEAN行为在拿真实数据跑通之前我会先构造一个最小合成信号验证CLEAN剥离精度。下面是一段调试脚本t (0:255) / 1000; sig 5*cos(2*pi*50*t) 0.5*cos(2*pi*120*t 1) 0.05*sin(2*pi*80*t); param.prf 1000; param.max_iter 30; param.threshold 0.2; sig_clean clean_1d(sig, param);前两个cos是模拟强杂波幅度为5和0.5第三个正弦是幅度0.05的弱目标。跑完后画abs(fft(sig))和abs(fft(sig_clean))观察50Hz和120Hz谱峰是否被压到-40dB以下80Hz目标是否保留。如果50Hz还有明显残峰首先怀疑频率落在栅栏中间需要增大FFT点数或用Goertzel做精化频率估计而不是盲目增加迭代次数。5.2 CLEAN 结果进入 CFAR 的接口CLEAN输出往往带窗函数伪幅度直接送CFAR会让门限偏低。我习惯在循环结束后执行一次逐单元补偿x_clean x_clean ./ hamming(size(slow_time, 2)); det_map cfar2d(abs(x_clean).^2, 8, 16, 1e-4);这里的cfar2d表示二维CFAR函数常见参数是保护单元8、训练单元16、虚警率1e-4。如果工程里没有现成CFAR也可以用一维OS-CFAR代替核心目的是验证杂波抑制后信杂比提升。信杂比改善因子可以这样算SCV_in max(abs(slow_time).^2) / median(abs(slow_time).^2); SCV_out max(abs(x_clean).^2) / median(abs(x_clean).^2); IF 10*log10(SCV_out / SCV_in);当IF大于10dB说明CLEAN确实压下了杂波能量如果接近0dB甚至为负多半是threshold或窗函数参数不对把目标能量连带削掉。用这个数值指导param.threshold和param.max_iter的取舍比只靠肉眼观察频谱更可靠。本文还有配套的精品资源点击获取