簡(jiǎn)介本資源是一套面向數(shù)字信號(hào)處理學(xué)習(xí)者與工程師的MATLAB噪聲生成工具集聚焦白噪聲、布朗噪聲Brown noise與粉紅噪聲Pink noise三類(lèi)典型有色噪聲的建模與仿真適用于系統(tǒng)建模、濾波器設(shè)計(jì)、聽(tīng)覺(jué)感知實(shí)驗(yàn)及生物醫(yī)學(xué)信號(hào)合成等場(chǎng)景。壓縮包共含11個(gè)文件其中6個(gè)為MATLAB函數(shù)腳本.m涵蓋核心噪聲生成器NoiseGenerator.m、時(shí)變AR模型TimeVariantAR.m、彩色噪聲統(tǒng)一接口ColoredNoise.m及多個(gè)測(cè)試用例test*.m5個(gè).mat文件為預(yù)置實(shí)測(cè)樣本數(shù)據(jù)如ECG、EMG、腦電BW等便于快速驗(yàn)證與對(duì)比分析。資源大小2.19MB結(jié)構(gòu)緊湊、即下即用。已有488人學(xué)習(xí)下載提供完整可運(yùn)行代碼、參數(shù)可調(diào)的模塊化設(shè)計(jì)、多噪聲類(lèi)型協(xié)同測(cè)試框架以及C#調(diào)用MATLAB編譯庫(kù)的集成提示顯著降低跨平臺(tái)噪聲仿真開(kāi)發(fā)門(mén)檻。 做了這么多年信號(hào)處理仿真我越來(lái)越確信一個(gè)問(wèn)題很多項(xiàng)目結(jié)果對(duì)不上不是模型寫(xiě)錯(cuò)了而是噪聲加得不對(duì)。白噪聲、Pink噪聲、Brown噪聲這三類(lèi)隨機(jī)信號(hào)是MATLAB仿真里最常用到的噪聲模型但大多數(shù)人對(duì)它們的理解停留在“有個(gè)randn就夠了”的程度。這份關(guān)于“matlab中生成各種噪聲的程序”的zip正好覆蓋了這三類(lèi)噪聲的生成邏輯把這套東西吃透能幫你在系統(tǒng)辨識(shí)、音頻測(cè)試、通信仿真、隨機(jī)游走建模里少走很多彎路。我會(huì)把每一類(lèi)噪聲的生成原理、MATLAB代碼、驗(yàn)證手段以及我在實(shí)際項(xiàng)目中踩過(guò)的坑一次說(shuō)清楚適合剛接觸隨機(jī)信號(hào)處理的學(xué)生也給正在做仿真工程的朋友當(dāng)個(gè)參考。1. 為什么仿真離不開(kāi)噪聲三種噪聲模型的頻譜差異與應(yīng)用場(chǎng)景1.1 白噪聲、Pink噪聲、Brown噪聲到底區(qū)別在哪先別急著看代碼把三個(gè)概念在腦子里立起來(lái)比什么都重要。很多時(shí)候代碼寫(xiě)完發(fā)現(xiàn)頻譜不對(duì)根子就在概念上理解偏了。白噪聲的功率譜密度在整個(gè)頻帶內(nèi)是平坦的也就是每個(gè)頻率點(diǎn)上的能量都相等。用生活里的聲音來(lái)類(lèi)比白噪聲聽(tīng)起來(lái)像雨點(diǎn)密密麻麻打在鐵皮屋頂上嘶嘶啦啦的能量均勻地鋪滿(mǎn)從低頻到高頻的所有頻段。MATLAB里的randn函數(shù)生成的高斯白噪聲就是最典型的白噪聲來(lái)源。Pink噪聲又叫粉紅噪聲它的功率譜密度與頻率成反比也就是說(shuō)頻率越高能量越小每倍頻程大約下降3 dB。它的名字來(lái)源于“粉紅”介于“白”與“紅”之間。聲音上更像遠(yuǎn)處的瀑布或者風(fēng)聲聽(tīng)起來(lái)比白噪聲“悶”一些低頻更飽滿(mǎn)。自然界里很多現(xiàn)象都表現(xiàn)出1/f頻譜的特征所以Pink噪聲在物理、生物、音頻工程里出現(xiàn)頻率極高。Brown噪聲通常也叫布朗噪聲或紅噪聲它的功率譜密度與頻率的平方成反比每倍頻程下降約6 dB。如果你聽(tīng)過(guò)沉悶的雷聲或者大型發(fā)動(dòng)機(jī)的低頻轟鳴那就是Brown噪聲給人的感覺(jué)。它的本質(zhì)是隨機(jī)游走過(guò)程也就是布朗運(yùn)動(dòng)的數(shù)學(xué)模型。注意區(qū)分一個(gè)細(xì)致的地方工程界說(shuō)的Brown噪聲和物理學(xué)里的紅噪聲在定義上略有差異但MATLAB社區(qū)里這類(lèi)zip程序通常指的都是功率按1/f^2衰減的隨機(jī)信號(hào)。用一張表格把三者的核心特征放一起看起來(lái)更直觀(guān)噪聲類(lèi)型功率譜密度幅值譜高頻衰減聲音感受白噪聲平坦平坦無(wú)衰減嘶嘶聲、雨聲Pink噪聲1/f1/sqrt(f)-3 dB/oct瀑布聲、風(fēng)聲Brown噪聲1/f^21/f-6 dB/oct雷聲、低鳴聲1.2 三種噪聲各自的應(yīng)用場(chǎng)景理解了差異再來(lái)看應(yīng)用場(chǎng)景就明白為什么這三類(lèi)噪聲的生成程序總是被打包在一起。白噪聲是系統(tǒng)辨識(shí)的基礎(chǔ)激勵(lì)信號(hào)因?yàn)樗谡麄€(gè)頻帶內(nèi)都能注入能量能夠激發(fā)系統(tǒng)的所有模態(tài)。在通信仿真里高斯白噪聲也是模擬信道加性噪聲的標(biāo)準(zhǔn)模型。它還是生成其他兩類(lèi)噪聲的“原材料”P(pán)ink和Brown噪聲都可以通過(guò)對(duì)白噪聲做變換得到所以白噪聲生成程序是整套方案的底座。Pink噪聲在音頻工程里地位特殊因?yàn)槿硕鷮?duì)頻率的感知大致是對(duì)數(shù)尺度而Pink噪聲每個(gè)倍頻程內(nèi)能量相等用它做測(cè)試信號(hào)最能模擬真實(shí)音樂(lè)和語(yǔ)音的頻譜分布揚(yáng)聲器測(cè)試、房間聲學(xué)測(cè)量、音頻設(shè)備頻響校準(zhǔn)都離不開(kāi)它。此外許多生物醫(yī)學(xué)信號(hào)如腦電、心電也有接近1/f的頻譜特征用Pink噪聲做背景模擬比白噪聲真實(shí)得多。Brown噪聲則大量用于隨機(jī)游走類(lèi)建模比如股價(jià)波動(dòng)、粒子擴(kuò)散、電子器件中的低頻噪聲閃爍噪聲的積分效應(yīng)等場(chǎng)景。在控制系統(tǒng)仿真里它經(jīng)常用來(lái)模擬傳感器漂移——想象一個(gè)陀螺儀的零偏隨著時(shí)間累積慢慢漂走用白噪聲積分得到的Brown噪聲去模擬就很貼切。在MATLAB里生成這三類(lèi)噪聲的代碼都不算長(zhǎng)但每類(lèi)噪聲都有自己的脾氣下面拆開(kāi)來(lái)說(shuō)。2. 白噪聲MATLAB里最基礎(chǔ)也最容易被用錯(cuò)的隨機(jī)源2.1 randn與rand的區(qū)別以及高斯白噪聲的功率含義生成白噪聲第一反應(yīng)通常是randn但很多人不清楚rand和randn到底該選哪個(gè)。rand生成的是均勻分布隨機(jī)數(shù)范圍在(0,1)之間它的功率譜雖然在頻帶內(nèi)也是平的但幅度分布是矩形的。randn生成的是標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)均值為0、方差為1幅度分布是鐘形的。工程上默認(rèn)的高斯白噪聲指的就是randn這一類(lèi)。為什么絕大多數(shù)仿真用高斯白噪聲而不是均勻白噪聲這背后是中心極限定理在起作用。實(shí)際物理系統(tǒng)中的噪聲來(lái)自大量獨(dú)立微小擾動(dòng)的疊加疊加的結(jié)果天然趨向高斯分布。通信系統(tǒng)里的熱噪聲就是典型的高斯白噪聲。所以除非研究的是量化噪聲這類(lèi)特殊場(chǎng)景否則默認(rèn)選randn沒(méi)錯(cuò)。這里有個(gè)特別容易忽略的概念randn生成的標(biāo)準(zhǔn)高斯序列方差是1單位是“幅度平方”也就是功率為1。如果把它直接當(dāng)成信號(hào)喂給系統(tǒng)它的功率譜密度并不是1而是和采樣率有關(guān)。這在量化信噪比的時(shí)候極其關(guān)鍵。% 生成一段高斯白噪聲 fs 44100; % 采樣率單位Hz duration 5; % 時(shí)長(zhǎng)單位s N fs * duration; % 樣本點(diǎn)數(shù) rng(42); % 固定隨機(jī)種子保證可重復(fù) x_white randn(1, N); % 查看基本統(tǒng)計(jì)量 disp([均值: , num2str(mean(x_white))]); disp([方差: , num2str(var(x_white))]);正常情況下均值會(huì)接近0方差接近1。如果你跑出來(lái)的均值和0差得比較多不用慌這是有限樣本的隨機(jī)波動(dòng)樣本數(shù)越大越接近理論值。但如果方差明顯偏離1那就要檢查是不是數(shù)組維度寫(xiě)錯(cuò)了。2.2 按信噪比需求生成指定功率的白噪聲實(shí)際項(xiàng)目里沒(méi)人直接生成一個(gè)“方差為1”的白噪聲就完事更常見(jiàn)的是“在信號(hào)上疊加信噪比為20 dB的噪聲”。這時(shí)候就要把白噪聲的幅度換算好。假設(shè)原始信號(hào)是s信號(hào)功率為P_s要達(dá)到信噪比SNR單位dB那么噪聲功率P_n滿(mǎn)足SNR 10 * log10(P_s / P_n) P_n P_s / 10^(SNR/10)而高斯白噪聲的方差就等于它的功率所以噪聲序列需要乘以sqrt(P_n)也就是噪聲標(biāo)準(zhǔn)差。% 生成指定信噪比的白噪聲并疊加到信號(hào)上 fs 1000; t 0:1/fs:1-1/fs; f0 50; s sin(2*pi*f0*t); % 原始信號(hào)幅值1 P_s mean(s.^2); % 信號(hào)功率 snr_db 20; P_n P_s / (10^(snr_db/10)); noise sqrt(P_n) * randn(1, length(t)); x s noise; snr_actual 10 * log10(P_s / mean(noise.^2)); disp([理論SNR: , num2str(snr_db), dB, 實(shí)測(cè)SNR: , num2str(snr_actual), dB]);實(shí)測(cè)SNR和理論值會(huì)有一點(diǎn)點(diǎn)偏差因?yàn)殡S機(jī)序列的方差本身是隨機(jī)量。如果項(xiàng)目對(duì)SNR誤差要求極嚴(yán)可以用randn生成后再對(duì)該段序列做精確的功率歸一化noise randn(1, length(t)); noise noise - mean(noise); % 先去掉直流分量 noise noise / sqrt(mean(noise.^2)); % 歸一化到單位功率 noise noise * sqrt(P_n); % 再按目標(biāo)功率縮放這樣做能確保這段噪聲的實(shí)際功率嚴(yán)格等于P_n代價(jià)是噪聲序列不再?lài)?yán)格符合獨(dú)立同分布假設(shè)因?yàn)闅w一化引入了微弱的樣本間相關(guān)性。實(shí)際上這個(gè)影響通?？梢院雎?。2.3 白噪聲生成中常見(jiàn)的三個(gè)坑第一個(gè)坑是忘記去掉均值。理論上randn的均值是0但在樣本數(shù)較少時(shí)均值可能會(huì)明顯偏離0這種直流偏置會(huì)直接影響后續(xù)的頻譜分析和系統(tǒng)辨識(shí)結(jié)果。短序列尤其明顯處理前先用mean減去均值是穩(wěn)妥的習(xí)慣。第二個(gè)坑是把方差和幅度搞混。很多人想生成“幅度為0.1的白噪聲”直接寫(xiě)0.1 * randn然后發(fā)現(xiàn)噪聲功率比預(yù)期大了或者小了。白噪聲的幅度是隨機(jī)變量沒(méi)有固定“幅度”只有分布參數(shù)。如果要控制的是功率就必須用方差來(lái)配算。第三個(gè)坑是不設(shè)隨機(jī)種子。仿真項(xiàng)目里如果每次跑出來(lái)結(jié)果都不一樣排查問(wèn)題根本無(wú)從下手。在腳本開(kāi)頭用rng函數(shù)固定種子能保證每次運(yùn)行結(jié)果完全一致這在調(diào)試對(duì)比時(shí)極其重要。等所有邏輯都調(diào)通了再移除種子或者讓種子跟隨時(shí)間變化這才是正確的使用節(jié)奏。3. Pink噪聲的頻域生成法一句話(huà)原理與完整代碼3.1 為什么Pink噪聲不能直接對(duì)白噪聲積分得到說(shuō)到生成Pink噪聲新手最容易想到的思路是“白噪聲做傅里葉變換乘上1/f加權(quán)再反變換回來(lái)”。這個(gè)思路方向是對(duì)的但有一個(gè)物理前提必須先說(shuō)清楚。從頻譜關(guān)系看白噪聲的功率譜是常數(shù)Pink噪聲的功率譜是1/f幅值譜是1/sqrt(f)。所以最簡(jiǎn)單的生成方式就是在頻域?qū)Π自肼暤姆底V做加權(quán)。Brown噪聲的幅值譜是1/f對(duì)白噪聲積分即可得到。而Pink需要的加權(quán)函數(shù)是1/sqrt(f)這個(gè)加權(quán)既不能靠時(shí)域積分一步到位也不能靠簡(jiǎn)單微分得到必須做分?jǐn)?shù)階的處理。這也是為什么Pink噪聲的生成比另外兩類(lèi)要麻煩一些。3.2 頻域加權(quán)法代碼實(shí)現(xiàn)與解釋頻域法是生成Pink噪聲最直觀(guān)的方式基本步驟就四步生成白噪聲序列做FFT在頻域上乘以1/sqrt(f)的幅值加權(quán)再I(mǎi)FFT回到時(shí)域。function x_pink pinknoise_freq(N, fs) % 頻域法生成Pink噪聲 % N: 樣本點(diǎn)數(shù), fs: 采樣率(僅用于頻率計(jì)算不影響生成結(jié)果) % 1. 生成復(fù)白噪聲譜保證相位隨機(jī) white randn(1, N) 1i * randn(1, N); % 2. 構(gòu)造頻率向量 f (0:N-1) * fs / N; % 3. 幅值加權(quán)alpha0.5 對(duì)應(yīng) Pink 噪聲 alpha 0.5; H zeros(1, N); H(2:end) f(2:end).^(-alpha); % 第一個(gè)點(diǎn)(直流)不參與加權(quán) H(1) 0; % 直流分量置零去掉偏移 % 4. 加權(quán)并反變換 X_pink white .* H; x_pink real(ifft(X_pink)); % 5. 歸一化到單位標(biāo)準(zhǔn)差 x_pink x_pink / std(x_pink); end這段代碼的細(xì)節(jié)值得仔細(xì)講。為什么要用復(fù)白噪聲而不是實(shí)白噪聲因?yàn)閺?fù)白噪聲能直接提供隨機(jī)幅度和隨機(jī)相位實(shí)白噪聲FFT之后的頻譜天然是共軛對(duì)稱(chēng)的頻域加權(quán)再反變換也能得到實(shí)信號(hào)但相位結(jié)構(gòu)會(huì)受限制生成效果稍差一些。直流分量為什么要置零因?yàn)?/sqrt(f)在f0處是無(wú)窮大如果不處理反變換出來(lái)的信號(hào)會(huì)帶上一個(gè)大直流偏置導(dǎo)致信號(hào)嚴(yán)重偏離零均值。在工程上這相當(dāng)于給系統(tǒng)注入了一個(gè)巨大的初始偏移后續(xù)處理全亂了。這個(gè)方法的優(yōu)點(diǎn)是好理解、寫(xiě)起來(lái)快但有一個(gè)明顯的隱患因?yàn)镕FT的長(zhǎng)度有限生成的時(shí)域信號(hào)在首尾邊界處會(huì)有不連續(xù)聽(tīng)起來(lái)像周期性重復(fù)的“循環(huán)噪聲”。對(duì)于需要長(zhǎng)時(shí)間連續(xù)噪聲的仿真場(chǎng)景這不是最優(yōu)解。3.3 頻域法的邊界問(wèn)題周期性偽影如何處理我在做音頻測(cè)試的時(shí)候就吃過(guò)這個(gè)虧。用上面這個(gè)函數(shù)生成了一段幾秒鐘的Pink噪聲直接拿來(lái)給揚(yáng)聲器做掃頻?；胤艜r(shí)明顯聽(tīng)到每隔幾秒就有一次“咔噠”聲起初以為是功放問(wèn)題后來(lái)定位到是信號(hào)本身的問(wèn)題。問(wèn)題出在頻域法生成的信號(hào)隱含了周期性假設(shè)IFFT默認(rèn)把這段序列當(dāng)成一個(gè)周期信號(hào)來(lái)處理首尾不連續(xù)就會(huì)在拼接處形成沖擊。處理辦法有兩個(gè)方向。第一個(gè)方向是在頻域生成更長(zhǎng)的序列比如目標(biāo)時(shí)長(zhǎng)的8倍生成后只取中間的一段丟棄首尾的過(guò)渡區(qū)域。這個(gè)方法的代價(jià)是計(jì)算量變大但實(shí)現(xiàn)簡(jiǎn)單。第二個(gè)方向是改用時(shí)域?yàn)V波器法。用一個(gè)白噪聲作為輸入經(jīng)過(guò)一組精心設(shè)計(jì)的濾波器網(wǎng)絡(luò)在時(shí)域上直接輸出Pink噪聲。這種方法的輸出是連續(xù)的沒(méi)有周期性偽影更適合實(shí)時(shí)或長(zhǎng)時(shí)間序列的場(chǎng)景。經(jīng)典的實(shí)現(xiàn)有Paul Kellet的算法通過(guò)一組IIR濾波器來(lái)逼近1/f頻響實(shí)現(xiàn)復(fù)雜一些但效果好。MATLAB的Audio Toolbox里也有現(xiàn)成的函數(shù)不過(guò)如果只是想快速驗(yàn)證算法頻域法配合裁剪已經(jīng)足夠應(yīng)付大部分項(xiàng)目。4. Brown噪聲積分法實(shí)現(xiàn)、直流漂移與工程修正4.1 Brown噪聲和隨機(jī)游走的關(guān)系Brown噪聲的生成是最有意思的因?yàn)樗臄?shù)學(xué)本質(zhì)就是隨機(jī)游走也叫布朗運(yùn)動(dòng)。想象一個(gè)粒子在水面上受到大量隨機(jī)碰撞它的位置每一步都在隨機(jī)變化但位移是累積的——下一步位置是上一步位置加上一個(gè)隨機(jī)擾動(dòng)。這個(gè)過(guò)程的功率譜密度就是1/f^2也就是Brown噪聲。MATLAB里生成Brown噪聲最直接的辦法就是把白噪聲做累計(jì)求和用cumsum函數(shù)一行代碼就能實(shí)現(xiàn)。這個(gè)操作的頻域含義是積分而積分在頻域上相當(dāng)于除以頻率所以白噪聲積分后幅值譜就變成了1/f功率譜變成了1/f^2完全符合Brown噪聲的定義。function x_brown brownnoise_cumsum(N) % 用白噪聲積分生成Brown噪聲 white randn(1, N); x_brown cumsum(white); % 修正直流和漂移 x_brown x_brown - mean(x_brown); % 歸一化到單位標(biāo)準(zhǔn)差 x_brown x_brown / std(x_brown); end4.2 用cumsum生成Brown噪聲的正確姿勢(shì)代碼雖然短但生成之后不能直接拿去用。cumsum出來(lái)的序列有一個(gè)特點(diǎn)它的方差會(huì)隨著樣本索引逐漸增大也就是說(shuō)信號(hào)在時(shí)域上看起來(lái)是一個(gè)“越走越寬”的過(guò)程。直觀(guān)表現(xiàn)是前幾百個(gè)樣本幅度還比較小后面就開(kāi)始大幅波動(dòng)甚至?xí)胶茈x譜的數(shù)值范圍。這在數(shù)學(xué)上完全合理因?yàn)殡S機(jī)游走的方差本來(lái)就隨時(shí)間線(xiàn)性增長(zhǎng)。但在工程仿真里一個(gè)幅度漂到幾百上千的噪聲序列往往不符合物理約束。所以在生成后必須做兩件事去均值和歸一化。去均值是為了消除隨機(jī)游走過(guò)程中的累積偏置。cumsum即使輸入的均值嚴(yán)格為0積累過(guò)程中也會(huì)因?yàn)殡S機(jī)波動(dòng)產(chǎn)生一個(gè)整體偏移不除掉會(huì)導(dǎo)致信號(hào)整體抬高。歸一化同樣關(guān)鍵不除以標(biāo)準(zhǔn)差你看到的Brown噪聲幅度范圍隨機(jī)游走的狀態(tài)而定可能這次是正負(fù)50下次是正負(fù)200完全不可控。除以標(biāo)準(zhǔn)差之后信號(hào)的統(tǒng)計(jì)特性才穩(wěn)定下來(lái)便于后續(xù)和信號(hào)疊加時(shí)做功率計(jì)算。4.3 為什么生成的Brown噪聲必須去趨勢(shì)或高通濾波去均值解決了直流偏移問(wèn)題但解決不了另一個(gè)更隱蔽的問(wèn)題——趨勢(shì)項(xiàng)。隨機(jī)游走序列在長(zhǎng)期演化中可能表現(xiàn)出看上去像線(xiàn)性增長(zhǎng)或下降的“趨勢(shì)”這是隨機(jī)過(guò)程的正常表現(xiàn)并非真的存在確定性趨勢(shì)。但對(duì)于頻譜分析這個(gè)慢變趨勢(shì)意味著超低頻分量被嚴(yán)重放大功率譜在極低頻段會(huì)異常高偏離理論上的1/f^2關(guān)系。我在用自制的Brown噪聲驗(yàn)證一個(gè)仿真系統(tǒng)時(shí)發(fā)現(xiàn)低頻段功率譜密度總是比理論值高出好幾個(gè)dB百思不得其解。后來(lái)用detrend函數(shù)把線(xiàn)性趨勢(shì)去掉再算頻譜結(jié)果就干凈多了。% 生成Brown噪聲并做線(xiàn)性去趨勢(shì) N 44100 * 10; white randn(1, N); x_brown cumsum(white); x_brown x_brown - mean(x_brown); x_brown detrend(x_brown, linear); % 去除線(xiàn)性趨勢(shì) x_brown x_brown / std(x_brown);如果對(duì)低頻段的要求更嚴(yán)格可以用高通濾波器把極低頻率成分直接濾掉。但要注意高通濾波的截止頻率不能設(shè)太高否則會(huì)把Brown噪聲本身的低頻特性削掉生成的東西就不算Brown噪聲了。通常選擇截止頻率在0.5 Hz到幾Hz之間取決于你的仿真場(chǎng)景需要保留多低的頻率成分。5. 生成完別急用如何用pwelch驗(yàn)證噪聲頻譜特性5.1 驗(yàn)證的完整步驟從代碼到頻譜圖生成噪聲之后第一個(gè)動(dòng)作應(yīng)該是畫(huà)頻譜而不是急著疊加到系統(tǒng)里跑仿真。我見(jiàn)過(guò)太多人拿著生成的噪聲用了半天最后發(fā)現(xiàn)頻譜不對(duì)整個(gè)實(shí)驗(yàn)白做。驗(yàn)證用MATLAB自帶的pwelch函數(shù)做功率譜密度估計(jì)一次性看清噪聲的品質(zhì)。% 驗(yàn)證三種噪聲的功率譜密度 fs 44100; N 44100 * 20; % 20秒信號(hào) % 生成三種噪聲 x_white randn(1, N); x_pink pinknoise_freq(N, fs); x_brown brownnoise_cumsum(N); % 使用pwelch估計(jì)功率譜密度 [pxx_w, f] pwelch(x_white, hann(8192), 4096, 8192, fs); [pxx_p, ~] pwelch(x_pink, hann(8192), 4096, 8192, fs); [pxx_b, ~] pwelch(x_brown, hann(8192), 4096, 8192, fs); % 繪制log-log坐標(biāo)下的功率譜 figure; loglog(f, pxx_w, LineWidth, 1.2); hold on; loglog(f, pxx_p, LineWidth, 1.2); loglog(f, pxx_b, LineWidth, 1.2); xlabel(頻率 (Hz)); ylabel(功率譜密度); legend(白噪聲, Pink噪聲, Brown噪聲); grid on;把三條曲線(xiàn)畫(huà)在對(duì)數(shù)-對(duì)數(shù)坐標(biāo)下特性一目了然。白噪聲的譜是水平的Pink噪聲的譜是斜率約-1的直線(xiàn)Brown噪聲的譜是斜率約-2的直線(xiàn)。注意這里說(shuō)“約”因?yàn)閜welch是估計(jì)值譜線(xiàn)會(huì)有波動(dòng)但整體趨勢(shì)應(yīng)該非常清晰。5.2 如何判斷生成對(duì)了斜率、頻段范圍與波動(dòng)范圍判斷生成質(zhì)量有三個(gè)指標(biāo)。第一個(gè)是斜率白噪聲斜率接近0Pink噪聲接近-1Brown噪聲接近-2。這個(gè)在loglog圖里可以直觀(guān)比較也可以用線(xiàn)性擬合來(lái)量化。第二個(gè)是頻段范圍。好的噪聲應(yīng)該在整個(gè)奈奎斯特頻帶內(nèi)都保持相應(yīng)的斜率特性。如果高頻段過(guò)早跌落說(shuō)明信號(hào)被某種平滑操作污染了如果低頻段異常上翹多半沒(méi)處理好直流或趨勢(shì)項(xiàng)。第三個(gè)是譜估計(jì)的波動(dòng)范圍。用pwelch做平均估計(jì)時(shí)譜線(xiàn)會(huì)在理論值周?chē)秳?dòng)抖動(dòng)幅度和所用的窗長(zhǎng)以及平均段數(shù)有關(guān)。如果抖動(dòng)特別劇烈說(shuō)明窗長(zhǎng)太短或者段數(shù)不夠適當(dāng)增加窗長(zhǎng)度可以讓譜線(xiàn)更平滑但會(huì)犧牲低頻分辨率。5.3 均值、方差、直方圖輔助檢查除了頻譜還有三個(gè)統(tǒng)計(jì)量值得隨手測(cè)一下均值應(yīng)該在0附近方差應(yīng)該接近歸一化目標(biāo)直方圖形狀則因噪聲類(lèi)型而異。白噪聲的直方圖應(yīng)該是鐘形的高斯曲線(xiàn)Pink噪聲因?yàn)榻?jīng)過(guò)了頻域整形時(shí)域分布仍然接近高斯但尾部略胖Brown噪聲由于積分效應(yīng)直方圖偏平緩但幅度范圍會(huì)明顯大于同長(zhǎng)度的白噪聲。如果直方圖出現(xiàn)明顯偏斜或雙峰說(shuō)明生成過(guò)程引入了意料之外的確定性結(jié)構(gòu)。這里再?gòu)?qiáng)調(diào)一次統(tǒng)計(jì)檢查和頻譜檢查要一起看只看頻譜容易被譜估計(jì)的參數(shù)選擇誤導(dǎo)只看統(tǒng)計(jì)量又發(fā)現(xiàn)不了頻譜形狀的問(wèn)題。兩個(gè)角度互為補(bǔ)充才能確認(rèn)噪聲確實(shí)是“對(duì)的那種”。6. 工程落地噪聲生成程序里那些文檔不會(huì)告訴你的細(xì)節(jié)6.1 隨機(jī)種子設(shè)置與可復(fù)現(xiàn)性前面提過(guò)rng設(shè)置種子對(duì)可復(fù)現(xiàn)性的意義這里專(zhuān)門(mén)展開(kāi)說(shuō)因?yàn)樗档靡粋€(gè)單獨(dú)的章節(jié)。仿真項(xiàng)目的典型痛苦場(chǎng)景是昨天跑的結(jié)果今天復(fù)現(xiàn)不出來(lái)了不是你改錯(cuò)了代碼而是隨機(jī)種子變了。在調(diào)試階段這是災(zāi)難。我自己的習(xí)慣是在腳本最開(kāi)頭集中設(shè)置一次種子并把種子值作為參量放在配置結(jié)構(gòu)體里。% 配置文件頭部統(tǒng)一管理隨機(jī)種子 cfg struct(); cfg.seed 20240601; rng(cfg.seed);這樣做的好處是哪天發(fā)現(xiàn)某個(gè)結(jié)果有問(wèn)題只要把種子改成當(dāng)時(shí)記錄的數(shù)值就能精確復(fù)現(xiàn)當(dāng)時(shí)的隨機(jī)序列逐步排查。等仿真全部完成需要做批量實(shí)驗(yàn)時(shí)再讓種子隨實(shí)驗(yàn)編號(hào)變化保證各批次之間結(jié)果獨(dú)立。6.2 采樣率、時(shí)長(zhǎng)與FFT點(diǎn)數(shù)的匹配關(guān)系生成噪聲時(shí)采樣率的選擇比你想象中更重要。采樣率決定了噪聲的最高頻率分量也就是奈奎斯特頻率。如果后續(xù)要和真實(shí)系統(tǒng)的采樣率對(duì)齊生成噪聲時(shí)就要直接使用目標(biāo)系統(tǒng)的采樣率避免生成之后再重采樣引入額外失真。時(shí)長(zhǎng)方面頻域法生成噪聲時(shí)FFT點(diǎn)數(shù)就是時(shí)域信號(hào)長(zhǎng)度。如果長(zhǎng)度恰好是2的冪FFT計(jì)算速度會(huì)快不少但現(xiàn)代MATLAB對(duì)非2的冪也能優(yōu)化得很好所以沒(méi)必要為了湊2的冪而強(qiáng)行截短信號(hào)。更重要的是FFT點(diǎn)數(shù)要足夠大才能在低頻段有足夠的頻率分辨率。比如想驗(yàn)證10 Hz以下的頻譜特性FFT點(diǎn)數(shù)至少要有幾百否則低頻段只有幾個(gè)離散頻點(diǎn)斜率根本看不出來(lái)。pwelch函數(shù)里的窗長(zhǎng)和重疊點(diǎn)數(shù)也要配套設(shè)置。窗長(zhǎng)決定了單段譜估計(jì)的頻率分辨率重疊點(diǎn)數(shù)越多譜越平滑但計(jì)算量越大。我的默認(rèn)搭配是窗長(zhǎng)8192、重疊4096、FFT點(diǎn)數(shù)8192對(duì)于大多數(shù)音頻和振動(dòng)信號(hào)都?jí)蛴谩?.3 疊加到信號(hào)上時(shí)的功率匹配噪聲生成的最終目的是疊加到系統(tǒng)里做測(cè)試。這時(shí)功率匹配的細(xì)節(jié)決定信噪比是否準(zhǔn)確。前面在白噪聲部分說(shuō)了計(jì)算公式這里補(bǔ)充一個(gè)工程技巧。如果噪聲和信號(hào)不是同時(shí)生成的比如信號(hào)是讀取的實(shí)測(cè)數(shù)據(jù)噪聲是獨(dú)立生成的那么疊加之前一定要先分別計(jì)算兩者的實(shí)際功率再做縮放。不要依賴(lài)之前的“理論值”因?yàn)閷?shí)測(cè)數(shù)據(jù)的功率計(jì)算出來(lái)往往和理論值有偏差。% 信號(hào)實(shí)測(cè)功率 P_s mean(signal.^2); % 噪聲先歸一化到單位功率 noise noise / sqrt(mean(noise.^2)); % 按目標(biāo)SNR疊加 snr_db 15; noise_scaled noise * sqrt(P_s / (10^(snr_db/10))); mixed signal noise_scaled; % 驗(yàn)證實(shí)際SNR snr_after 10 * log10(mean(signal.^2) / mean(noise_scaled.^2));這樣處理之后就算信號(hào)本身的功率和預(yù)期相差不少最終混疊信號(hào)的信噪比仍然是準(zhǔn)確的。6.4 不同算法生成結(jié)果的差異與選型建議最后把三種噪聲的生成算法做個(gè)橫向?qū)Ρ确奖隳惆错?xiàng)目需求選型。噪聲類(lèi)型推薦生成方式優(yōu)點(diǎn)缺點(diǎn)適用場(chǎng)景白噪聲r(shí)andn 功率縮放最簡(jiǎn)單精確可控?zé)o通用、系統(tǒng)辨識(shí)、通信仿真Pink噪聲頻域加權(quán)法原理清晰、實(shí)現(xiàn)快有周期性偽影快速驗(yàn)證、離線(xiàn)分析Pink噪聲時(shí)域?yàn)V波器法連續(xù)無(wú)偽影、實(shí)時(shí)性好代碼復(fù)雜長(zhǎng)時(shí)間音頻測(cè)試、實(shí)時(shí)系統(tǒng)Brown噪聲cumsum積分法一行代碼、數(shù)學(xué)意義明確有低頻漂移隨機(jī)游走、漂移模擬Brown噪聲頻域加權(quán)法幅頻特性精確邊界效應(yīng)對(duì)頻譜形狀要求高的場(chǎng)景我個(gè)人的建議是仿真驗(yàn)證階段先用最簡(jiǎn)單的方法把原理跑通、頻譜驗(yàn)證做好再根據(jù)邊界條件和實(shí)時(shí)性要求決定要不要換更復(fù)雜的算法。很多項(xiàng)目停留在驗(yàn)證階段就結(jié)束了用頻域法完全足夠。一旦要上實(shí)時(shí)系統(tǒng)或者做長(zhǎng)時(shí)間連續(xù)測(cè)試再投入精力換成時(shí)域?yàn)V波器實(shí)現(xiàn)。另外補(bǔ)充一點(diǎn)就算算法不變不同的隨機(jī)種子生成的噪聲波形也完全不同但統(tǒng)計(jì)特性一致。所以對(duì)比實(shí)驗(yàn)時(shí)盡量用同一個(gè)種子生成的噪聲作為對(duì)照組減少隨機(jī)因素對(duì)結(jié)論的干擾。做噪聲生成這件事表面看只是幾行MATLAB代碼實(shí)際上牽扯到對(duì)隨機(jī)過(guò)程頻譜特性的理解、對(duì)FFT邊界效應(yīng)的警惕、對(duì)功率匹配的精確計(jì)算。把這些細(xì)節(jié)都拿捏住你手里的zip程序就不再是幾段能跑的代碼而是一套可以放心依賴(lài)的仿真基礎(chǔ)設(shè)施。尤其建議你在自己的常用腳本里封裝好這三個(gè)函數(shù)每次項(xiàng)目需要不同特性的隨機(jī)激勵(lì)時(shí)直接調(diào)用把采樣率、時(shí)長(zhǎng)、種子作為輸入?yún)?shù)從此再也不用為“噪聲到底對(duì)不對(duì)”這件事操心。本文還有配套的精品資源點(diǎn)擊獲取