简介本资源是一套面向数字信号处理学习者与工程师的MATLAB噪声生成工具集聚焦白噪声、布朗噪声Brown noise与粉红噪声Pink noise三类典型有色噪声的建模与仿真适用于系统建模、滤波器设计、听觉感知实验及生物医学信号合成等场景。压缩包共11个文件含6个核心MATLAB函数.m——如NoiseGenerator、ColoredNoise、TimeVariantAR等支持参数化噪声生成与时变特性模拟另含5个.mat数据文件如SampleECG1.mat、EM.mat等提供实测生理信号样本用于噪声叠加与对比验证。资源大小为2.19MB结构紧凑、即下即用。已有488人学习下载配套多个test*.m测试脚本如testNoiseGenerator.m、testSyntheticNoiseGenerator.m覆盖完整调用流程、参数说明与可视化示例便于快速理解算法原理并迁移至C#等外部平台通过MATLAB Compiler封装调用。 做信号仿真、音频处理或者随机过程分析的朋友迟早都会碰上一个需求在MATLAB里生成指定类型的噪声。标题里这个zip包我整理过一份类似的里面囊括了白噪声、粉红噪声pink noise和布朗噪声brown noise的生成程序还配了频谱验证脚本。这玩意儿看着不起眼但真到用的时候很多人栽在“怎么生成对了”和“怎么验证对了”这两件事上。这篇就把三种噪声的原理、MATLAB实现、参数踩坑一次说透你可直接拿去对照自己的代码改。1. 噪声是什么为什么你的仿真离不开它1.1 三种噪声的频谱性格差异先把概念摆正。噪声在信号处理里不是“没用的东西”而是“随机信号”的代名词。白噪声、粉红噪声、布朗噪声的核心区别在功率谱密度的斜率上。白噪声white noise功率谱密度在全频带内平坦也就是说每个频率分量贡献的功率一样。就像白光照亮所有可见光频率一样所以叫“白”。粉红噪声pink noise功率谱密度与频率成反比也就是每倍频程下降3 dB。低频能量明显高于高频听起来比白噪声“闷”一些像瀑布声、雨声。布朗噪声brown noise功率谱密度与频率的平方成反比每倍频程下降6 dB。能量集中在更低频声音听起来像沉闷的轰鸣也常被称为红噪声或随机游走噪声。这张表建议存下来写代码前先想清楚自己要哪一个。噪声类型功率谱斜率每倍频程衰减典型特征典型应用白噪声00 dB全频段均匀通信信道加噪、系统辨识激励粉红噪声1/f3 dB低频偏强音频测试、自然界背景声模拟布朗噪声1/f²6 dB能量集中低频随机游走、金融时序、物理噪声建模1.2 光会调函数不够得明白背后的数学MATLAB里生成白噪声只需要一行randn但生成粉红噪声和布朗噪声就没那么直接了因为不存在一个叫pinknoise的内置函数直接给你完整序列R2018b之后有dsp.ColoredNoise系统对象但很多老项目还在用R2016a/R2018a或者不想依赖额外工具箱。所以自己写生成器之前得先把三种噪声的数学关系吃透。白噪声的离散定义是一系列独立同分布的随机变量均值通常为0方差为常数。频域上它的功率谱密度PSD是平直的。实现最基础的是randn生成的是标准正态分布序列均值为0方差为1。粉红噪声的本质是“让白噪声通过一个衰减斜率为-3 dB/oct的滤波器”。因为功率谱密度掉得慢不能简单用一个一阶高通或者低通直接拼出来得用近似方法。布朗噪声更直接它本质上是对白噪声的累积和cumulative sum也就是随机游走。一个白噪声序列积分一次功率谱密度就变成原来的 1/f²。理解这三层关系之后你就能明白为什么 MATLAB 代码看起来五花八门——因为每种噪声的“正统做法”本来就不同不存在一个万能公式通吃三种。2. 核心算法拆解与代码实现2.1 白噪声两行代码的事但细节有讲究白噪声生成简单但很多人忽略了一个细节rand和randn的区别。rand生成的是均匀分布randn生成的是正态分布。通信系统仿真里加性高斯白噪声AWGN必须用randn因为中心极限定理决定了真实物理噪声叠加之后趋近正态分布。你要是用rand做高斯信道信噪比算出来全是错的。% 白噪声生成示例 Fs 44100; % 采样率单位Hz duration 2; % 时长单位秒 nSamples Fs * duration; whiteNoise randn(nSamples, 1); % 标准正态分布白噪声 % 验证均值应接近0方差应接近1 fprintf(均值: %.4f\n, mean(whiteNoise)); fprintf(方差: %.4f\n, var(whiteNoise));有一个细节值得说下用randn生成的时候默认种子是随机变化的每次运行结果都不一样。做算法对比实验时你需要可复现的噪声序列那就得在开头写rng(42)把随机数生成器的种子固定。rng(42); % 固定随机数种子保证结果可复现 whiteNoise randn(nSamples, 1);我个人习惯是把所有噪声生成函数写成独立m文件输入参数是“长度”和“采样率”如果需要的话输出是一列double数据这样在Simulink或外部脚本里调用特别顺手。2.2 粉红噪声从Voss-McCartney算法到FFT频域法粉红噪声是三种里最麻烦的。它的功率谱需要满足 1/f 的形态但时域上又不能简单地滤出来。生成方法有两种主流方案Voss-McCartney 算法和 FFT 频域法。Voss-McCartney 算法原理比较巧妙它用多层白噪声源叠加每一层的更新频率递减类似二进制计数器分频最终叠加结果的功率谱就接近 1/f。这种算法计算快适合实时生成。但它的频谱在低频段会有轻微波纹而且实现起来需要维护多个状态变量。function pinkNoise pinkNoiseVoss(nSamples) % Voss-McCartney 粉红噪声生成算法 % 参考实现16层白噪声源叠加 numSources 16; sources zeros(numSources, 1); pinkNoise zeros(nSamples, 1); for n 1:nSamples % 用二进制计数器决定哪些层更新 mask 1; for i 1:numSources if bitand(n, mask) ~ 0 sources(i) randn; end mask mask * 2; end pinkNoise(n) sum(sources) / numSources; end endFFT频域法的思路更直观既然我们要的粉红噪声在频域的幅度谱是 1/sqrt(f)那就直接在频域构造出这个形状加上随机相位再做逆傅里叶变换。function pinkNoise pinkNoiseFFT(nSamples, Fs) % 频域法生成粉红噪声 % 原理幅度谱设为 1/sqrt(f)相位随机化 f (0:nSamples-1) * Fs / nSamples; f(1) 1; % 避免直流分量为0设成1保证幅度 amplitudeSpectrum 1 ./ sqrt(f); amplitudeSpectrum(isinf(amplitudeSpectrum)) 0; % 处理f0附近 % 生成具有指定幅度的复数频谱 randomPhases exp(1i * 2 * pi * rand(nSamples, 1)); spectrum amplitudeSpectrum .* randomPhases; % 保证共轭对称使逆变换结果是实信号 for k 2:floor(nSamples/2) spectrum(nSamples - k 2) conj(spectrum(k)); end spectrum(1) real(spectrum(1)); % 直流分量实数 pinkNoise real(ifft(spectrum)); % 归一化到原始白噪声的方差水平 pinkNoise pinkNoise / std(pinkNoise); end注意频域法有个重要细节为了得到实信号频谱必须满足共轭对称。很多人第一次写直接把1./sqrt(f)乘个随机相位再去ifft出来的信号带虚部或者波形明显不对就是因为忘了做共轭对称处理。如果用的是 R2018b 之后的MATLAB可以直接用内置系统对象省去不少麻烦hpink dsp.ColoredNoise(Color, pink, ... SamplesPerFrame, nSamples, ... NumChannels, 1); pinkNoise hpink();但我还是建议自己实现一遍因为很多场景下你需要精确控制频谱形状内置对象不一定满足你的定制需求。2.3 布朗噪声积分白噪声但要处理漂移布朗噪声是最容易被“随便写写”搞砸的一种。它的数学定义是对白噪声做积分离散情况下就是累加% 布朗噪声生成 - 最朴素的版本 whiteNoise randn(nSamples, 1); brownNoise cumsum(whiteNoise);跑完这段代码画出波形你会发现一个大问题信号明显在漂移。因为累积和对低频没有约束直流分量可以无限累积导致信号要么冲到正很大要么跌到负很大。这不是数学问题而是物理上真实布朗粒子位移就长这样——但你在做信号处理仿真时往往需要的是方差稳定、均值大致稳定的“带限”布朗噪声。处理方法也不难生成之后减去线性趋势或者用一个一阶高通滤波器把超低频漂移滤掉。% 布朗噪声生成 - 去漂移版 whiteNoise randn(nSamples, 1); brownNoise cumsum(whiteNoise); % 去趋势减去最小二乘拟合的线性趋势 t (0:nSamples-1); p polyfit(t, brownNoise, 1); brownNoiseDetrended brownNoise - polyval(p, t); % 归一化到方差为1 brownNoiseDetrended brownNoiseDetrended / std(brownNoiseDetrended);另一种做法是直接构造滤波器。布朗噪声的频谱是 1/f²对应在离散域可以用一阶积分器近似% 一阶积分器生成布朗噪声 b 1; a [1, -0.999]; brownNoise filter(b, a, whiteNoise);这个方法的优点是不会出现线性漂移那么明显的问题但滤波器的极点非常靠近单位圆a的第二个系数接近-1数值稳定性要留意。实际上filter在处理长序列时如果极点太靠近单位圆会出现幅值异常的样本点所以如果用这个方法建议生成后做一次幅度裁剪或者标准化。3. 实操参数设置与应用场景3.1 为什么你的噪声频谱总是不达标很多朋友生成噪声后拿pwelch一画频谱发现斜率不对或者高频段翘起来。这里十有八九是pwelch的参数设置问题而不是生成算法的问题。pwelch的默认分段、窗函数、重叠率直接影响频谱估计的平滑度和频率分辨率。以粉红噪声为例正确画法是这样的% 验证粉红噪声频谱 [pinkSpectrum, freq] pwelch(pinkNoise, ... hamming(1024), ... % 窗函数建议hamming或hanning 512, ... % 重叠样本数 [], ... % 使用默认FFT点数 Fs); % 采样率 % 斜率验证log-log图上斜率应约为 -1 loglog(freq, pinkSpectrum); xlabel(频率 (Hz)); ylabel(功率谱密度); title(粉红噪声频谱验证斜率 ≈ -1); grid on;几个容易被坑的点窗函数不要用矩形窗旁瓣泄漏会让频谱高频段出现假抬高。首选 hamming 或 hanning。重叠率一般设 50%窗口长度的一半太低会导致频谱估计方差大画出来毛刺多。pwelch输出的频率轴是单边的但功率谱密度是按单边计算的所以画 loglog 图时低频段的点数会很少斜率看起来可能“掉得快”这是正常现象。如果噪声序列不够长建议先做重复实验求平均或者把序列分段多次计算再平均。还有一个常见的认知误区验证频率特性时横轴到底用线性坐标还是对数坐标白噪声用线性坐标看很平直但粉红噪声和布朗噪声必须用 log-log 坐标才能看出斜率的线性关系。你用线性坐标看粉红噪声会感觉“低频巨高、高频几乎为零”容易误判为代码写错了。3.2 实际仿真时如何校准噪声幅度生成之后的噪声序列幅度大概率不是你想要的值。拿通信系统仿真举例你想要的是“0 dBm 功率的噪声”但randn出来的序列方差是1功率就是1归一化阻抗下。这时候需要做功率校准。噪声的功率等于方差零均值情况下所以校准的思路是function scaledNoise ScaleNoiseToPower(noise, targetPower) currentPower var(noise); scalingFactor sqrt(targetPower / currentPower); scaledNoise noise * scalingFactor; end如果你的目标不是功率而是信噪比SNR那更常见的是先确定信号功率再反推噪声功率。举个例子一个正弦信号幅度为 A其功率约等于 A²/2单频信号如果要求 SNR 10 dB噪声功率就是信号功率除以 10^(SNR/10)。这个计算在仿真脚本里很常用。% 给正弦信号添加指定SNR的白噪声 Fs 1000; t 0:1/Fs:1; signalAmp 1; signal signalAmp * sin(2*pi*50*t); signalPower signalAmp^2 / 2; SNR_dB 10; noisePower signalPower / (10^(SNR_dB/10)); noise randn(size(t)); noise noise / sqrt(var(noise)) * sqrt(noisePower); received signal noise;这个例子里有个细节值得注意噪声先标准化到方差1再乘 sqrt(noisePower)而不是直接乘 noisePower。因为方差的单位是功率的单位要得到指定功率的噪声缩放因子应该是功率的平方根不是功率本身。这个错我见过不少人犯。3.3 三种噪声在典型场景中的实战用法噪声生成不是孤立的技术活它服务于具体场景。我在这里列几个自己在项目中实际用过的场景方便你对照思考。通信系统仿真中白噪声是信道噪声的基本模型。加性高斯白噪声是Channels工具箱和自写误码率仿真中最常见的背景噪声。粉红噪声则常见于音频分析和声学测试扬声器频率响应测试里用粉红噪声做激励信号比白噪声更合理因为粉红噪声在高频段的能量低不容易过载功放同时低频段能量充足能有效激励低频单元。布朗噪声在物理仿真里是个高频需求。比如模拟粒子的随机游走轨迹、股票价格的对数收益率累积过程还有光学系统里的频率噪声积分相位噪声。金融方向的朋友可能更熟悉“随机游走”这个说法布朗噪声本质就是它。做这类仿真时序列长度往往要很长比如模拟10万步的随机游走这时候cumsum的性能就很关键。生成一个多通道噪声矩阵也是个常见需求。比如做阵列信号处理时你需要一个 N×M 的噪声矩阵每一列是独立的白噪声/粉红噪声序列且各列之间不相关。这个需求只需要在生成函数外层加个循环注意每次循环重置随机种子否则各通道的噪声会完全一样。% 多通道独立粉红噪声生成 numChannels 4; nSamples 10000; Fs 1000; multiPink zeros(nSamples, numChannels); for ch 1:numChannels rng(100 ch); % 每个通道使用不同种子 multiPink(:, ch) pinkNoiseFFT(nSamples, Fs); end4. 常见问题与排查技巧实录4.1 问题一粉红噪声频谱的高频段翘起来了现象log-log 图上频率超过奈奎斯特频率的一半后PSD 开始往上翘整条曲线呈“U”型。原因大多数自实现的粉红噪声算法尤其是 Voss-McCartney 算法在高频段没有严格做到 1/f 衰减。Voss 算法叠加的白噪声源数量有限通常16层最高频层的更新频率受限于采样率导致高频段能量没被压下去出现谱隆起。解决办法改用 FFT 频域法它对频谱形状的控制更精确。如果必须用 Voss 算法增加层数到32层同时加上一个简单的一阶低通滤波器截止频率设在奈奎斯特频率的0.8倍把边缘效应压掉。在验证脚本里只看奈奎斯特频率之前 80% 的频段不要拿边缘节点说事。4.2 问题二布朗噪声的均值和方差一直在变现象每次生成的布朗噪声均值差异很大方差甚至能差出几个数量级。原因cumsum不做归一化累积和的方差是随时间线性增长的理论上布朗运动的方差正比于时间。所以长度不同的序列方差天然就不一样。这不是bug是数学特性。解决办法生成后做标准化处理减去均值、除以标准差。如果你需要多批次生成的噪声方差一致那么标准化步骤不能省。另外如果做蒙特卡洛仿真注意每一批都要单独标准化不能全批次混在一起算均值方差。4.3 问题三生成的噪声长度和预期不一致现象函数输出了 nSamples 个点但实际信号时长不是想要的。原因很多人在写生成函数时把采样率和时长搞混了。比如时长为1秒、采样率44100 Hz需要的点数是44100个但如果中途改变了采样率之前的fix(Fs * duration)计算就失效了。解决办法在函数入口处加断言直接检查长度和采样率的匹配关系assert(length(noise) round(Fs * duration), ... 噪声长度与采样率/时长不匹配);写代码时建议统一使用“采样率 时长”作为外部接口参数内部计算出样本数而不是让调用者直接传样本数这样能避免很多低级错误。4.4 问题四随机数种子设置不对导致对比实验失效现象做降噪算法对比时同一段信号加不同的噪声算法如白噪声 vs 粉红噪声但两种条件下测试结果差异很大无法判断是算法差异还是噪声差异。原因没有固定随机种子两组实验用的是不同的随机序列。虽然都是噪声但具体到某一段序列统计特性是有波动的这个波动跟真正的算法差异混在一起。解决办法在仿真脚本开头统一设置随机种子并在每次生成噪声前再设一次保证完全可复现。我自己的习惯是% 在每个实验开始前固定种子 rng(20240101 experimentIndex);这样只要 experimentIndex 不变结果就一定可复现。做完对比之后再取消固定种子做多次蒙特卡洛统计平均得到更接近真实的性能结论。4.5 问题五长序列生成时内存爆了现象生成长达几百万点的粉红噪声时MATLAB卡死或报“请求的数组超出最大数组大小”。原因FFT频域法需要构造一个跟输出长度相同的复数频谱数组再做 ifft。如果输出长度是10^7那频谱数组占用的内存就是 2 × 10^7 × 8字节double精度实部虚部各占8字节≈ 160 MB加上其他临时变量很容易冲破内存限制。解决办法优先用 Voss 算法或滤波器法生成这类方法内存占用是 O(n) 且系数小得多。如果非要频域法就分块生成。先定义总长度拆成每块 100 万点生成拼接时注意块与块之间要有重叠或交叉淡入淡出否则块边界会有咔哒声音频场景尤其明显。% 分块生成长序列粉红噪声示例 blockLength 1000000; numBlocks 10; pinkLong zeros(blockLength * numBlocks, 1); for k 1:numBlocks block pinkNoiseFFT(blockLength, Fs); pinkLong((k-1)*blockLength 1 : k*blockLength) block; end分块生成的代价是每个块之间的初始相位是独立的导致的后果是整条序列的功率谱在低频段会有细微不连续。如果你的应用对低频非常敏感那还得用重叠-相加overlap-add结构来平滑过渡。5. 关于程序包结构的建议回到最开始那个zip包如果你打算自己整理这样一个工具包我建议按下面的目录结构放拿起来就能用noise_generator/ ├── generateWhiteNoise.m ├── generatePinkNoise.m ├── generateBrownNoise.m ├── verifyNoiseSpectrum.m ├── examples/ │ ├── demo_white.m │ ├── demo_pink.m │ └── demo_brown.m └── README.mdgenerateWhiteNoise.m是最简单的封装generatePinkNoise.m里建议同时保留 Voss 和 FFT 两种实现并提供一个method参数切换。verifyNoiseSpectrum.m是验证脚本调用pwelch画频谱这是整个包里最有价值的部分——生成一堆噪声但不验证频谱特性等于白做。README 里写清楚每个函数的入参、出参、依赖工具箱其实只用基础MATLAB就行不依赖额外工具箱以及每种噪声的理论频谱斜率。这样有同行拿到你的包三分钟就能跑通。我个人在整理这类工具包时的体会是噪声生成本身不难难的是想清楚这些噪声在什么场景下用、长什么样才算对。所以这篇文章最终的落脚点不是代码而是验证意识——拿到任何一段噪声序列先画频谱、看均值、看方差确认它符合理论预期再做后续处理。这个习惯能帮你省下大量排查问题的时间。本文还有配套的精品资源点击获取 MATLAB噪声生成指南：白噪声、粉红噪声与布朗噪声的原理与实现
简介本资源是一套面向数字信号处理学习者与工程师的MATLAB噪声生成工具集聚焦白噪声、布朗噪声Brown noise与粉红噪声Pink noise三类典型有色噪声的建模与仿真适用于系统建模、滤波器设计、听觉感知实验及生物医学信号合成等场景。压缩包共11个文件含6个核心MATLAB函数.m——如NoiseGenerator、ColoredNoise、TimeVariantAR等支持参数化噪声生成与时变特性模拟另含5个.mat数据文件如SampleECG1.mat、EM.mat等提供实测生理信号样本用于噪声叠加与对比验证。资源大小为2.19MB结构紧凑、即下即用。已有488人学习下载配套多个test*.m测试脚本如testNoiseGenerator.m、testSyntheticNoiseGenerator.m覆盖完整调用流程、参数说明与可视化示例便于快速理解算法原理并迁移至C#等外部平台通过MATLAB Compiler封装调用。 做信号仿真、音频处理或者随机过程分析的朋友迟早都会碰上一个需求在MATLAB里生成指定类型的噪声。标题里这个zip包我整理过一份类似的里面囊括了白噪声、粉红噪声pink noise和布朗噪声brown noise的生成程序还配了频谱验证脚本。这玩意儿看着不起眼但真到用的时候很多人栽在“怎么生成对了”和“怎么验证对了”这两件事上。这篇就把三种噪声的原理、MATLAB实现、参数踩坑一次说透你可直接拿去对照自己的代码改。1. 噪声是什么为什么你的仿真离不开它1.1 三种噪声的频谱性格差异先把概念摆正。噪声在信号处理里不是“没用的东西”而是“随机信号”的代名词。白噪声、粉红噪声、布朗噪声的核心区别在功率谱密度的斜率上。白噪声white noise功率谱密度在全频带内平坦也就是说每个频率分量贡献的功率一样。就像白光照亮所有可见光频率一样所以叫“白”。粉红噪声pink noise功率谱密度与频率成反比也就是每倍频程下降3 dB。低频能量明显高于高频听起来比白噪声“闷”一些像瀑布声、雨声。布朗噪声brown noise功率谱密度与频率的平方成反比每倍频程下降6 dB。能量集中在更低频声音听起来像沉闷的轰鸣也常被称为红噪声或随机游走噪声。这张表建议存下来写代码前先想清楚自己要哪一个。噪声类型功率谱斜率每倍频程衰减典型特征典型应用白噪声00 dB全频段均匀通信信道加噪、系统辨识激励粉红噪声1/f3 dB低频偏强音频测试、自然界背景声模拟布朗噪声1/f²6 dB能量集中低频随机游走、金融时序、物理噪声建模1.2 光会调函数不够得明白背后的数学MATLAB里生成白噪声只需要一行randn但生成粉红噪声和布朗噪声就没那么直接了因为不存在一个叫pinknoise的内置函数直接给你完整序列R2018b之后有dsp.ColoredNoise系统对象但很多老项目还在用R2016a/R2018a或者不想依赖额外工具箱。所以自己写生成器之前得先把三种噪声的数学关系吃透。白噪声的离散定义是一系列独立同分布的随机变量均值通常为0方差为常数。频域上它的功率谱密度PSD是平直的。实现最基础的是randn生成的是标准正态分布序列均值为0方差为1。粉红噪声的本质是“让白噪声通过一个衰减斜率为-3 dB/oct的滤波器”。因为功率谱密度掉得慢不能简单用一个一阶高通或者低通直接拼出来得用近似方法。布朗噪声更直接它本质上是对白噪声的累积和cumulative sum也就是随机游走。一个白噪声序列积分一次功率谱密度就变成原来的 1/f²。理解这三层关系之后你就能明白为什么 MATLAB 代码看起来五花八门——因为每种噪声的“正统做法”本来就不同不存在一个万能公式通吃三种。2. 核心算法拆解与代码实现2.1 白噪声两行代码的事但细节有讲究白噪声生成简单但很多人忽略了一个细节rand和randn的区别。rand生成的是均匀分布randn生成的是正态分布。通信系统仿真里加性高斯白噪声AWGN必须用randn因为中心极限定理决定了真实物理噪声叠加之后趋近正态分布。你要是用rand做高斯信道信噪比算出来全是错的。% 白噪声生成示例 Fs 44100; % 采样率单位Hz duration 2; % 时长单位秒 nSamples Fs * duration; whiteNoise randn(nSamples, 1); % 标准正态分布白噪声 % 验证均值应接近0方差应接近1 fprintf(均值: %.4f\n, mean(whiteNoise)); fprintf(方差: %.4f\n, var(whiteNoise));有一个细节值得说下用randn生成的时候默认种子是随机变化的每次运行结果都不一样。做算法对比实验时你需要可复现的噪声序列那就得在开头写rng(42)把随机数生成器的种子固定。rng(42); % 固定随机数种子保证结果可复现 whiteNoise randn(nSamples, 1);我个人习惯是把所有噪声生成函数写成独立m文件输入参数是“长度”和“采样率”如果需要的话输出是一列double数据这样在Simulink或外部脚本里调用特别顺手。2.2 粉红噪声从Voss-McCartney算法到FFT频域法粉红噪声是三种里最麻烦的。它的功率谱需要满足 1/f 的形态但时域上又不能简单地滤出来。生成方法有两种主流方案Voss-McCartney 算法和 FFT 频域法。Voss-McCartney 算法原理比较巧妙它用多层白噪声源叠加每一层的更新频率递减类似二进制计数器分频最终叠加结果的功率谱就接近 1/f。这种算法计算快适合实时生成。但它的频谱在低频段会有轻微波纹而且实现起来需要维护多个状态变量。function pinkNoise pinkNoiseVoss(nSamples) % Voss-McCartney 粉红噪声生成算法 % 参考实现16层白噪声源叠加 numSources 16; sources zeros(numSources, 1); pinkNoise zeros(nSamples, 1); for n 1:nSamples % 用二进制计数器决定哪些层更新 mask 1; for i 1:numSources if bitand(n, mask) ~ 0 sources(i) randn; end mask mask * 2; end pinkNoise(n) sum(sources) / numSources; end endFFT频域法的思路更直观既然我们要的粉红噪声在频域的幅度谱是 1/sqrt(f)那就直接在频域构造出这个形状加上随机相位再做逆傅里叶变换。function pinkNoise pinkNoiseFFT(nSamples, Fs) % 频域法生成粉红噪声 % 原理幅度谱设为 1/sqrt(f)相位随机化 f (0:nSamples-1) * Fs / nSamples; f(1) 1; % 避免直流分量为0设成1保证幅度 amplitudeSpectrum 1 ./ sqrt(f); amplitudeSpectrum(isinf(amplitudeSpectrum)) 0; % 处理f0附近 % 生成具有指定幅度的复数频谱 randomPhases exp(1i * 2 * pi * rand(nSamples, 1)); spectrum amplitudeSpectrum .* randomPhases; % 保证共轭对称使逆变换结果是实信号 for k 2:floor(nSamples/2) spectrum(nSamples - k 2) conj(spectrum(k)); end spectrum(1) real(spectrum(1)); % 直流分量实数 pinkNoise real(ifft(spectrum)); % 归一化到原始白噪声的方差水平 pinkNoise pinkNoise / std(pinkNoise); end注意频域法有个重要细节为了得到实信号频谱必须满足共轭对称。很多人第一次写直接把1./sqrt(f)乘个随机相位再去ifft出来的信号带虚部或者波形明显不对就是因为忘了做共轭对称处理。如果用的是 R2018b 之后的MATLAB可以直接用内置系统对象省去不少麻烦hpink dsp.ColoredNoise(Color, pink, ... SamplesPerFrame, nSamples, ... NumChannels, 1); pinkNoise hpink();但我还是建议自己实现一遍因为很多场景下你需要精确控制频谱形状内置对象不一定满足你的定制需求。2.3 布朗噪声积分白噪声但要处理漂移布朗噪声是最容易被“随便写写”搞砸的一种。它的数学定义是对白噪声做积分离散情况下就是累加% 布朗噪声生成 - 最朴素的版本 whiteNoise randn(nSamples, 1); brownNoise cumsum(whiteNoise);跑完这段代码画出波形你会发现一个大问题信号明显在漂移。因为累积和对低频没有约束直流分量可以无限累积导致信号要么冲到正很大要么跌到负很大。这不是数学问题而是物理上真实布朗粒子位移就长这样——但你在做信号处理仿真时往往需要的是方差稳定、均值大致稳定的“带限”布朗噪声。处理方法也不难生成之后减去线性趋势或者用一个一阶高通滤波器把超低频漂移滤掉。% 布朗噪声生成 - 去漂移版 whiteNoise randn(nSamples, 1); brownNoise cumsum(whiteNoise); % 去趋势减去最小二乘拟合的线性趋势 t (0:nSamples-1); p polyfit(t, brownNoise, 1); brownNoiseDetrended brownNoise - polyval(p, t); % 归一化到方差为1 brownNoiseDetrended brownNoiseDetrended / std(brownNoiseDetrended);另一种做法是直接构造滤波器。布朗噪声的频谱是 1/f²对应在离散域可以用一阶积分器近似% 一阶积分器生成布朗噪声 b 1; a [1, -0.999]; brownNoise filter(b, a, whiteNoise);这个方法的优点是不会出现线性漂移那么明显的问题但滤波器的极点非常靠近单位圆a的第二个系数接近-1数值稳定性要留意。实际上filter在处理长序列时如果极点太靠近单位圆会出现幅值异常的样本点所以如果用这个方法建议生成后做一次幅度裁剪或者标准化。3. 实操参数设置与应用场景3.1 为什么你的噪声频谱总是不达标很多朋友生成噪声后拿pwelch一画频谱发现斜率不对或者高频段翘起来。这里十有八九是pwelch的参数设置问题而不是生成算法的问题。pwelch的默认分段、窗函数、重叠率直接影响频谱估计的平滑度和频率分辨率。以粉红噪声为例正确画法是这样的% 验证粉红噪声频谱 [pinkSpectrum, freq] pwelch(pinkNoise, ... hamming(1024), ... % 窗函数建议hamming或hanning 512, ... % 重叠样本数 [], ... % 使用默认FFT点数 Fs); % 采样率 % 斜率验证log-log图上斜率应约为 -1 loglog(freq, pinkSpectrum); xlabel(频率 (Hz)); ylabel(功率谱密度); title(粉红噪声频谱验证斜率 ≈ -1); grid on;几个容易被坑的点窗函数不要用矩形窗旁瓣泄漏会让频谱高频段出现假抬高。首选 hamming 或 hanning。重叠率一般设 50%窗口长度的一半太低会导致频谱估计方差大画出来毛刺多。pwelch输出的频率轴是单边的但功率谱密度是按单边计算的所以画 loglog 图时低频段的点数会很少斜率看起来可能“掉得快”这是正常现象。如果噪声序列不够长建议先做重复实验求平均或者把序列分段多次计算再平均。还有一个常见的认知误区验证频率特性时横轴到底用线性坐标还是对数坐标白噪声用线性坐标看很平直但粉红噪声和布朗噪声必须用 log-log 坐标才能看出斜率的线性关系。你用线性坐标看粉红噪声会感觉“低频巨高、高频几乎为零”容易误判为代码写错了。3.2 实际仿真时如何校准噪声幅度生成之后的噪声序列幅度大概率不是你想要的值。拿通信系统仿真举例你想要的是“0 dBm 功率的噪声”但randn出来的序列方差是1功率就是1归一化阻抗下。这时候需要做功率校准。噪声的功率等于方差零均值情况下所以校准的思路是function scaledNoise ScaleNoiseToPower(noise, targetPower) currentPower var(noise); scalingFactor sqrt(targetPower / currentPower); scaledNoise noise * scalingFactor; end如果你的目标不是功率而是信噪比SNR那更常见的是先确定信号功率再反推噪声功率。举个例子一个正弦信号幅度为 A其功率约等于 A²/2单频信号如果要求 SNR 10 dB噪声功率就是信号功率除以 10^(SNR/10)。这个计算在仿真脚本里很常用。% 给正弦信号添加指定SNR的白噪声 Fs 1000; t 0:1/Fs:1; signalAmp 1; signal signalAmp * sin(2*pi*50*t); signalPower signalAmp^2 / 2; SNR_dB 10; noisePower signalPower / (10^(SNR_dB/10)); noise randn(size(t)); noise noise / sqrt(var(noise)) * sqrt(noisePower); received signal noise;这个例子里有个细节值得注意噪声先标准化到方差1再乘 sqrt(noisePower)而不是直接乘 noisePower。因为方差的单位是功率的单位要得到指定功率的噪声缩放因子应该是功率的平方根不是功率本身。这个错我见过不少人犯。3.3 三种噪声在典型场景中的实战用法噪声生成不是孤立的技术活它服务于具体场景。我在这里列几个自己在项目中实际用过的场景方便你对照思考。通信系统仿真中白噪声是信道噪声的基本模型。加性高斯白噪声是Channels工具箱和自写误码率仿真中最常见的背景噪声。粉红噪声则常见于音频分析和声学测试扬声器频率响应测试里用粉红噪声做激励信号比白噪声更合理因为粉红噪声在高频段的能量低不容易过载功放同时低频段能量充足能有效激励低频单元。布朗噪声在物理仿真里是个高频需求。比如模拟粒子的随机游走轨迹、股票价格的对数收益率累积过程还有光学系统里的频率噪声积分相位噪声。金融方向的朋友可能更熟悉“随机游走”这个说法布朗噪声本质就是它。做这类仿真时序列长度往往要很长比如模拟10万步的随机游走这时候cumsum的性能就很关键。生成一个多通道噪声矩阵也是个常见需求。比如做阵列信号处理时你需要一个 N×M 的噪声矩阵每一列是独立的白噪声/粉红噪声序列且各列之间不相关。这个需求只需要在生成函数外层加个循环注意每次循环重置随机种子否则各通道的噪声会完全一样。% 多通道独立粉红噪声生成 numChannels 4; nSamples 10000; Fs 1000; multiPink zeros(nSamples, numChannels); for ch 1:numChannels rng(100 ch); % 每个通道使用不同种子 multiPink(:, ch) pinkNoiseFFT(nSamples, Fs); end4. 常见问题与排查技巧实录4.1 问题一粉红噪声频谱的高频段翘起来了现象log-log 图上频率超过奈奎斯特频率的一半后PSD 开始往上翘整条曲线呈“U”型。原因大多数自实现的粉红噪声算法尤其是 Voss-McCartney 算法在高频段没有严格做到 1/f 衰减。Voss 算法叠加的白噪声源数量有限通常16层最高频层的更新频率受限于采样率导致高频段能量没被压下去出现谱隆起。解决办法改用 FFT 频域法它对频谱形状的控制更精确。如果必须用 Voss 算法增加层数到32层同时加上一个简单的一阶低通滤波器截止频率设在奈奎斯特频率的0.8倍把边缘效应压掉。在验证脚本里只看奈奎斯特频率之前 80% 的频段不要拿边缘节点说事。4.2 问题二布朗噪声的均值和方差一直在变现象每次生成的布朗噪声均值差异很大方差甚至能差出几个数量级。原因cumsum不做归一化累积和的方差是随时间线性增长的理论上布朗运动的方差正比于时间。所以长度不同的序列方差天然就不一样。这不是bug是数学特性。解决办法生成后做标准化处理减去均值、除以标准差。如果你需要多批次生成的噪声方差一致那么标准化步骤不能省。另外如果做蒙特卡洛仿真注意每一批都要单独标准化不能全批次混在一起算均值方差。4.3 问题三生成的噪声长度和预期不一致现象函数输出了 nSamples 个点但实际信号时长不是想要的。原因很多人在写生成函数时把采样率和时长搞混了。比如时长为1秒、采样率44100 Hz需要的点数是44100个但如果中途改变了采样率之前的fix(Fs * duration)计算就失效了。解决办法在函数入口处加断言直接检查长度和采样率的匹配关系assert(length(noise) round(Fs * duration), ... 噪声长度与采样率/时长不匹配);写代码时建议统一使用“采样率 时长”作为外部接口参数内部计算出样本数而不是让调用者直接传样本数这样能避免很多低级错误。4.4 问题四随机数种子设置不对导致对比实验失效现象做降噪算法对比时同一段信号加不同的噪声算法如白噪声 vs 粉红噪声但两种条件下测试结果差异很大无法判断是算法差异还是噪声差异。原因没有固定随机种子两组实验用的是不同的随机序列。虽然都是噪声但具体到某一段序列统计特性是有波动的这个波动跟真正的算法差异混在一起。解决办法在仿真脚本开头统一设置随机种子并在每次生成噪声前再设一次保证完全可复现。我自己的习惯是% 在每个实验开始前固定种子 rng(20240101 experimentIndex);这样只要 experimentIndex 不变结果就一定可复现。做完对比之后再取消固定种子做多次蒙特卡洛统计平均得到更接近真实的性能结论。4.5 问题五长序列生成时内存爆了现象生成长达几百万点的粉红噪声时MATLAB卡死或报“请求的数组超出最大数组大小”。原因FFT频域法需要构造一个跟输出长度相同的复数频谱数组再做 ifft。如果输出长度是10^7那频谱数组占用的内存就是 2 × 10^7 × 8字节double精度实部虚部各占8字节≈ 160 MB加上其他临时变量很容易冲破内存限制。解决办法优先用 Voss 算法或滤波器法生成这类方法内存占用是 O(n) 且系数小得多。如果非要频域法就分块生成。先定义总长度拆成每块 100 万点生成拼接时注意块与块之间要有重叠或交叉淡入淡出否则块边界会有咔哒声音频场景尤其明显。% 分块生成长序列粉红噪声示例 blockLength 1000000; numBlocks 10; pinkLong zeros(blockLength * numBlocks, 1); for k 1:numBlocks block pinkNoiseFFT(blockLength, Fs); pinkLong((k-1)*blockLength 1 : k*blockLength) block; end分块生成的代价是每个块之间的初始相位是独立的导致的后果是整条序列的功率谱在低频段会有细微不连续。如果你的应用对低频非常敏感那还得用重叠-相加overlap-add结构来平滑过渡。5. 关于程序包结构的建议回到最开始那个zip包如果你打算自己整理这样一个工具包我建议按下面的目录结构放拿起来就能用noise_generator/ ├── generateWhiteNoise.m ├── generatePinkNoise.m ├── generateBrownNoise.m ├── verifyNoiseSpectrum.m ├── examples/ │ ├── demo_white.m │ ├── demo_pink.m │ └── demo_brown.m └── README.mdgenerateWhiteNoise.m是最简单的封装generatePinkNoise.m里建议同时保留 Voss 和 FFT 两种实现并提供一个method参数切换。verifyNoiseSpectrum.m是验证脚本调用pwelch画频谱这是整个包里最有价值的部分——生成一堆噪声但不验证频谱特性等于白做。README 里写清楚每个函数的入参、出参、依赖工具箱其实只用基础MATLAB就行不依赖额外工具箱以及每种噪声的理论频谱斜率。这样有同行拿到你的包三分钟就能跑通。我个人在整理这类工具包时的体会是噪声生成本身不难难的是想清楚这些噪声在什么场景下用、长什么样才算对。所以这篇文章最终的落脚点不是代码而是验证意识——拿到任何一段噪声序列先画频谱、看均值、看方差确认它符合理论预期再做后续处理。这个习惯能帮你省下大量排查问题的时间。本文还有配套的精品资源点击获取