ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

白噪声与有色噪声的频谱特征:功率谱密度仿真分析指南

白噪声与有色噪声的频谱特征:功率谱密度仿真分析指南 做信号处理仿真的朋友大概率都跟噪声打过交道。把一段噪声信号扔进傅里叶变换得到的频谱是平平一条底噪还是低频端明显翘起来的斜坡直接决定你在频谱分析里能不能看清真实信号。这是“傅里叶变换与频谱分析”系列的第九篇专门把噪声信号的频谱特征讲透白噪声、有色噪声怎么建模怎么用仿真软件求出它们的功率谱密度以及拿到频谱图之后如何反推问题。不管你是学生、刚入行的算法工程师还是做雷达信号处理、通信系统验证的老手这套噪声仿真的思路都值得花半小时完整跑一遍。1. 噪声的数学模型先搞懂“白”和“有色”再动手1.1 白噪声为什么“白”它在频域里有多平白噪声的“白”是从白光借来的说法。白光的频谱里各波长能量均匀而白噪声在各个频率上的功率密度也基本是常数。数学上最常用的是高斯白噪声即每个采样时刻的幅度独立地从正态分布中抽取均值通常记作0方差记作σ²。离散域里它的功率谱密度近似为常数σ²/fsfs是采样率也就是说把频谱画出来应该是一条随频率变化的水平线。这里有个关键点理想白噪声在数学上并不存在因为它要求带宽无限大能量无限大。工程中说的“白噪声”都是带限白噪声只要在关心的频带内功率谱平坦就可以当成白噪声处理。ADC量化噪声就是一个非常典型的带限白噪声例子量化误差在奈奎斯特频带内近似均匀分布频谱平坦得让人放心。理解白噪声频谱为什么平要从自相关函数看。白噪声序列各点之间互不相关自相关函数在零点是一个冲击其余位置都是零。维纳-辛钦定理告诉我们功率谱密度等于自相关函数的傅里叶变换冲击对应的频谱就是常数。所以“白噪声频谱平”这件事不是靠直觉猜的而是有严格的数学推导兜底。1.2 有色噪声粉红噪声和布朗噪声为什么低频堆能量与白噪声相对的是有色噪声名字同样来自光学类比意思是不同频率上的“能量颜色”不均匀。最常见的是粉红噪声功率谱密度大约与频率成反比S(f) ∝ 1/f画在对数坐标里是一条每十倍频程下降3dB的斜线。更“重”一点的是布朗噪声功率谱密度与频率的平方成反比S(f) ∝ 1/f²每十倍频程下降6dB。粉红噪声在工程里非常常见半导体器件的闪烁噪声也就是1/f噪声、机械系统的低频漂移、部分传感器输出的慢变扰动都带有明显的1/f特征。布朗噪声则更像是“随机游走”的结果典型例子包括陀螺仪的零偏漂移、时钟信号的长期相位漂移。这两类噪声放到频域看低频端会明显抬升严重的时候直接把真实信号淹没在低频斜坡里。我最早做传感器信号处理时吃过这个亏拿一段带漂移的加速度计数据直接做FFT低频那一大片抬升还以为是机械振动后来换了对数坐标一看斜率接近-6dB/oct才意识到是布朗噪声在作怪。所以看见低频异常隆起先别急着怀疑是真实信号要想想噪声的颜色对不对。1.3 时域随机性到频域特征的桥时域里“随机”和频域里“频谱形状”之间有一座桥就是自相关函数。白噪声自相关短当前值跟过去未来都不沾边所以频谱平有色噪声自相关长当前值受历史影响大低频成分自然就多。可以用一个生活化的类比理解自相关函数反映的是随机过程对“过去”的记忆长度。白噪声像失忆的人每一秒都是全新的布朗噪声像记性特别好的人每一刻的状态都是过去所有状态累加出来的。记忆越长变化越慢低频能量越集中。仿真噪声信号的时候脑子里装着这张“记忆-频谱”对应图后面看到任何奇怪的谱形都不会心里发虚。2. 工具选型与参数设计噪声仿真链路怎么搭最稳2.1 MATLAB还是Python仿真链路怎么选做噪声信号仿真目前主流还是MATLAB和Python二选一。MATLAB的优势在于Signal Processing Toolbox里现成函数非常全randn生成高斯白噪声periodogram算周期图pwelch做Welch平均谱dsp.ColoredNoise甚至可以直接生成指定颜色的噪声序列对快速验证想法极其友好。Python的优势在于免费、开源numpy生成随机序列scipy.signal里的welch和firwin完全能覆盖MATLAB的功能而且算法后期嵌入业务系统更方便。很多人问Multisim、SPICE这类电路仿真软件能不能做噪声分析我的看法是它们和这里讨论的不是一个层次。电路仿真关心的是电阻热噪声、晶体管沟道噪声在具体电路里的传递属于器件级和电路级仿真而本文讨论的是把噪声作为随机信号在算法层面做频谱特征分析。两者都是“噪声仿真”但对象和方法完全不同别混在一起选型。我的习惯是教学和方案论证用MATLAB因为画图方便、交互顺手一旦需要批量跑数据和做算法移植就切到Python。这不是谁替代谁的问题而是两个阶段的工具各有长短。2.2 采样率和采样点数仿真的第一个坎噪声仿真的频率分辨率由采样率和点数共同决定分辨率Δf fs/N。举个例子采样率fs取1000Hz点数N取10000那么频率分辨率是0.1Hz能看清0到500Hz范围内的频谱细节。如果只取1024点分辨率大约0.98Hz低频细节基本糊成一团。随机信号的频谱还有一个特点观测时间越短谱线波动越大。这是因为随机过程单次实现受到样本波动影响数据越短估计方差越大。所以噪声仿真我通常建议N取8192以上想要平滑的谱线就用Welch平均法把长数据切成多段再平均。这里补一个实操经验采样率不要刚好卡在信号最高频率的两倍。理论上奈奎斯特条件允许但真实频谱会有过渡带和泄漏留个20%到50%的余量会让结果舒服很多。想观察50Hz工频干扰附近的情况我会把fs定在500Hz到1kHz之间既保证细节又不会让数据量失控。3. 核心实操白噪声频谱分析的完整复现链路3.1 仿真参数与信号模型设计这一节用一个具体例子把整个链路跑通。参数定得保守但实用采样率fs 1000Hz随机序列长度N 10000随机种子固定为42生成高斯白噪声的标准差取0.5。选标准差0.5是为了让信号幅度在±1.5V左右跟很多采集卡输入量程对得上方便后续套用真实场景。生成白噪声之后第一步是去均值。虽然randn理论均值是0但实际一帧数据的样本均值往往不是严格0直流分量会在频谱的0频处形成一个巨大的尖峰对数坐标下会压制低频段的显示效果。去均值这个动作虽小却能避免后续读谱时被直流泄漏干扰判断。3.2 MATLAB实现从randn到功率谱密度用MATLAB实现的代码非常短fs 1000; % 采样率 1000 Hz N 10000; % 样本点数 rng(42); % 固定随机种子保证可复现 x randn(1, N) * 0.5; % 高斯白噪声标准差0.5 x x - mean(x); % 去直流 % 周期图法估计功率谱密度 [pxx, f] periodogram(x, [], N, fs); plot(f, 10*log10(pxx), LineWidth, 1); xlabel(Frequency (Hz)); ylabel(Power/Frequency (dB/Hz)); grid on;代码背后有几件事需要解释。periodogram默认使用矩形窗相当于直接对整段数据做FFT再取模平方除以N和fs得到PSD。随机噪声的周期图会有明显的毛刺这是正常的方差表现不代表频谱有问题。如果觉得谱线太晃眼可以把periodogram换成pwelch把数据分段求平均曲线会平滑很多。实测下来这段代码输出的PSD会在以-3dB/Hz附近波动波动范围大概±2dB整个0到500Hz范围内没有系统性倾斜。这就是白噪声在频域应该有的样子。注意对数纵坐标非常关键噪声的功率谱动态范围动辄几十dB线性坐标根本看不清楚全貌。3.3 叠加正弦信号用SNR控制谱峰与底噪的高度差单独看白噪声只是热身工程里真正常用的是“正弦信号白噪声”模型用来模拟接收机中信号叠加底噪的场景。先解释一下信噪比SNR的定义信号功率除以噪声功率。假设噪声标准差是0.5噪声功率就是0.25。要让SNR等于20dB即功率比100倍信号功率应该是25正弦信号幅值就是sqrt(2)*5≈7.07。fs 1000; N 10000; rng(42); noise 0.5 * randn(1, N); t (0:N-1) / fs; snr_target 20; % 目标SNR 20dB signal_amp sqrt(2 * 0.25 * 10^(snr_target/10)); x signal_amp * sin(2*pi*50*t) noise; x x - mean(x); [pxx, f] periodogram(x, [], N, fs); plot(f, 10*log10(pxx)); xlabel(Frequency (Hz)); ylabel(Power/Frequency (dB/Hz));看这个频谱时50Hz处会出现一根明显的谱线高度比底噪高出大约20dB。这其实提供了一种从频谱反推SNR的办法读出谱峰与附近底噪的差值再考虑等效噪声带宽修正就能估算接收链路的信噪比。这个方法在实测信号分析里非常实用后面第5章会详细说修正细节。3.4 用Python复现同一帧噪声的频谱如果偏好Python代码同样不复杂。只需要numpy生成数据scipy.signal.welch计算谱。import numpy as np from scipy.signal import welch fs 1000 N 10000 rng np.random.default_rng(42) noise 0.5 * rng.standard_normal(N) t np.arange(N) / fs snr_target 20 signal_amp np.sqrt(2 * 0.25 * 10**(snr_target/10)) x signal_amp * np.sin(2*np.pi*50*t) noise x x - np.mean(x) f, pxx welch(x, fs, nperseg1024, noverlap512) plt.semilogy(f, pxx)两套工具基于同一组随机种子能得到几乎一致的结果。不同之处在于welch默认就做了分段平均谱线比periodogram平滑很多。这说明工具选择不影响物理结论影响的是观感和估计方差。我通常会用MATLAB快出图用Python做数据批处理两者配合效率最高。4. 进阶实操有色噪声仿真与频谱对比4.1 粉红噪声的生成滤波器逼近与工具箱方案粉红噪声没有像randn这样简单的一行随机函数可以直接生成最常用的方法是用滤波器对白噪声整形。MATLAB里有两个便捷路线一是直接用dsp.ColoredNoise系统对象指定Color为pink即可二是用一组经典的IIR滤波器系数逼近1/f谱形。这里给一组我在多个项目里用过的粉红噪声滤波器系数来自Paul Kellet的经典设计实测在10Hz到10kHz范围内能保持大约±0.5dB的1/f斜率b [0.049922035, -0.095993537, 0.050612437, -0.004408786]; a [1, -2.494956002, 2.017265875, -0.522189400]; N 10000; fs 1000; rng(42); x filter(b, a, 0.5 * randn(N, 1)); x x - mean(x); [pxx, f] periodogram(x, [], N, fs); loglog(f, pxx);用loglog坐标看这个频谱会在宽的频率范围内看到一条接近直线的斜坡斜率约-3dB/oct。这就是粉红噪声的标志。需要注意的是滤波器会有瞬态响应前几十个点的输出不能直接用我习惯去掉前200个点再分析。4.2 布朗噪声与带限噪声随机游走和滤波白噪声布朗噪声的生成方式更简单直接对白噪声做累加也就是随机游走x cumsum(randn(N, 1)); % 布朗噪声 x x - mean(x); % 去掉随机游走带来的直流漂移布朗噪声的频谱在低频端非常陡接近-6dB/oct。但由于随机游走会导致幅度不断扩散长时间序列的数值很容易溢出做频谱分析前必须做归一化。我在实战中会先用std(x)把信号缩放到合理范围再算PSD。带限噪声则是先设计一个带通FIR滤波器然后把白噪声送进去。比如用fir1设计一个100到200Hz的带通滤波器滤波后的白噪声只在带内平坦带外迅速衰减。这种噪声用于模拟窄带干扰和特定频段的射频底噪非常方便。4.3 三种噪声频谱特征对照把白噪声、粉红噪声、布朗噪声放在一起对比特征非常鲜明噪声类型PSD斜率生成方式典型工程来源白噪声0 dB/octrandn / standard_normalADC量化噪声、热噪声近似粉红噪声-3 dB/octIIR滤波器整形、dsp.ColoredNoise闪烁噪声、低频漂移布朗噪声-6 dB/octcumsum(randn)随机游走、陀螺零偏漂移看图经验同样重要。白噪声是平底粉红噪声呈现明显的低端斜坡布朗噪声低频端几乎要“上天”。如果仿真时看到低频端幅度大到淹没了整个纵轴范围八成是布朗噪声没归一化先压缩幅度再画图。5. 频谱图上到底看什么功率谱密度与关键指标5.1 为什么噪声要看PSD而不是FFT幅度谱很多新手一开始直接对噪声序列做FFT取幅度画出来发现谱线毛刺极多而且同一个噪声序列取不同长度的数据幅度值飘忽不定。这是因为随机信号的FFT幅度本身就随样本变化而且幅度没有对频率分辨率归一化物理意义不明确。噪声分析应该用功率谱密度PSD单位是功率每赫兹W/Hz或dBm/Hz。PSD的意义是单位频率宽度内携带的功率是多少。这样不同采样率、不同点数下算出来的结果可以直接比较。我对同一个白噪声序列做过验证用pwelch分两段独立计算PSD两条曲线几乎重合而直接FFT幅度则差异很大。所以凡是涉及随机信号优先用PSD说话。5.2 底噪电平、信噪比与灵敏度快速估算拿到一张横轴频率、纵轴dBm/Hz的频谱图如何快速估算整个系统的信噪比这里有一个常用公式某带宽BHz内的噪声总功率 PSD值 10*log10(B)。举个例子如果频谱底噪是-60dBm/Hz信号占用带宽1kHz那么噪声总功率是-6030-30dBm。假设信号功率是-5dBm信噪比就是25dB。这个换算在工程中价值很大。比如调试雷达信号处理板时从采集数据的频谱底噪可以反推接收链路噪声系数是否正常底噪每抬升3dB灵敏度就恶化3dB。一张频谱图配合简单加减法比任何仿真参数都更直观。5.3 等效噪声带宽与窗函数修正噪声功率谱密度计算里有一个很容易被忽略的修正量等效噪声带宽ENBW。加窗会影响单根谱线吸收噪声的能力非矩形窗会让谱线展宽吸收的噪声功率变多。ENBW的计算公式是ENBW fs * sum(w²) / (sum(w))²其中w是窗函数序列。常见窗函数的ENBW修正值以FFT bin为单位大致是矩形窗1.0汉明窗1.36汉宁窗1.5布莱克曼窗1.73。这意味着如果用汉宁窗估计底噪测得的功率谱要把每个bin的结果除以1.5才能还原真实功率密度。许多人在实测中觉得底噪比理论值高往往就是漏掉了这个修正。6. 常见问题与排查技巧实录6.1 频谱泄漏好好的谱线怎么变成了一座小山纯正弦信号经过FFT后如果频率没有恰好落在某个bin中心能量就会泄漏到相邻频率谱线变成一座带“裙边”的小山底噪也会被抬高。最直接的症状是加窗之前谱峰旁边有一串旁瓣像梳子一样加窗之后旁瓣消失但谱峰变宽。解决办法很简单根据分析目标选窗函数。要精确测幅值选平顶窗要测频率位置选汉宁窗而分析宽带噪声时矩形窗加Welch平均就可以了。我的经验是如果目标信号是离散谱线就用汉宁窗如果目标是宽带噪声或调制信号就别加太重的窗否则会扭曲真实的谱形。6.2 每次仿出来的频谱都不一样怎么复现randn每次运行都会产生不同的随机序列所以频谱图上毛刺位置总在变。这不是BUG是随机信号的固有性质。为了结果可复现必须在仿真开始前固定随机种子MATLAB里用rng(42)Python里用np.random.default_rng(42)。固定种子的好处是你写报告时的每一条谱线、每一个指标同事可以原样跑出来核对。另一个细节是去均值。哪怕只残留0.01的直流分量在0Hz处的谱峰也会比噪声底噪高几十dB对数坐标下低频段全被压扁。我的习惯是任何噪声信号在计算PSD之前都先减掉均值成本一次减法收益是整个低频段都干净。6.3 Welch分段参数该设多大Welch平均能有效降低谱估计方差但参数设置直接影响结果。nperseg太小分段多、方差低但频率分辨率差nperseg太大分辨率好了方差又上去了。工程上我一般取nperseg N/8到N/4之间重叠率50%是经典配置75%更平滑但多一倍计算量。参数推荐范围说明npersegN/8 ~ N/4折中分辨率和方差noverlap50% ~ 75%越高谱线越平滑nfft等于nperseg补零不提升真实分辨率还有一个坑用filter给白噪声整形时初始瞬态会污染前几百个点。这些点参与PSD计算会让低频段多出莫名其妙的突起。处理办法是丢弃滤波器输出的前一段数据或者用filter函数的初始状态参数正确初始化。我通常直接丢弃前200个点简单省事。6.4 实战案例一块信号处理板底噪“抬底”是怎么查出来的最后分享一个发生在身边的排查经历。有一次拿到一块信号处理板的采集数据理论底噪应该是-70dBm/Hz实测却抬高了6dB。直接看频谱全频段都高同时低频段还有明显的斜坡50Hz附近带着一圈小边带。当时第一反应是怀疑有电源噪声串进了模拟前端。我在仿真里用“50Hz正弦粉红噪声白噪声”叠加成模型跑出来的频谱形态跟实测数据几乎是一个模子印出来的。这从侧面验证了判断低频斜坡来自信号调理链路的1/f噪声或电源漂移50Hz边带来自工频耦合。后来在板卡上加了一组去耦电容、优化了模拟地走线底噪恢复到理论值附近。这次排查最值钱的体会是噪声频谱分析不是纸上谈兵。仿真里建立起来的白噪声、粉红噪声、布朗噪声形态记忆到了实测现场就是判断故障来源的“指纹库”。你先在仿真里把各种噪声的谱形都见过真遇到问题才不会被一片抬升的底噪搞得无从下手。
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进