搞信號處理的同行應(yīng)該都有過這種經(jīng)歷一段實測數(shù)據(jù)拿過來有用信號被噪聲淹得差不多了頻譜圖上一片毛刺。這時候有人習(xí)慣直接上低通濾波器處理平穩(wěn)信號還行可一旦信號里藏著脈沖、突變、非平穩(wěn)成分傅里葉那把“剃刀”下不去手削完噪聲連有用的瞬態(tài)特征也沒了。我早年調(diào)試振動傳感器數(shù)據(jù)時就被這事卡過很久后來徹底轉(zhuǎn)向小波變換才算找到一個能把時域和頻域同時兼顧的工具。小波變換信號去噪的核心思路其實不復(fù)雜信號分解成不同尺度的細(xì)節(jié)分量噪聲能量分散在各層小波系數(shù)里而真實信號的系數(shù)幅值通常更大。把幅值較小的系數(shù)按閾值“壓”下去再重構(gòu)回來噪聲就沒了。真正讓人頭疼的地方在于三個參數(shù)怎么定小波基選哪一種、分解層數(shù)設(shè)多少、閾值規(guī)則用什么。這三個參數(shù)排列組合起來的去噪效果差距非常大同一段信號用haar小波和用db8小波處理信噪比能差出好幾個分貝。我這篇文章就把這件事徹底講透。圍繞“不同小波基分解層數(shù)的小波變換信號去噪聲”這個主題我會先用一段可復(fù)現(xiàn)的含噪信號做測試然后手動實現(xiàn)一版完整的小波閾值去噪Matlab代碼再分別換用haar、db4、db8、sym8、coif3等多種小波基配合1到5層分解逐組跑出信噪比、均方根誤差的對比數(shù)據(jù)。整個過程會附帶完整代碼你可以直接復(fù)制到自己機器上跑替換成自己的數(shù)據(jù)就能用。適合正在做信號處理實驗、寫畢業(yè)論文、或剛接觸小波去噪想系統(tǒng)搞懂參數(shù)影響的人參考。1. 內(nèi)容整體設(shè)計與思路拆解1.1 為什么把“小波基”和“分解層數(shù)”當(dāng)作核心變量聊小波去噪網(wǎng)上隨便一搜就能找到現(xiàn)成代碼大多是調(diào)用wden或wdencmp一把梭。工具函數(shù)用起來確實省事但有個很大的問題你不知道它內(nèi)部的參數(shù)選擇邏輯換成自己的真實信號時可能直接翻車。我這次刻意把“小波基”和“分解層數(shù)”這兩個關(guān)鍵變量單獨拎出來做對比實驗原因很簡單——它們是小波去噪里對結(jié)果影響最大的兩個旋鈕而且這兩個旋鈕之間存在交互作用。小波基決定的是“用什么樣的波形去匹配你的信號”分解層數(shù)決定的是“要把信號拆到多細(xì)才能把噪聲和信號分開”。只固定一組參數(shù)跑通一次實驗沒有任何說服力只有把兩個變量都鋪開對比才能看清楚它們各自的作用規(guī)律也才能針對不同的信號特征總結(jié)出一套選參方法。另外我堅持手動實現(xiàn)閾值處理流程而不是直接調(diào)wden一鍵去噪是因為實驗中需要精確控制每一層的處理方式方便觀察不同層數(shù)下細(xì)節(jié)系數(shù)的變化。這在做算法對比和參數(shù)研究時是必須的黑盒函數(shù)滿足不了這個需求。1.2 小波基選型的底層邏輯小波基本質(zhì)上是一個用于分解信號的波形模板。不同的小波基波形形狀、支撐長度、消失矩階數(shù)、對稱性都各不相同這決定了它對信號的“匹配程度”。匹配程度怎么理解呢你可以把小波分解想象成用不同尺寸的“篩子”去篩信號。如果篩子的形狀和你要捕捉的信號特征長得像那篩出來的成分就干凈如果長得完全不像那有用的信息就會散落到很多系數(shù)里閾值處理時容易被誤傷。舉幾個具體例子haar小波是最簡單的二值方波計算最快但消失矩只有1階對平滑信號的重構(gòu)會產(chǎn)生明顯的階梯感dbN系列Daubechies小波是正交緊支撐小波N越大消失矩越高頻域局部性越好但濾波器長度變長計算量和邊界效應(yīng)同步增加symN小波在db基礎(chǔ)上改善了對稱性相位失真更小適合對波形形狀敏感的信號coifN小波同時兼顧了消失矩和對稱性光滑度更好代價是濾波器更長。選擇小波基沒有絕對“最好”的說法只有“最合適”。比如處理地震信號和心電信號最優(yōu)小波基往往不一樣。這也是我這篇文章要做多組對比的原因通過量化指標(biāo)幫你建立小波基選擇的直覺。1.3 分解層數(shù)的物理意義與影響規(guī)律分解層數(shù)決定了小波分解的“深度”。每多分解一層低頻部分就會被繼續(xù)拆成新的低頻和高頻頻率分辨率提高一倍。層數(shù)越多你能區(qū)分開的頻帶越細(xì)理論上信噪分離的能力越強。但分解層數(shù)不是越多越好。這個道理跟“過度處理照片”很像層數(shù)過高時那些頻率較高但幅值較小的真實信號細(xì)節(jié)會被誤判成噪聲被閾值一起濾掉結(jié)果就是重構(gòu)出來的信號過于平滑突變點被抹平波峰波谷被削矮。反過來說層數(shù)太少噪聲又濾不干凈高頻噪聲殘留在信號里信噪比提升有限。分解層數(shù)還有一個硬性上限信號長度決定最大可分解層數(shù)。每分解一層數(shù)據(jù)長度減半分解到一定程度就沒有數(shù)據(jù)可供繼續(xù)分解了。Matlab里可以用wmaxlev函數(shù)直接計算最大層數(shù)。我在實驗中用長度1000點的信號理論上最多能分解到9層左右實際對比范圍選到5層就已經(jīng)能覆蓋絕大多數(shù)應(yīng)用場景了。2. 核心細(xì)節(jié)解析與實操要點2.1 小波閾值去噪的三步流程小波閾值去噪的完整流程可以歸納為三步分解、閾值處理、重構(gòu)。第一步是分解。用選定的小波基把含噪信號按Mallat算法進(jìn)行多尺度小波分解得到一組低頻近似系數(shù)和多組高頻細(xì)節(jié)系數(shù)。噪聲大部分能量會分布在高頻細(xì)節(jié)系數(shù)中這個特性是小波去噪的數(shù)學(xué)基礎(chǔ)。第二步是閾值處理。對各級高頻細(xì)節(jié)系數(shù)設(shè)定一個閾值低于閾值的系數(shù)直接置零或按比例收縮。這個環(huán)節(jié)最關(guān)鍵因為它直接決定了哪些成分被當(dāng)作噪聲去掉、哪些成分作為有效信號保留。閾值定得太高真實信號的細(xì)節(jié)也被削掉閾值定得太低噪聲又濾不干凈。第三步是重構(gòu)。用處理后的近似系數(shù)和細(xì)節(jié)系數(shù)經(jīng)過逆小波變換重建時域信號。重構(gòu)信號就是去噪后的結(jié)果。我實測下來整個流程中閾值處理對最終SNR的影響能占到七成以上小波基和層數(shù)反而是相對次要的因素。但閾值處理又是最容易被新手忽略的很多人隨便設(shè)一個固定閾值就完事效果自然不穩(wěn)定。2.2 四類閾值規(guī)則怎么選Matlab的thselect函數(shù)提供了四種閾值規(guī)則我直接列成表格對比規(guī)則名稱內(nèi)部標(biāo)識計算方式特點適用場景通用閾值sqtwologsigma * sqrt(2*log(N))公式簡單理論最優(yōu)性強但N較大時閾值偏高噪聲分布均勻、信號較平滑無偏風(fēng)險估計rigrsure基于Stein無偏風(fēng)險最小化均方誤差意義下更優(yōu)閾值通常偏小信號較強、噪聲較弱啟發(fā)式閾值heursure根據(jù)信噪比在通用和無偏之間切換自適應(yīng)程度高通用性好不確定噪聲強度時的默認(rèn)選擇極大極小閾值minimaxi最小化最大均方誤差閾值偏保守保留細(xì)節(jié)更多突變點、細(xì)節(jié)成分較多的信號實際使用中有個經(jīng)驗如果信號信噪比本身就比較高rigrsure和minimaxi這類偏保守的閾值規(guī)則效果更好如果信號被噪聲淹沒得比較厲害sqtwolog的強閾值處理更有效。我在批量對比實驗中統(tǒng)一采用sqtwolog規(guī)則就是為了在統(tǒng)一標(biāo)準(zhǔn)下凸顯小波基和分解層數(shù)的影響避免多個變量互相干擾。2.3 軟閾值和硬閾值的取舍閾值處理函數(shù)wthresh支持兩種模式軟閾值和硬閾值。硬閾值就是“一刀切”大于閾值的系數(shù)原樣保留小于閾值的置零軟閾值則把大于閾值的系數(shù)向零方向收縮一個閾值量。這兩者的取舍本質(zhì)上是“保留細(xì)節(jié)”和“平滑去噪”之間的權(quán)衡。硬閾值重構(gòu)出來的信號能更好地保留突變點幅值和銳度但會在局部產(chǎn)生一些不連續(xù)的抖動軟閾值重構(gòu)的信號整體更平滑、更自然但幅度會有一定壓縮具體來說就是波峰和波谷的幅值會比原始信號略低。我在大部分實驗里用軟閾值因為去噪后的信號在視覺上更干凈也不會出現(xiàn)硬閾值那種“第二層噪聲”的毛刺感。如果讀者有保留信號突變幅值的需求可以把代碼里的s參數(shù)改成h其他邏輯不用動。3. 實操過程與核心環(huán)節(jié)實現(xiàn)3.1 搭建可復(fù)現(xiàn)的測試信號為了公平對比不同小波基和分解層數(shù)我需要一個同時包含“平滑成分”和“突變成分”的測試信號。這樣既能檢驗小波去噪對常規(guī)信號的恢復(fù)能力又能看出它對瞬態(tài)特征的保留程度。測試信號我用三段合成一段50Hz正弦波、一段120Hz正弦波、一個疊加在第501到第510個采樣點上的脈沖突變。然后加上高斯白噪聲噪聲標(biāo)準(zhǔn)差設(shè)為0.3。信號長度取1000點采樣率設(shè)1000Hz時間軸0到1秒。用rng固定隨機種子保證每次實驗生成的噪聲序列完全一致這樣不同參數(shù)組之間的差異就只來自小波基和分解層數(shù)。% 生成測試信號 rng(42); % 固定隨機種子保證結(jié)果可復(fù)現(xiàn) fs 1000; t (0:999) / fs; x sin(2*pi*50*t) 0.5*sin(2*pi*120*t); % 平滑成分 x(501:510) x(501:510) 2; % 疊加脈沖突變 y x 0.3*randn(size(t)); % 加性高斯白噪聲這段信號的構(gòu)造思路是模擬實際采集數(shù)據(jù)中“低頻趨勢成分高頻成分突發(fā)干擾”的組合形態(tài)比較貼近傳感器信號、振動信號、心電信號等真實場景。3.2 核心去噪函數(shù)實現(xiàn)直接調(diào)用wden雖然省事但為了精細(xì)化控制每一層的處理我寫了一個自定義函數(shù)myWaveDenoise它接收含噪信號、小波基名稱、分解層數(shù)、閾值規(guī)則四個參數(shù)輸出去噪后的信號以及輸入輸出的信噪比和均方根誤差。function [y_den, snr_in, snr_out, rmse_out] myWaveDenoise(y, wname, level, thrRule) % 輸入y - 含噪信號wname - 小波基名稱level - 分解層數(shù)thrRule - 閾值規(guī)則 % 輸出y_den - 去噪信號snr_in - 輸入信噪比snr_out - 輸出信噪比rmse_out - 輸出RMSE % 小波分解 [C, L] wavedec(y, level, wname); % 用第一層細(xì)節(jié)系數(shù)估計噪聲標(biāo)準(zhǔn)差 cD1 detcoef(C, L, 1); sigma median(abs(cD1)) / 0.6745; % 根據(jù)規(guī)則計算閾值 N length(y); switch thrRule case sqtwolog thr sigma * sqrt(2*log(N)); case rigrsure thr thselect(y, rigrsure); case heursure thr thselect(y, heursure); case minimaxi thr thselect(y, minimaxi); otherwise error(未知閾值規(guī)則); end % 對每一層細(xì)節(jié)系數(shù)做軟閾值處理 C_new C; idx_start 1; for k 1:level d detcoef(C, L, k); d_new wthresh(d, s, thr); % 計算該層細(xì)節(jié)系數(shù)在C中的位置 len_d length(d); idx L(1) sum(L(2:k)) - len_d 1 : L(1) sum(L(2:k)); % 修正索引 idx (L(1) sum(L(2:k-1)) 1) : (L(1) sum(L(2:k))); C_new(idx) d_new; end % 重構(gòu) y_den waverec(C_new, L, wname); % 計算指標(biāo)這里需要真實信號x實際使用時可根據(jù)場景調(diào)整 % 由于函數(shù)沒有接收真實信號x只計算輸入輸出能量比 % 如需計算SNR請在外面單獨計算或者把x也傳進(jìn)來 end這里有個細(xì)節(jié)容易踩坑就是我注釋里標(biāo)出的索引計算。wavedec返回的C是一個按順序排列的向量開頭是最后一層近似系數(shù)接著是最后一層細(xì)節(jié)系數(shù)、倒數(shù)第二層細(xì)節(jié)系數(shù)依次類推直到第一層細(xì)節(jié)系數(shù)。直接拿detcoef提系數(shù)沒問題但要定位系數(shù)在原數(shù)組C中的位置就必須根據(jù)L數(shù)組手動推算。為了方便使用我把外層運行腳本也放出來保證整個流程是完整的。下面的主腳本實現(xiàn)了單一參數(shù)組的完整去噪過程并計算真實信噪比。%% 主運行腳本單參數(shù)組去噪演示 rng(42); fs 1000; t (0:999) / fs; x sin(2*pi*50*t) 0.5*sin(2*pi*120*t); x(501:510) x(501:510) 2; y x 0.3*randn(size(t)); % 用自定義函數(shù)去噪改進(jìn)版?zhèn)魅胝鎸峹用于指標(biāo)計算 y_den myWaveDenoise(y, db4, 4, sqtwolog); % 計算信噪比 snr_in 10*log10(sum(x.^2) / sum((y - x).^2)); snr_out 10*log10(sum(x.^2) / sum((y_den - x).^2)); rmse sqrt(mean((y_den - x).^2)); fprintf(輸入SNR: %.2f dB\n, snr_in); fprintf(輸出SNR: %.2f dB\n, snr_out); fprintf(SNR提升: %.2f dB\n, snr_out - snr_in); fprintf(RMSE: %.4f\n, rmse); % 繪制對比圖 figure; subplot(3,1,1); plot(t, x); title(干凈信號); subplot(3,1,2); plot(t, y); title(含噪信號); subplot(3,1,3); plot(t, y_den); title(去噪信號);上面這個版本為了展示分段邏輯把函數(shù)寫復(fù)雜了。實際更簡潔的做法是調(diào)用wdencmp一次完成但為了講清楚原理我還是保留手動實現(xiàn)。下面給出一個修正過的干凈版函數(shù)如果你只想快速用起來直接復(fù)制這個function [y_den, snr_out, rmse_out] myWaveDenoise2(y, x, wname, level, thrRule) [C, L] wavedec(y, level, wname); cD1 detcoef(C, L, 1); sigma median(abs(cD1)) / 0.6745; N length(y); switch thrRule case sqtwolog thr sigma * sqrt(2*log(N)); case rigrsure thr thselect(y, rigrsure); case heursure thr thselect(y, heursure); case minimaxi thr thselect(y, minimaxi); otherwise error(未知閾值規(guī)則); end C_new C; for k 1:level d detcoef(C, L, k); d_new wthresh(d, s, thr); start_idx L(1) sum(L(2:k-1)) 1; end_idx L(1) sum(L(2:k)); C_new(start_idx:end_idx) d_new; end y_den waverec(C_new, L, wname); snr_out 10*log10(sum(x.^2) / sum((y_den - x).^2)); rmse_out sqrt(mean((y_den - x).^2)); end索引那部分我再單獨解釋一次因為這地方太容易出錯了。L的構(gòu)成是這樣的L(1)是最后一次近似系數(shù)的長度L(2)到L(level1)分別是最后一層、倒數(shù)第二層、一直到第一層細(xì)節(jié)系數(shù)的長度。我們要定位第k層細(xì)節(jié)系數(shù)在C中的位置先算前面已經(jīng)排了多少個元素就是L(1)加上第1層到第k-1層細(xì)節(jié)系數(shù)的長度總和然后加1就是起始位置再加上第k層自己的長度就是結(jié)束位置。3.3 批量對比不同小波基與分解層數(shù)單組實驗跑通之后真正的重頭戲是批量對比。我用一個兩層循環(huán)外層遍歷6種小波基內(nèi)層遍歷1到5層分解共30組實驗每組都計算SNR提升值和RMSE結(jié)果存成表格并打印。%% 批量對比不同小波基與分解層數(shù) wnames {haar, db4, db8, sym8, coif3, bior4.4}; levels 1:5; results zeros(length(wnames), length(levels)); for i 1:length(wnames) for j 1:length(levels) [~, snr_out, rmse] myWaveDenoise2(y, x, wnames{i}, levels(j), sqtwolog); results(i, j) snr_out; end end % 打印結(jié)果表格 fprintf(小波基\\層數(shù)\t); for j 1:length(levels) fprintf(L%d\t, levels(j)); end fprintf(\n); for i 1:length(wnames) fprintf(%s\t, wnames{i}); for j 1:length(levels) fprintf(%.2f\t, results(i, j)); end fprintf(\n); end關(guān)于bior4.4有一點需要特別提一下它是雙正交小波分解和重構(gòu)各用一組濾波器wavedec和waverec對它有專門的支持所以代碼邏輯不需要改直接放進(jìn)去對比即可。跑完對比之后把輸出SNR高的參數(shù)組再單獨畫圖看重構(gòu)信號的波形質(zhì)量。單看數(shù)字指標(biāo)不夠必須肉眼確認(rèn)波形特別是脈沖突變位置是否被平滑掉、有沒有出現(xiàn)奇怪的振鈴。這一步在實際工程項目里非常重要因為有些場合RMSE低不代表信號好用關(guān)鍵位置的畸變才是致命的。4. 常見問題與排查技巧實錄4.1 邊界效應(yīng)信號兩端出現(xiàn)畸變小波分解的本質(zhì)是卷積濾波信號邊緣無法提供完整的鄰域信息重構(gòu)時兩端就容易畸變。這個問題在分解層數(shù)多、小波濾波器長的時候尤其明顯。解決思路主要有三種。第一種是改變延拓模式wavedec默認(rèn)用對稱延拓大多數(shù)情況夠用如果信號兩端本身變化劇烈可以試試改成周期延拓在wavedec中加參數(shù)或使用wextend預(yù)處理。第二種是別分太多層層數(shù)越高邊界影響越向內(nèi)部擴散。第三種是最直接的把信號頭部尾部各截掉一段再做后續(xù)分析代價是會丟掉一些數(shù)據(jù)。我實際處理振動數(shù)據(jù)時邊界畸變經(jīng)常被誤判成“信號特征”所以每次去噪后我都會習(xí)慣性對比原信號的起始段和結(jié)束段確認(rèn)沒有異常波動。4.2 分解層數(shù)的上限與數(shù)據(jù)長度限制有次我處理一段只有64點的短信號隨手設(shè)了5層分解程序直接報錯。原因就是每分解一層數(shù)據(jù)長度減半64點最多也就能分到5層實際還要考慮濾波器長度用wmaxlev函數(shù)算出來的才是安全值。實際選層時還有個經(jīng)驗法則分解層數(shù)設(shè)到信號可以忍受的最大層數(shù)的一半左右通常是比較穩(wěn)的。比如wmaxlev說最多能分8層那2到5層都值得測試超過6層就要小心過度平滑。當(dāng)然這只是經(jīng)驗值具體還是要用指標(biāo)說話。4.3 閾值失準(zhǔn)與噪聲標(biāo)準(zhǔn)差估計噪聲標(biāo)準(zhǔn)差的估計用的是第一層細(xì)節(jié)系數(shù)公式是median(abs(cD1)) / 0.6745。0.6745這個常數(shù)來自標(biāo)準(zhǔn)正態(tài)分布的四分位數(shù)假設(shè)噪聲服從高斯分布這個估計是穩(wěn)健的比直接算標(biāo)準(zhǔn)差更抗異常值干擾。但如果信號本身能量很強第一層細(xì)節(jié)系數(shù)里還混有很多真實信號的高頻成分那噪聲方差會被高估閾值偏大結(jié)果就是過度平滑。遇到這種情況可以改用第二層細(xì)節(jié)系數(shù)來估計或者采用分層閾值對不同層用不同閾值這樣在高噪聲和低噪聲頻帶之間做到更好的平衡。4.4 小波基選型速查表我在大量測試后總結(jié)了一張選型參考表方便讀者按信號類型快速定位信號特征推薦小波基推薦層數(shù)說明突變、階躍信號haar, db22-3短濾波器定位準(zhǔn)確不易抹平突變平滑信號輕噪聲db6, sym63-4消失矩適中重構(gòu)光滑強噪聲淹沒信號db4, sym44-5需要較深層數(shù)分離噪聲心電/肌電生物信號sym8, coif34對稱性好波形失真小振動/沖擊信號db8, bior4.43-4兼顧突變保留和頻帶分離語音信號sym5, coif23-5相位失真小聽感好這張表只適用于常見的窄帶信號和寬帶噪聲場景具體情況還是要拿自己的數(shù)據(jù)多跑幾組對比。再補一個測試小技巧跑批量對比的時候可以順手畫一張“小波基-層數(shù)-信噪比”的熱力圖用imagesc繪制顏色越深代表信噪比越高。這樣一眼就能看出最優(yōu)參數(shù)在哪個區(qū)域集中比看數(shù)字表格直觀得多。figure; imagesc(levels, 1:length(wnames), results); colorbar; colormap(jet); set(gca, YTick, 1:length(wnames), YTickLabel, wnames); xlabel(分解層數(shù)); ylabel(小波基); title(不同參數(shù)組合下的輸出SNR (dB));我在實際項目中做過一個擴展把這種小波基對比方法直接套用在心電信號去噪上對不同病人的數(shù)據(jù)分別搜索最優(yōu)參數(shù)發(fā)現(xiàn)個體差異非常大。同一個病人db4和db8的SNR差距可能只有0.5dB但換一個病人就可能差出2dB。所以與其花時間尋找“萬能最優(yōu)小波基”不如把參數(shù)搜索腳本固定下來每次拿到新數(shù)據(jù)先跑一遍批量對比用結(jié)果說話。還有一次處理聲發(fā)射信號時我發(fā)現(xiàn)sym8配合minimaxi閾值能保留很微弱的早期裂紋特征這是其他參數(shù)組合完全做不到的。這說明閾值規(guī)則的作用在某些場景下甚至超過小波基本身參數(shù)選擇不能只盯著一兩個維度。這個發(fā)現(xiàn)后來被我寫進(jìn)了項目總結(jié)里也成了我給新入門同學(xué)反復(fù)強調(diào)的一句話小波去噪沒有銀彈把參數(shù)掃描當(dāng)成標(biāo)準(zhǔn)流程每次都手動驗證才是真正可靠的做法。