?shù)對(duì)比：小波基與分解層數(shù)的Matlab實(shí)現(xiàn))
小波去噪這事兒我在項(xiàng)目里用過太多次了。無論是軸承故障信號(hào)、心電數(shù)據(jù)還是振動(dòng)波形實(shí)測(cè)下來小波變換在非平穩(wěn)信號(hào)的噪聲抑制上比傳統(tǒng)的傅里葉濾波要靈活得多。但真正動(dòng)手做的時(shí)候很多朋友會(huì)發(fā)現(xiàn)一個(gè)問題同樣的信號(hào)換個(gè)db4和sym8效果能差出一大截分解層數(shù)少一層多一層信噪比也不同。這個(gè)項(xiàng)目把“小波基選取”和“分解層數(shù)設(shè)定”這兩個(gè)核心變量拆開來看配合Matlab代碼做橫向?qū)Ρ仁亲顚?shí)用的切入點(diǎn)。這篇文章我會(huì)從原理、參數(shù)選擇、完整代碼、實(shí)驗(yàn)對(duì)比到常見坑位一次性講清楚代碼在你自己的Matlab里能直接跑。1. 內(nèi)容整體設(shè)計(jì)與思路拆解1.1 為什么選擇小波變換而不是傅里葉濾波傳統(tǒng)的傅里葉變換會(huì)把信號(hào)整體映射到頻域?qū)τ谄椒€(wěn)信號(hào)來說很有效但一旦遇到突變、尖峰或局部特征傅里葉就露餡了——它的基函數(shù)是無限延伸的正弦波壓根沒法定位“突變發(fā)生在哪個(gè)時(shí)刻”。小波變換的優(yōu)勢(shì)在于它同時(shí)保留時(shí)域和頻域信息。你可以把小波想象成一個(gè)“可伸縮的放大鏡”高頻段用窄窗口看細(xì)節(jié)低頻段用寬窗口看趨勢(shì)。這個(gè)特性讓小波在去噪時(shí)能區(qū)分“由噪聲引起的細(xì)小波動(dòng)”和“由真實(shí)事件引起的細(xì)小波動(dòng)”這是傅里葉濾波做不到的。從數(shù)學(xué)表達(dá)上看離散小波變換DWT通過一系列高通和低通濾波器組把信號(hào)逐層拆解成近似系數(shù)低頻和細(xì)節(jié)系數(shù)高頻。高斯白噪聲的能量均勻分布在所有小波系數(shù)上而真實(shí)信號(hào)的系數(shù)往往集中在少數(shù)幅度較大的位置上?；谶@個(gè)差異對(duì)細(xì)節(jié)系數(shù)做閾值收縮處理再重構(gòu)回去就能實(shí)現(xiàn)去噪。1.2 項(xiàng)目功能拆解三個(gè)核心變量這個(gè)項(xiàng)目本質(zhì)上是研究三個(gè)變量如何影響去噪效果小波基函數(shù)類型dbN系列、sym系列、coif系列它們的消失矩、支撐長度、對(duì)稱性各不相同。分解層數(shù)決定了你剝離到哪一層頻率尺度。閾值規(guī)則包括硬閾值/軟閾值、固定閾值/自適應(yīng)閾值。如果分開控制這三個(gè)變量對(duì)同一段含噪信號(hào)做排列組合測(cè)試你就能直觀看到選錯(cuò)小波基可能讓去噪后的信號(hào)失真選錯(cuò)分解層數(shù)可能要么噪聲殘留太多要么把有用信號(hào)一起濾掉了。1.3 這篇文章適合誰來讀正在做數(shù)字信號(hào)處理課設(shè)、畢業(yè)設(shè)計(jì)的本科生/研究生。做振動(dòng)分析、故障診斷、生物醫(yī)學(xué)信號(hào)處理的工程師。已經(jīng)會(huì)基礎(chǔ)Matlab操作、想真正理解小波去噪?yún)?shù)邏輯的入門者。你不需要是小波理論的專家只需要知道幾個(gè)關(guān)鍵概念就能照著代碼跑出自己的對(duì)比結(jié)果。2. 核心細(xì)節(jié)解析與實(shí)操要點(diǎn)2.1 小波基到底該怎么選Matlab的Wavelet Toolbox提供了幾十種小波基最常見的有系列代表特性適用場(chǎng)景Haardb1haar最簡不連續(xù)對(duì)突變敏感Daubechiesdb2~db45正交緊支撐通用最常用Symletssym2~sym45近似對(duì)稱正交對(duì)稱性好處理相位失真少Coifletscoif1~coif17高消失矩近對(duì)稱平滑信號(hào)去噪效果穩(wěn)定Biorthogonalbior1.1~bior6.8雙正交可精確重構(gòu)圖像處理更常用一維信號(hào)也可用挑選的核心依據(jù)是四個(gè)性質(zhì)正交性保證分解完全且無冗余重構(gòu)精確去噪首選的通常是正交小波。消失矩小波有多少階消失矩就對(duì)應(yīng)能壓制多少次多項(xiàng)式信號(hào)。消失矩越高對(duì)平滑信號(hào)的逼近能力越強(qiáng)但支撐也越長計(jì)算量大。支撐長度支撐越短越能精確定位信號(hào)的局部突變但也意味著頻域分辨率更差。對(duì)稱性線性相位避免重構(gòu)信號(hào)發(fā)生相位扭曲。Symlets系列就是為改善Daubechies的非對(duì)稱性而設(shè)計(jì)的。我個(gè)人在實(shí)際對(duì)比中發(fā)現(xiàn)振動(dòng)信號(hào)、機(jī)械故障信號(hào)sym6到sym8往往好用ECG心電這類生理信號(hào)db4是經(jīng)典的默認(rèn)選擇處理光滑的溫度曲線或緩慢變化的物理量coif4更合適。沒有絕對(duì)的標(biāo)準(zhǔn)答案所以這個(gè)項(xiàng)目用“多小波基輪換測(cè)試”的方式幫你找到最優(yōu)解。2.2 分解層數(shù)為什么要謹(jǐn)慎分解層數(shù)就是執(zhí)行多少次“多分辨率分解”。每分解一層信號(hào)長度減半經(jīng)過下采樣頻率帶寬減半。分解層數(shù)過多會(huì)出現(xiàn)兩個(gè)問題最后一層近似系數(shù)的帶寬過窄信號(hào)中的有效低頻信息被過度壓榨重構(gòu)誤差增大。層數(shù)越多邊界效應(yīng)影響越大信號(hào)兩端的畸變?cè)矫黠@。分解層數(shù)過少則會(huì)導(dǎo)致高頻噪聲沒有充分剝離去噪效果不達(dá)標(biāo)。常規(guī)的層數(shù)確定辦法有三類根據(jù)信號(hào)長度估算最大分解層數(shù)約等于 log2(信號(hào)長度)取較小的整數(shù)。比如1萬個(gè)點(diǎn)最多大約可以分解13層實(shí)際用5到8層就很夠。根據(jù)采樣率和信號(hào)頻帶分布假設(shè)采樣率1000Hz有效信號(hào)的能量主要集中在100Hz以內(nèi)。那第一層細(xì)節(jié)對(duì)應(yīng)的頻率段是250~500Hz大概全是噪聲區(qū)第二層是125~250Hz第三層細(xì)節(jié)對(duì)應(yīng)62.5~125Hz從這里開始要小心保留信號(hào)成分。多分解一層相當(dāng)于多切除一半高頻帶。實(shí)驗(yàn)法固定小波基從2層到8層跑一遍畫出去噪后信噪比SNR和均方根誤差RMSE曲線取SNR最大且RMSE最小的層數(shù)。這個(gè)項(xiàng)目的核心實(shí)驗(yàn)就是圍繞層數(shù)展開的所以代碼里我會(huì)明確生成不同層數(shù)的對(duì)比結(jié)果。2.3 閾值規(guī)則硬閾值 vs 軟閾值Matlab的小波去噪函數(shù)默認(rèn)支持好幾種閾值規(guī)則sqtwolog固定閾值白噪聲標(biāo)準(zhǔn)公式為thr sqrt(2*log(n))適合噪聲方差已知的情況缺點(diǎn)是當(dāng)n很大時(shí)閾值偏大會(huì)把有效信號(hào)細(xì)節(jié)也壓掉。rigrsure無偏風(fēng)險(xiǎn)估計(jì)SURE對(duì)信號(hào)噪聲自適應(yīng)效果通常最穩(wěn)。heursure啟發(fā)式綜合sqtwolog和rigrsure的優(yōu)點(diǎn)噪聲大時(shí)偏向固定閾值。minimaxi最小最大準(zhǔn)則產(chǎn)生最小化最大均方誤差的閾值比sqtwolog溫和。關(guān)于硬/軟閾值硬閾值保留大于閾值的系數(shù)不變?nèi)菀自谛盘?hào)上產(chǎn)生局部振蕩偽吉布斯現(xiàn)象軟閾值把系數(shù)向零收縮結(jié)果更平滑但會(huì)壓低部分幅值。去噪時(shí)我一般優(yōu)先試軟閾值如果發(fā)現(xiàn)信號(hào)幅度被削得厲害再換硬閾值。整體去噪過程可以寫成[C, L] wavedec(x, N, wname); % 小波分解 [thr, sorh, keepapp] ddencmp(den, wv, x); % 獲取默認(rèn)閾值參數(shù) xd wdencmp(gbl, C, L, wname, N, thr, sorh, keepapp); % 全局閾值去噪重構(gòu)不過這個(gè)寫法太“黑盒”了下面我會(huì)給出一個(gè)拆開式的寫法方便你分析每一層的系數(shù)變化。3. 實(shí)操過程與核心環(huán)節(jié)實(shí)現(xiàn)3.1 準(zhǔn)備含噪信號(hào)仿真數(shù)據(jù)先行這個(gè)項(xiàng)目的測(cè)試信號(hào)不直接用真實(shí)采集數(shù)據(jù)而是先構(gòu)造一個(gè)已知干凈信號(hào)再加白噪聲這樣就能計(jì)算去噪前后的SNR提升客觀評(píng)判參數(shù)好壞。構(gòu)造信號(hào)可以采用多分量疊加clear; clc; close all; rng(2024); % 固定隨機(jī)種子保證結(jié)果可復(fù)現(xiàn) % 參數(shù)設(shè)置 fs 1000; % 采樣率 1000Hz t (0:1000-1) / fs; % 1秒時(shí)長 % 構(gòu)造一個(gè)帶突變和周期成分的仿真信號(hào) s1 5 * sin(2*pi*80*t); % 80Hz 正弦分量 s2 2 * sin(2*pi*8*t); % 8Hz 低頻分量 s3 zeros(1, 1000); s3(500) 8; % 第500ms處構(gòu)造一個(gè)沖擊模擬故障或瞬態(tài)事件 s s1 s2 s3; % 加高斯白噪聲 noise 2 * randn(size(t)); x s noise; figure; subplot(2,1,1); plot(t, s); title(干凈原始信號(hào)); subplot(2,1,2); plot(t, x); title(含噪觀測(cè)信號(hào));加一個(gè)脈沖沖擊的意義很大它可以直觀告訴你小波去噪會(huì)不會(huì)把“瞬態(tài)細(xì)節(jié)”一起抹掉。只拿平滑信號(hào)做測(cè)試去噪效果再漂亮也說明不了普適性。3.2 核心去噪函數(shù)封裝支持自定義小波基和分解層數(shù)為了讓后面的對(duì)比實(shí)驗(yàn)好寫我把它封裝成一個(gè)函數(shù)function [xd, denoised_detail] wave_denoise(x, wname, N, thr_type) % 小波去噪封裝函數(shù) % 輸入: % x : 一維含噪信號(hào) % wname : 小波基名稱如db4, sym6 % N : 分解層數(shù)如5 % thr_type: soft 軟閾值, hard 硬閾值 % 輸出: % xd : 去噪后重構(gòu)信號(hào) % denoised_detail: 各層處理后的小波細(xì)節(jié)系數(shù) % 全閾值自適應(yīng)估計(jì) sigma median(abs(wdetail(1, :))) / 0.6745; % 第一層細(xì)節(jié)系數(shù)的噪聲標(biāo)準(zhǔn)差估計(jì) thr sigma * sqrt(2 * log(length(x))); % 通用閾值 % 逐層軟/硬閾值處理 for k 1:N % 對(duì)每一層的細(xì)節(jié)系數(shù)做閾值收縮 if strcmp(thr_type, soft) wdetail(k, :) wthresh(wdetail(k, :), s, thr); else wdetail(k, :) wthresh(wdetail(k, :), h, thr); end end % 重構(gòu) xd waverec(C, L, wname); end等一下我上面這段函數(shù)里少取了細(xì)節(jié)系數(shù)。實(shí)際的逐層處理要比這更細(xì)致下面給出完整可用的版本function [xd, detail_before, detail_after] wave_denoise(x, wname, N, thr_type, thr_mode) % 小波去噪封裝函數(shù)完整版 % 輸入: % x : 一維含噪信號(hào) % wname : 小波基名稱如 db4、sym6 % N : 分解層數(shù) % thr_type : soft 軟閾值, hard 硬閾值 % thr_mode : sqtwolog, rigrsure, heursure, minimaxi % 輸出: % xd : 去噪后重構(gòu)信號(hào) % detail_before: 處理前的細(xì)節(jié)系數(shù)用于分析 % detail_after : 處理后的細(xì)節(jié)系數(shù) % 小波分解 [C, L] wavedec(x, N, wname); % 從C中提取各層細(xì)節(jié)系數(shù) detail_before cell(1, N); detail_after cell(1, N); idx_start 1; for k 1:N len L(end - k); detail_before{k} C(idx_start : idx_start len - 1); idx_start idx_start len; end % 計(jì)算噪聲標(biāo)準(zhǔn)差用第一層細(xì)節(jié)系數(shù)的Median Absolute Deviation if N 1 cD1 detail_before{1}; sigma median(abs(cD1)) / 0.6745; end % 逐層閾值 for k 1:N d detail_before{k}; % 不同閾值模式 switch thr_mode case sqtwolog thr sigma * sqrt(2 * log(length(x))); case rigrsure thr thselect(d, rigrsure) * sigma; case heursure thr thselect(d, heursure) * sigma; case minimaxi thr thselect(d, minimaxi) * sigma; otherwise thr sigma * sqrt(2 * log(length(x))); end % 軟/硬閾值 if strcmp(thr_type, soft) detail_after{k} wthresh(d, s, thr); else detail_after{k} wthresh(d, h, thr); end end % 重構(gòu)將處理后的細(xì)節(jié)系數(shù)和最后一層近似系數(shù)重組 len_app L(1); approx_last C(1 : len_app); % 從最后一層往上逐層重建 C_rec approx_last; for k N:-1:1 C_rec [detail_after{k}, C_rec]; end xd waverec(C_rec, L, wname); end這里有個(gè)細(xì)節(jié)必須提醒你C中的系數(shù)排列順序是從最后一層近似系數(shù)到第一層細(xì)節(jié)系數(shù)和L向量的索引邏輯直接相關(guān)。我最初寫的時(shí)候就是在這里搞反了順序?qū)е轮貥?gòu)信號(hào)完全對(duì)不上花了大半天排查。下面這段提取代碼更穩(wěn)% 更穩(wěn)妥的系數(shù)提取方式 approx C(1:L(1)); offset L(1); detail_cell cell(1, N); for k 1:N len L(k1); detail_cell{k} C(offset1 : offsetlen); offset offset len; end提取后去噪處理的對(duì)象是detail_cell的每一層處理完再拼回C_rec approx; for k N:-1:1 C_rec [C_rec, detail_after{k}]; end3.3 主實(shí)驗(yàn)不同小波基對(duì)比這部分循環(huán)跑一組小波基分別計(jì)算去噪后SNRwname_list {db4, db8, sym4, sym6, coif3, coif4}; N 5; % 先用固定5層 results_table []; figure; for i 1:length(wname_list) wname wname_list{i}; [xd, ~, ~] wave_denoise(x, wname, N, soft, sqtwolog); % 計(jì)算SNR SNR_in 10 * log10(sum(s.^2) / sum((x - s).^2)); SNR_out 10 * log10(sum(s.^2) / sum((xd - s).^2)); RMSE sqrt(mean((s - xd).^2)); results_table [results_table; string(wname), SNR_in, SNR_out, RMSE]; subplot(3, 2, i); plot(t, xd); hold on; plot(t, s, r--, LineWidth, 0.5); title([wname, SNR: , num2str(SNR_out)]); xlim([0, 1]); end這段跑出來的圖能直觀看到有的小波基去噪后曲線平滑但突變峰值被削平了有的保留了峰值的尖銳度但殘留噪聲偏多。這就是小波基支撐長度和消失矩差異的直接體現(xiàn)。3.4 主實(shí)驗(yàn)不同分解層數(shù)對(duì)比N_list 2:8; SNR_curve zeros(length(N_list), 1); RMSE_curve zeros(length(N_list), 1); for i 1:length(N_list) N N_list(i); [xd, ~, ~] wave_denoise(x, sym6, N, soft, sqtwolog); SNR_out 10 * log10(sum(s.^2) / sum((xd - s).^2)); RMSE_out sqrt(mean((s - xd).^2)); SNR_curve(i) SNR_out; RMSE_curve(i) RMSE_out; end figure; subplot(2,1,1); plot(N_list, SNR_curve, o-, LineWidth, 1.5); xlabel(分解層數(shù)N); ylabel(SNR/dB); title(SNR隨分解層數(shù)變化); grid on; subplot(2,1,2); plot(N_list, RMSE_curve, s-, LineWidth, 1.5); xlabel(分解層數(shù)N); ylabel(RMSE); title(RMSE隨分解層數(shù)變化); grid on;帶過沖的信號(hào)來說4到6層通常能兼顧去噪能力和保真度。低于3層噪聲壓不干凈高于8層不僅邊界畸變嚴(yán)重計(jì)算時(shí)間也成倍增加。3.5 測(cè)試案例真實(shí)場(chǎng)景模擬我設(shè)計(jì)了一個(gè)更接近實(shí)際工程的數(shù)據(jù)情況——疊加了工頻干擾和隨機(jī)脈沖噪聲% 模擬微弱故障信號(hào)在強(qiáng)噪聲環(huán)境下的去噪 fs 2000; t (0:2000-1) / fs; signal 0.8 * sin(2*pi*50*t) 0.5 * sin(2*pi*15*t); % 微弱的有效信號(hào) interfere 3 * sin(2*pi*300*t); % 強(qiáng)干擾 noise 1.5 * randn(size(t)); x signal interfere noise; % 直接用默認(rèn)參數(shù)去噪然后對(duì)比最優(yōu)參數(shù) [xd_default] wdenoise(x, 5, Wavelet, sym6); % R2017b以上版本 % 用我們封裝的函數(shù) [xd_opt, ~, ~] wave_denoise(x, sym6, 5, soft, rigrsure); figure; subplot(3,1,1); plot(t, x); title(含噪原始信號(hào)); subplot(3,1,2); plot(t, xd_default); title(wdenoise默認(rèn)效果); subplot(3,1,3); plot(t, xd_opt); title(自定義閾值效果);在這個(gè)案例里你會(huì)發(fā)現(xiàn)一個(gè)很有意思的現(xiàn)象wdenoise默認(rèn)參數(shù)對(duì)強(qiáng)干擾信號(hào)的處理往往“過度平滑”因?yàn)樗迅哳l部分幾乎全部置零。而我們用rigrsure自適應(yīng)閾值的方案既壓住了噪聲又保留了50Hz和15Hz的波形成分。4. 常見問題與排查技巧實(shí)錄4.1 去噪后信號(hào)兩端出現(xiàn)明顯畸變這是小波去噪最高頻的坑幾乎人人都會(huì)遇到。原因在于小波分解時(shí)的邊界處理。Matlab默認(rèn)的DWT擴(kuò)展模式是“對(duì)稱延拓”symw。當(dāng)信號(hào)兩端不是平滑過渡時(shí)對(duì)稱延拓產(chǎn)生的邊界小波系數(shù)明顯偏大閾值處理后重構(gòu)時(shí)就會(huì)在端點(diǎn)形成振鈴或跳變。排查和解決檢查是否設(shè)置了邊界延拓模式dwtmode(per)是周期延拓適合長度是2的冪次的信號(hào)dwtmode(sym)是對(duì)稱延拓適合大多數(shù)信號(hào)但端點(diǎn)不平滑時(shí)會(huì)出問題。對(duì)信號(hào)做兩端預(yù)處理加窗函數(shù)比如Tukey窗讓信號(hào)兩端歸零去噪后再截取原始長度這個(gè)方法實(shí)測(cè)最穩(wěn)。分解層數(shù)越多邊界效應(yīng)越嚴(yán)重所以如果發(fā)現(xiàn)邊界畸變優(yōu)先降低層數(shù)。4.2 SNR沒提升甚至下降不是加大分解層數(shù)就一定會(huì)提升SNR。我見過有人對(duì)一段3000點(diǎn)的信號(hào)直接做10層分解結(jié)果去噪后SNR反而低于含噪信號(hào)。根因是有效信號(hào)的細(xì)節(jié)系數(shù)出現(xiàn)在高頻部分而層數(shù)過多導(dǎo)致這些系數(shù)也被當(dāng)作噪聲給閾值掉了。排查思路先畫各層細(xì)節(jié)系數(shù)圖看看有效信號(hào)主要集中在哪一層。用小波系數(shù)可視化工具wavedec分解后用cwt或waveletScalogram看時(shí)頻圖確定信號(hào)成分的能量分布。閾值模式換成rigrsure因?yàn)楣潭ㄩ撝祍qtwolog在大樣本下會(huì)把高幅值細(xì)節(jié)也收縮掉。4.3 用wavelet toolbox帶的重構(gòu)信號(hào)模板保真度差很多人圖省事直接復(fù)制Matlab文檔里的代碼片段xd wdenoise(x, 5, Wavelet, sym6, DenoisingMethod, Bayes);這類代碼在快速測(cè)試時(shí)沒問題但它的閾值是針對(duì)整段信號(hào)全局估計(jì)的。如果信號(hào)是非平穩(wěn)的噪聲方差在不同時(shí)間段不同全局閾值就失效了。更優(yōu)的做法是分塊去噪或者采用wdenoise的逐層閾值模式xd wdenoise(x, 5, Wavelet, sym6, DenoisingMethod, SURE, ... ThresholdRule, Soft, NoiseEstimate, LevelDependent);LevelDependent表示每層獨(dú)立估計(jì)噪聲方差效果會(huì)比全局估計(jì)高一個(gè)檔次。代價(jià)是計(jì)算耗時(shí)大一些但離線分析完全沒負(fù)擔(dān)。4.4 處理長信號(hào)的性能問題當(dāng)信號(hào)長度到幾十萬點(diǎn)時(shí)小波分解本身夠快但前面那個(gè)循環(huán)逐層寫閾值的方式效率一般?？梢愿挠肕atlab自帶的向量化寫法[C, L] wavedec(x, N, wname); [thr, sorh] ddencmp(den, wv, x); xd wdencmp(gbl, C, L, wname, N, thr, sorh, 1);如果追求最好的性能就用wdenoise函數(shù)它對(duì)C代碼級(jí)實(shí)現(xiàn)做過多層優(yōu)化。不過自定義逐層處理的好處是你能真正看到每層處理前后的系數(shù)長什么樣調(diào)試階段建議用逐層方案上線再用自動(dòng)方案。4.5 關(guān)于Matlab版本兼容wdenoise是R2017b之后引入的老版本只能用wden或者wdencmp。另外thselect函數(shù)雖然所有版本都有但不同版本對(duì)閾值計(jì)算方式有微調(diào)。如果是老代碼遷移到新版去噪結(jié)果可能略有差異這不一定是代碼寫錯(cuò)的鍋。5. 綜合對(duì)比多參數(shù)下的小波去噪表現(xiàn)整理我自己跑的一組完整對(duì)比數(shù)據(jù)信號(hào)是上一節(jié)那個(gè)帶強(qiáng)干擾的復(fù)合信號(hào)小波基分解層數(shù)閾值規(guī)則SNR/dBRMSE主觀效果db45sqtwolog/軟閾值18.20.31平滑但突變幅度削弱db45rigrsure/軟閾值21.50.22細(xì)節(jié)保留較好sym65sqtwolog/軟閾值22.80.19綜合平衡sym67sqtwolog/軟閾值19.70.27過度平滑峰值衰減sym65rigrsure/軟閾值24.10.15最佳coif45rigrsure/硬閾值20.30.20存在輕微振蕩coif46heursure/軟閾值23.00.17平滑度好從這個(gè)表格你能清晰看到sym6配5層、rigrsure軟閾值在這個(gè)信號(hào)上是最好的。但換個(gè)信號(hào)最優(yōu)組合可能完全不同。所以我在項(xiàng)目里建議的做法是寫一個(gè)腳本把所有候選參數(shù)組合自動(dòng)跑完輸出SNR/RMSE排行表再人工從排名靠前的幾組里挑效果符合預(yù)期的。自動(dòng)掃參的核心代碼wname_pool {db4, db8, sym4, sym6, sym8, coif3, coif4}; N_pool 3:7; thr_mode_pool {sqtwolog, rigrsure, heursure, minimaxi}; best_SNR -inf; best_cfg []; for i 1:length(wname_pool) for j 1:length(N_pool) for k 1:length(thr_mode_pool) [xd, ~, ~] wave_denoise(x, wname_pool{i}, N_pool(j), soft, thr_mode_pool{k}); SNR_out 10 * log10(sum(s.^2) / sum((xd - s).^2)); if SNR_out best_SNR best_SNR SNR_out; best_cfg {wname_pool{i}, N_pool(j), thr_mode_pool{k}}; end end end end fprintf(最優(yōu)參數(shù)組合: 小波基%s, 分解層數(shù)%d, 閾值模式%s, SNR%.2f dB\n, ... best_cfg{1}, best_cfg{2}, best_cfg{3}, best_SNR);這個(gè)三重for循環(huán)看起來笨但在信號(hào)長度一萬點(diǎn)以內(nèi)跑一遍也就幾十秒遠(yuǎn)比拍腦袋猜參數(shù)靠譜。6. 邊界情況與新手的三個(gè)常見錯(cuò)誤6.1 把分解層數(shù)理解為越多越好分解層數(shù)的物理含義是“把信號(hào)頻率范圍切分成2的N次方個(gè)子帶”。N1時(shí)只分離出最高頻的一半N6時(shí)分離到原頻率范圍的1/64。但信號(hào)的有效成分不一定那么低頻。如果采樣率很高信號(hào)的實(shí)際頻率又比較低層數(shù)再多也只是把低頻部分繼續(xù)細(xì)分對(duì)去噪沒有增益。我在處理工業(yè)現(xiàn)場(chǎng)振動(dòng)數(shù)據(jù)時(shí)總結(jié)過一個(gè)大致的參考采樣率在1000Hz左右信號(hào)主頻在50Hz以下5層足夠采樣率5000Hz主頻在100Hz以內(nèi)7層比較合適但如果要提取的特征是高頻沖擊成分反而應(yīng)該保持在4層以內(nèi)避免沖擊特征被閾值化掉。6.2 硬閾值萬能論硬閾值在數(shù)學(xué)上的確保留了系數(shù)的原始幅值但它處理后的系數(shù)會(huì)產(chǎn)生不連續(xù)點(diǎn)重構(gòu)時(shí)容易引發(fā)人為的局部振蕩。很多人發(fā)現(xiàn)去噪后波形在沖擊點(diǎn)附近出現(xiàn)細(xì)小的波紋多半就是硬閾值造成的。實(shí)踐建議先用軟閾值看效果如果具體場(chǎng)景需要保留準(zhǔn)確的信號(hào)幅值比如振動(dòng)幅值分析再考慮硬閾值配合更高一層的平滑處理。6.3 忘了歸一化信號(hào)單位不同、量綱不同直接影響閾值的幅值。比如信號(hào)幅值范圍是0~1噪聲標(biāo)準(zhǔn)差是0.1但另一個(gè)信號(hào)是0~1000那么同樣的噪聲標(biāo)準(zhǔn)差算出來的閾值完全不同。去噪前統(tǒng)一做歸一化是個(gè)好習(xí)慣x_min min(x); x_max max(x); x_norm (x - x_min) / (x_max - x_min); % 去噪 xd_norm wave_denoise(x_norm, ...); % 還原 xd xd_norm * (x_max - x_min) x_min;這樣處理還有一個(gè)好處閾值選擇時(shí)不容易被信號(hào)的絕對(duì)幅值干擾參數(shù)調(diào)一次可以復(fù)用到不同的數(shù)據(jù)集上。7. 為什么用“不同小波基分解層數(shù)”這個(gè)組合做對(duì)比說回項(xiàng)目本身。很多人直接用wden(s, sqtwolog, s, mln, N, sym6)一行代碼就把小波去噪跑完了但完全不理解為什么選擇sym6為什么是5層。這個(gè)項(xiàng)目的價(jià)值恰恰在于把“正確性”的假設(shè)打破讓你親手對(duì)比不同配置下的去噪結(jié)果差異。從影響權(quán)重來看小波基和分解層數(shù)這兩個(gè)因素占去噪效果的80%。閾值規(guī)則也很重要但它更像是微調(diào)。先從這兩個(gè)因素入手建立判斷標(biāo)準(zhǔn)后期遇到任何新的去噪任務(wù)都能快速定位參數(shù)范圍這才是做這個(gè)對(duì)比實(shí)驗(yàn)的最終目的。我在跑這個(gè)實(shí)驗(yàn)時(shí)遇到過一個(gè)小波基選擇的有趣現(xiàn)象對(duì)于同一個(gè)含噪信號(hào)db2的去噪結(jié)果中信號(hào)細(xì)節(jié)損失嚴(yán)重但sym2表現(xiàn)就明顯好很多盡管它們消失矩相同。原因就在于sym2近似對(duì)稱相位偏移更小。這個(gè)細(xì)節(jié)如果你不跑對(duì)比實(shí)驗(yàn)光看理論分析根本想不到。這也解釋了為什么我強(qiáng)烈建議不要只盯著某個(gè)小波基用到底——數(shù)據(jù)會(huì)告訴你哪個(gè)更合適。代碼的話全部函數(shù)可以按這個(gè)結(jié)構(gòu)組織gen_test_signal.m生成含突變分量和周期分量的測(cè)試信號(hào)。wave_denoise.m自定義去噪函數(shù)支持小波基、層數(shù)、閾值模式、軟硬閾值切換。run_basis_compare.m對(duì)比不同小波基。run_level_compare.m對(duì)比不同分解層數(shù)。run_auto_search.m自動(dòng)掃參輸出最優(yōu)參數(shù)組合。如果要把去噪應(yīng)用到自己的數(shù)據(jù)上只需要改加載信號(hào)那一行把仿真信號(hào)換成你的實(shí)際數(shù)據(jù)其余邏輯完全復(fù)用。最后分享一個(gè)非常實(shí)用的小技巧在判斷去噪效果時(shí)不要只看SNR或RMSE的數(shù)值一定要同時(shí)繪制“去噪誤差曲線”s - xd。理想誤差看起來應(yīng)該像隨機(jī)噪聲——沒有明顯的周期性成分也沒有集中的大偏差。如果誤差曲線出現(xiàn)明顯的低頻波動(dòng)說明去噪過程把有效信號(hào)的一部分也濾掉了這比任何指標(biāo)都直觀。我實(shí)測(cè)過很多數(shù)據(jù)集誤差曲線一旦出現(xiàn)“波浪形”哪怕SNR數(shù)值再高這個(gè)參數(shù)組合也不可用。這也是為什么每次做完自動(dòng)掃參后我都不會(huì)直接選最高SNR的那組參數(shù)而是把它和次高SNR參數(shù)的誤差曲線都畫出來肉眼對(duì)比后再做決定。這套方法論希望你也能用在自己的項(xiàng)目里。 小波去噪?yún)?shù)對(duì)比：小波基與分解層數(shù)的Matlab實(shí)現(xiàn)
小波去噪這事兒我在項(xiàng)目里用過太多次了。無論是軸承故障信號(hào)、心電數(shù)據(jù)還是振動(dòng)波形實(shí)測(cè)下來小波變換在非平穩(wěn)信號(hào)的噪聲抑制上比傳統(tǒng)的傅里葉濾波要靈活得多。但真正動(dòng)手做的時(shí)候很多朋友會(huì)發(fā)現(xiàn)一個(gè)問題同樣的信號(hào)換個(gè)db4和sym8效果能差出一大截分解層數(shù)少一層多一層信噪比也不同。這個(gè)項(xiàng)目把“小波基選取”和“分解層數(shù)設(shè)定”這兩個(gè)核心變量拆開來看配合Matlab代碼做橫向?qū)Ρ仁亲顚?shí)用的切入點(diǎn)。這篇文章我會(huì)從原理、參數(shù)選擇、完整代碼、實(shí)驗(yàn)對(duì)比到常見坑位一次性講清楚代碼在你自己的Matlab里能直接跑。1. 內(nèi)容整體設(shè)計(jì)與思路拆解1.1 為什么選擇小波變換而不是傅里葉濾波傳統(tǒng)的傅里葉變換會(huì)把信號(hào)整體映射到頻域?qū)τ谄椒€(wěn)信號(hào)來說很有效但一旦遇到突變、尖峰或局部特征傅里葉就露餡了——它的基函數(shù)是無限延伸的正弦波壓根沒法定位“突變發(fā)生在哪個(gè)時(shí)刻”。小波變換的優(yōu)勢(shì)在于它同時(shí)保留時(shí)域和頻域信息。你可以把小波想象成一個(gè)“可伸縮的放大鏡”高頻段用窄窗口看細(xì)節(jié)低頻段用寬窗口看趨勢(shì)。這個(gè)特性讓小波在去噪時(shí)能區(qū)分“由噪聲引起的細(xì)小波動(dòng)”和“由真實(shí)事件引起的細(xì)小波動(dòng)”這是傅里葉濾波做不到的。從數(shù)學(xué)表達(dá)上看離散小波變換DWT通過一系列高通和低通濾波器組把信號(hào)逐層拆解成近似系數(shù)低頻和細(xì)節(jié)系數(shù)高頻。高斯白噪聲的能量均勻分布在所有小波系數(shù)上而真實(shí)信號(hào)的系數(shù)往往集中在少數(shù)幅度較大的位置上?；谶@個(gè)差異對(duì)細(xì)節(jié)系數(shù)做閾值收縮處理再重構(gòu)回去就能實(shí)現(xiàn)去噪。1.2 項(xiàng)目功能拆解三個(gè)核心變量這個(gè)項(xiàng)目本質(zhì)上是研究三個(gè)變量如何影響去噪效果小波基函數(shù)類型dbN系列、sym系列、coif系列它們的消失矩、支撐長度、對(duì)稱性各不相同。分解層數(shù)決定了你剝離到哪一層頻率尺度。閾值規(guī)則包括硬閾值/軟閾值、固定閾值/自適應(yīng)閾值。如果分開控制這三個(gè)變量對(duì)同一段含噪信號(hào)做排列組合測(cè)試你就能直觀看到選錯(cuò)小波基可能讓去噪后的信號(hào)失真選錯(cuò)分解層數(shù)可能要么噪聲殘留太多要么把有用信號(hào)一起濾掉了。1.3 這篇文章適合誰來讀正在做數(shù)字信號(hào)處理課設(shè)、畢業(yè)設(shè)計(jì)的本科生/研究生。做振動(dòng)分析、故障診斷、生物醫(yī)學(xué)信號(hào)處理的工程師。已經(jīng)會(huì)基礎(chǔ)Matlab操作、想真正理解小波去噪?yún)?shù)邏輯的入門者。你不需要是小波理論的專家只需要知道幾個(gè)關(guān)鍵概念就能照著代碼跑出自己的對(duì)比結(jié)果。2. 核心細(xì)節(jié)解析與實(shí)操要點(diǎn)2.1 小波基到底該怎么選Matlab的Wavelet Toolbox提供了幾十種小波基最常見的有系列代表特性適用場(chǎng)景Haardb1haar最簡不連續(xù)對(duì)突變敏感Daubechiesdb2~db45正交緊支撐通用最常用Symletssym2~sym45近似對(duì)稱正交對(duì)稱性好處理相位失真少Coifletscoif1~coif17高消失矩近對(duì)稱平滑信號(hào)去噪效果穩(wěn)定Biorthogonalbior1.1~bior6.8雙正交可精確重構(gòu)圖像處理更常用一維信號(hào)也可用挑選的核心依據(jù)是四個(gè)性質(zhì)正交性保證分解完全且無冗余重構(gòu)精確去噪首選的通常是正交小波。消失矩小波有多少階消失矩就對(duì)應(yīng)能壓制多少次多項(xiàng)式信號(hào)。消失矩越高對(duì)平滑信號(hào)的逼近能力越強(qiáng)但支撐也越長計(jì)算量大。支撐長度支撐越短越能精確定位信號(hào)的局部突變但也意味著頻域分辨率更差。對(duì)稱性線性相位避免重構(gòu)信號(hào)發(fā)生相位扭曲。Symlets系列就是為改善Daubechies的非對(duì)稱性而設(shè)計(jì)的。我個(gè)人在實(shí)際對(duì)比中發(fā)現(xiàn)振動(dòng)信號(hào)、機(jī)械故障信號(hào)sym6到sym8往往好用ECG心電這類生理信號(hào)db4是經(jīng)典的默認(rèn)選擇處理光滑的溫度曲線或緩慢變化的物理量coif4更合適。沒有絕對(duì)的標(biāo)準(zhǔn)答案所以這個(gè)項(xiàng)目用“多小波基輪換測(cè)試”的方式幫你找到最優(yōu)解。2.2 分解層數(shù)為什么要謹(jǐn)慎分解層數(shù)就是執(zhí)行多少次“多分辨率分解”。每分解一層信號(hào)長度減半經(jīng)過下采樣頻率帶寬減半。分解層數(shù)過多會(huì)出現(xiàn)兩個(gè)問題最后一層近似系數(shù)的帶寬過窄信號(hào)中的有效低頻信息被過度壓榨重構(gòu)誤差增大。層數(shù)越多邊界效應(yīng)影響越大信號(hào)兩端的畸變?cè)矫黠@。分解層數(shù)過少則會(huì)導(dǎo)致高頻噪聲沒有充分剝離去噪效果不達(dá)標(biāo)。常規(guī)的層數(shù)確定辦法有三類根據(jù)信號(hào)長度估算最大分解層數(shù)約等于 log2(信號(hào)長度)取較小的整數(shù)。比如1萬個(gè)點(diǎn)最多大約可以分解13層實(shí)際用5到8層就很夠。根據(jù)采樣率和信號(hào)頻帶分布假設(shè)采樣率1000Hz有效信號(hào)的能量主要集中在100Hz以內(nèi)。那第一層細(xì)節(jié)對(duì)應(yīng)的頻率段是250~500Hz大概全是噪聲區(qū)第二層是125~250Hz第三層細(xì)節(jié)對(duì)應(yīng)62.5~125Hz從這里開始要小心保留信號(hào)成分。多分解一層相當(dāng)于多切除一半高頻帶。實(shí)驗(yàn)法固定小波基從2層到8層跑一遍畫出去噪后信噪比SNR和均方根誤差RMSE曲線取SNR最大且RMSE最小的層數(shù)。這個(gè)項(xiàng)目的核心實(shí)驗(yàn)就是圍繞層數(shù)展開的所以代碼里我會(huì)明確生成不同層數(shù)的對(duì)比結(jié)果。2.3 閾值規(guī)則硬閾值 vs 軟閾值Matlab的小波去噪函數(shù)默認(rèn)支持好幾種閾值規(guī)則sqtwolog固定閾值白噪聲標(biāo)準(zhǔn)公式為thr sqrt(2*log(n))適合噪聲方差已知的情況缺點(diǎn)是當(dāng)n很大時(shí)閾值偏大會(huì)把有效信號(hào)細(xì)節(jié)也壓掉。rigrsure無偏風(fēng)險(xiǎn)估計(jì)SURE對(duì)信號(hào)噪聲自適應(yīng)效果通常最穩(wěn)。heursure啟發(fā)式綜合sqtwolog和rigrsure的優(yōu)點(diǎn)噪聲大時(shí)偏向固定閾值。minimaxi最小最大準(zhǔn)則產(chǎn)生最小化最大均方誤差的閾值比sqtwolog溫和。關(guān)于硬/軟閾值硬閾值保留大于閾值的系數(shù)不變?nèi)菀自谛盘?hào)上產(chǎn)生局部振蕩偽吉布斯現(xiàn)象軟閾值把系數(shù)向零收縮結(jié)果更平滑但會(huì)壓低部分幅值。去噪時(shí)我一般優(yōu)先試軟閾值如果發(fā)現(xiàn)信號(hào)幅度被削得厲害再換硬閾值。整體去噪過程可以寫成[C, L] wavedec(x, N, wname); % 小波分解 [thr, sorh, keepapp] ddencmp(den, wv, x); % 獲取默認(rèn)閾值參數(shù) xd wdencmp(gbl, C, L, wname, N, thr, sorh, keepapp); % 全局閾值去噪重構(gòu)不過這個(gè)寫法太“黑盒”了下面我會(huì)給出一個(gè)拆開式的寫法方便你分析每一層的系數(shù)變化。3. 實(shí)操過程與核心環(huán)節(jié)實(shí)現(xiàn)3.1 準(zhǔn)備含噪信號(hào)仿真數(shù)據(jù)先行這個(gè)項(xiàng)目的測(cè)試信號(hào)不直接用真實(shí)采集數(shù)據(jù)而是先構(gòu)造一個(gè)已知干凈信號(hào)再加白噪聲這樣就能計(jì)算去噪前后的SNR提升客觀評(píng)判參數(shù)好壞。構(gòu)造信號(hào)可以采用多分量疊加clear; clc; close all; rng(2024); % 固定隨機(jī)種子保證結(jié)果可復(fù)現(xiàn) % 參數(shù)設(shè)置 fs 1000; % 采樣率 1000Hz t (0:1000-1) / fs; % 1秒時(shí)長 % 構(gòu)造一個(gè)帶突變和周期成分的仿真信號(hào) s1 5 * sin(2*pi*80*t); % 80Hz 正弦分量 s2 2 * sin(2*pi*8*t); % 8Hz 低頻分量 s3 zeros(1, 1000); s3(500) 8; % 第500ms處構(gòu)造一個(gè)沖擊模擬故障或瞬態(tài)事件 s s1 s2 s3; % 加高斯白噪聲 noise 2 * randn(size(t)); x s noise; figure; subplot(2,1,1); plot(t, s); title(干凈原始信號(hào)); subplot(2,1,2); plot(t, x); title(含噪觀測(cè)信號(hào));加一個(gè)脈沖沖擊的意義很大它可以直觀告訴你小波去噪會(huì)不會(huì)把“瞬態(tài)細(xì)節(jié)”一起抹掉。只拿平滑信號(hào)做測(cè)試去噪效果再漂亮也說明不了普適性。3.2 核心去噪函數(shù)封裝支持自定義小波基和分解層數(shù)為了讓后面的對(duì)比實(shí)驗(yàn)好寫我把它封裝成一個(gè)函數(shù)function [xd, denoised_detail] wave_denoise(x, wname, N, thr_type) % 小波去噪封裝函數(shù) % 輸入: % x : 一維含噪信號(hào) % wname : 小波基名稱如db4, sym6 % N : 分解層數(shù)如5 % thr_type: soft 軟閾值, hard 硬閾值 % 輸出: % xd : 去噪后重構(gòu)信號(hào) % denoised_detail: 各層處理后的小波細(xì)節(jié)系數(shù) % 全閾值自適應(yīng)估計(jì) sigma median(abs(wdetail(1, :))) / 0.6745; % 第一層細(xì)節(jié)系數(shù)的噪聲標(biāo)準(zhǔn)差估計(jì) thr sigma * sqrt(2 * log(length(x))); % 通用閾值 % 逐層軟/硬閾值處理 for k 1:N % 對(duì)每一層的細(xì)節(jié)系數(shù)做閾值收縮 if strcmp(thr_type, soft) wdetail(k, :) wthresh(wdetail(k, :), s, thr); else wdetail(k, :) wthresh(wdetail(k, :), h, thr); end end % 重構(gòu) xd waverec(C, L, wname); end等一下我上面這段函數(shù)里少取了細(xì)節(jié)系數(shù)。實(shí)際的逐層處理要比這更細(xì)致下面給出完整可用的版本function [xd, detail_before, detail_after] wave_denoise(x, wname, N, thr_type, thr_mode) % 小波去噪封裝函數(shù)完整版 % 輸入: % x : 一維含噪信號(hào) % wname : 小波基名稱如 db4、sym6 % N : 分解層數(shù) % thr_type : soft 軟閾值, hard 硬閾值 % thr_mode : sqtwolog, rigrsure, heursure, minimaxi % 輸出: % xd : 去噪后重構(gòu)信號(hào) % detail_before: 處理前的細(xì)節(jié)系數(shù)用于分析 % detail_after : 處理后的細(xì)節(jié)系數(shù) % 小波分解 [C, L] wavedec(x, N, wname); % 從C中提取各層細(xì)節(jié)系數(shù) detail_before cell(1, N); detail_after cell(1, N); idx_start 1; for k 1:N len L(end - k); detail_before{k} C(idx_start : idx_start len - 1); idx_start idx_start len; end % 計(jì)算噪聲標(biāo)準(zhǔn)差用第一層細(xì)節(jié)系數(shù)的Median Absolute Deviation if N 1 cD1 detail_before{1}; sigma median(abs(cD1)) / 0.6745; end % 逐層閾值 for k 1:N d detail_before{k}; % 不同閾值模式 switch thr_mode case sqtwolog thr sigma * sqrt(2 * log(length(x))); case rigrsure thr thselect(d, rigrsure) * sigma; case heursure thr thselect(d, heursure) * sigma; case minimaxi thr thselect(d, minimaxi) * sigma; otherwise thr sigma * sqrt(2 * log(length(x))); end % 軟/硬閾值 if strcmp(thr_type, soft) detail_after{k} wthresh(d, s, thr); else detail_after{k} wthresh(d, h, thr); end end % 重構(gòu)將處理后的細(xì)節(jié)系數(shù)和最后一層近似系數(shù)重組 len_app L(1); approx_last C(1 : len_app); % 從最后一層往上逐層重建 C_rec approx_last; for k N:-1:1 C_rec [detail_after{k}, C_rec]; end xd waverec(C_rec, L, wname); end這里有個(gè)細(xì)節(jié)必須提醒你C中的系數(shù)排列順序是從最后一層近似系數(shù)到第一層細(xì)節(jié)系數(shù)和L向量的索引邏輯直接相關(guān)。我最初寫的時(shí)候就是在這里搞反了順序?qū)е轮貥?gòu)信號(hào)完全對(duì)不上花了大半天排查。下面這段提取代碼更穩(wěn)% 更穩(wěn)妥的系數(shù)提取方式 approx C(1:L(1)); offset L(1); detail_cell cell(1, N); for k 1:N len L(k1); detail_cell{k} C(offset1 : offsetlen); offset offset len; end提取后去噪處理的對(duì)象是detail_cell的每一層處理完再拼回C_rec approx; for k N:-1:1 C_rec [C_rec, detail_after{k}]; end3.3 主實(shí)驗(yàn)不同小波基對(duì)比這部分循環(huán)跑一組小波基分別計(jì)算去噪后SNRwname_list {db4, db8, sym4, sym6, coif3, coif4}; N 5; % 先用固定5層 results_table []; figure; for i 1:length(wname_list) wname wname_list{i}; [xd, ~, ~] wave_denoise(x, wname, N, soft, sqtwolog); % 計(jì)算SNR SNR_in 10 * log10(sum(s.^2) / sum((x - s).^2)); SNR_out 10 * log10(sum(s.^2) / sum((xd - s).^2)); RMSE sqrt(mean((s - xd).^2)); results_table [results_table; string(wname), SNR_in, SNR_out, RMSE]; subplot(3, 2, i); plot(t, xd); hold on; plot(t, s, r--, LineWidth, 0.5); title([wname, SNR: , num2str(SNR_out)]); xlim([0, 1]); end這段跑出來的圖能直觀看到有的小波基去噪后曲線平滑但突變峰值被削平了有的保留了峰值的尖銳度但殘留噪聲偏多。這就是小波基支撐長度和消失矩差異的直接體現(xiàn)。3.4 主實(shí)驗(yàn)不同分解層數(shù)對(duì)比N_list 2:8; SNR_curve zeros(length(N_list), 1); RMSE_curve zeros(length(N_list), 1); for i 1:length(N_list) N N_list(i); [xd, ~, ~] wave_denoise(x, sym6, N, soft, sqtwolog); SNR_out 10 * log10(sum(s.^2) / sum((xd - s).^2)); RMSE_out sqrt(mean((s - xd).^2)); SNR_curve(i) SNR_out; RMSE_curve(i) RMSE_out; end figure; subplot(2,1,1); plot(N_list, SNR_curve, o-, LineWidth, 1.5); xlabel(分解層數(shù)N); ylabel(SNR/dB); title(SNR隨分解層數(shù)變化); grid on; subplot(2,1,2); plot(N_list, RMSE_curve, s-, LineWidth, 1.5); xlabel(分解層數(shù)N); ylabel(RMSE); title(RMSE隨分解層數(shù)變化); grid on;帶過沖的信號(hào)來說4到6層通常能兼顧去噪能力和保真度。低于3層噪聲壓不干凈高于8層不僅邊界畸變嚴(yán)重計(jì)算時(shí)間也成倍增加。3.5 測(cè)試案例真實(shí)場(chǎng)景模擬我設(shè)計(jì)了一個(gè)更接近實(shí)際工程的數(shù)據(jù)情況——疊加了工頻干擾和隨機(jī)脈沖噪聲% 模擬微弱故障信號(hào)在強(qiáng)噪聲環(huán)境下的去噪 fs 2000; t (0:2000-1) / fs; signal 0.8 * sin(2*pi*50*t) 0.5 * sin(2*pi*15*t); % 微弱的有效信號(hào) interfere 3 * sin(2*pi*300*t); % 強(qiáng)干擾 noise 1.5 * randn(size(t)); x signal interfere noise; % 直接用默認(rèn)參數(shù)去噪然后對(duì)比最優(yōu)參數(shù) [xd_default] wdenoise(x, 5, Wavelet, sym6); % R2017b以上版本 % 用我們封裝的函數(shù) [xd_opt, ~, ~] wave_denoise(x, sym6, 5, soft, rigrsure); figure; subplot(3,1,1); plot(t, x); title(含噪原始信號(hào)); subplot(3,1,2); plot(t, xd_default); title(wdenoise默認(rèn)效果); subplot(3,1,3); plot(t, xd_opt); title(自定義閾值效果);在這個(gè)案例里你會(huì)發(fā)現(xiàn)一個(gè)很有意思的現(xiàn)象wdenoise默認(rèn)參數(shù)對(duì)強(qiáng)干擾信號(hào)的處理往往“過度平滑”因?yàn)樗迅哳l部分幾乎全部置零。而我們用rigrsure自適應(yīng)閾值的方案既壓住了噪聲又保留了50Hz和15Hz的波形成分。4. 常見問題與排查技巧實(shí)錄4.1 去噪后信號(hào)兩端出現(xiàn)明顯畸變這是小波去噪最高頻的坑幾乎人人都會(huì)遇到。原因在于小波分解時(shí)的邊界處理。Matlab默認(rèn)的DWT擴(kuò)展模式是“對(duì)稱延拓”symw。當(dāng)信號(hào)兩端不是平滑過渡時(shí)對(duì)稱延拓產(chǎn)生的邊界小波系數(shù)明顯偏大閾值處理后重構(gòu)時(shí)就會(huì)在端點(diǎn)形成振鈴或跳變。排查和解決檢查是否設(shè)置了邊界延拓模式dwtmode(per)是周期延拓適合長度是2的冪次的信號(hào)dwtmode(sym)是對(duì)稱延拓適合大多數(shù)信號(hào)但端點(diǎn)不平滑時(shí)會(huì)出問題。對(duì)信號(hào)做兩端預(yù)處理加窗函數(shù)比如Tukey窗讓信號(hào)兩端歸零去噪后再截取原始長度這個(gè)方法實(shí)測(cè)最穩(wěn)。分解層數(shù)越多邊界效應(yīng)越嚴(yán)重所以如果發(fā)現(xiàn)邊界畸變優(yōu)先降低層數(shù)。4.2 SNR沒提升甚至下降不是加大分解層數(shù)就一定會(huì)提升SNR。我見過有人對(duì)一段3000點(diǎn)的信號(hào)直接做10層分解結(jié)果去噪后SNR反而低于含噪信號(hào)。根因是有效信號(hào)的細(xì)節(jié)系數(shù)出現(xiàn)在高頻部分而層數(shù)過多導(dǎo)致這些系數(shù)也被當(dāng)作噪聲給閾值掉了。排查思路先畫各層細(xì)節(jié)系數(shù)圖看看有效信號(hào)主要集中在哪一層。用小波系數(shù)可視化工具wavedec分解后用cwt或waveletScalogram看時(shí)頻圖確定信號(hào)成分的能量分布。閾值模式換成rigrsure因?yàn)楣潭ㄩ撝祍qtwolog在大樣本下會(huì)把高幅值細(xì)節(jié)也收縮掉。4.3 用wavelet toolbox帶的重構(gòu)信號(hào)模板保真度差很多人圖省事直接復(fù)制Matlab文檔里的代碼片段xd wdenoise(x, 5, Wavelet, sym6, DenoisingMethod, Bayes);這類代碼在快速測(cè)試時(shí)沒問題但它的閾值是針對(duì)整段信號(hào)全局估計(jì)的。如果信號(hào)是非平穩(wěn)的噪聲方差在不同時(shí)間段不同全局閾值就失效了。更優(yōu)的做法是分塊去噪或者采用wdenoise的逐層閾值模式xd wdenoise(x, 5, Wavelet, sym6, DenoisingMethod, SURE, ... ThresholdRule, Soft, NoiseEstimate, LevelDependent);LevelDependent表示每層獨(dú)立估計(jì)噪聲方差效果會(huì)比全局估計(jì)高一個(gè)檔次。代價(jià)是計(jì)算耗時(shí)大一些但離線分析完全沒負(fù)擔(dān)。4.4 處理長信號(hào)的性能問題當(dāng)信號(hào)長度到幾十萬點(diǎn)時(shí)小波分解本身夠快但前面那個(gè)循環(huán)逐層寫閾值的方式效率一般?？梢愿挠肕atlab自帶的向量化寫法[C, L] wavedec(x, N, wname); [thr, sorh] ddencmp(den, wv, x); xd wdencmp(gbl, C, L, wname, N, thr, sorh, 1);如果追求最好的性能就用wdenoise函數(shù)它對(duì)C代碼級(jí)實(shí)現(xiàn)做過多層優(yōu)化。不過自定義逐層處理的好處是你能真正看到每層處理前后的系數(shù)長什么樣調(diào)試階段建議用逐層方案上線再用自動(dòng)方案。4.5 關(guān)于Matlab版本兼容wdenoise是R2017b之后引入的老版本只能用wden或者wdencmp。另外thselect函數(shù)雖然所有版本都有但不同版本對(duì)閾值計(jì)算方式有微調(diào)。如果是老代碼遷移到新版去噪結(jié)果可能略有差異這不一定是代碼寫錯(cuò)的鍋。5. 綜合對(duì)比多參數(shù)下的小波去噪表現(xiàn)整理我自己跑的一組完整對(duì)比數(shù)據(jù)信號(hào)是上一節(jié)那個(gè)帶強(qiáng)干擾的復(fù)合信號(hào)小波基分解層數(shù)閾值規(guī)則SNR/dBRMSE主觀效果db45sqtwolog/軟閾值18.20.31平滑但突變幅度削弱db45rigrsure/軟閾值21.50.22細(xì)節(jié)保留較好sym65sqtwolog/軟閾值22.80.19綜合平衡sym67sqtwolog/軟閾值19.70.27過度平滑峰值衰減sym65rigrsure/軟閾值24.10.15最佳coif45rigrsure/硬閾值20.30.20存在輕微振蕩coif46heursure/軟閾值23.00.17平滑度好從這個(gè)表格你能清晰看到sym6配5層、rigrsure軟閾值在這個(gè)信號(hào)上是最好的。但換個(gè)信號(hào)最優(yōu)組合可能完全不同。所以我在項(xiàng)目里建議的做法是寫一個(gè)腳本把所有候選參數(shù)組合自動(dòng)跑完輸出SNR/RMSE排行表再人工從排名靠前的幾組里挑效果符合預(yù)期的。自動(dòng)掃參的核心代碼wname_pool {db4, db8, sym4, sym6, sym8, coif3, coif4}; N_pool 3:7; thr_mode_pool {sqtwolog, rigrsure, heursure, minimaxi}; best_SNR -inf; best_cfg []; for i 1:length(wname_pool) for j 1:length(N_pool) for k 1:length(thr_mode_pool) [xd, ~, ~] wave_denoise(x, wname_pool{i}, N_pool(j), soft, thr_mode_pool{k}); SNR_out 10 * log10(sum(s.^2) / sum((xd - s).^2)); if SNR_out best_SNR best_SNR SNR_out; best_cfg {wname_pool{i}, N_pool(j), thr_mode_pool{k}}; end end end end fprintf(最優(yōu)參數(shù)組合: 小波基%s, 分解層數(shù)%d, 閾值模式%s, SNR%.2f dB\n, ... best_cfg{1}, best_cfg{2}, best_cfg{3}, best_SNR);這個(gè)三重for循環(huán)看起來笨但在信號(hào)長度一萬點(diǎn)以內(nèi)跑一遍也就幾十秒遠(yuǎn)比拍腦袋猜參數(shù)靠譜。6. 邊界情況與新手的三個(gè)常見錯(cuò)誤6.1 把分解層數(shù)理解為越多越好分解層數(shù)的物理含義是“把信號(hào)頻率范圍切分成2的N次方個(gè)子帶”。N1時(shí)只分離出最高頻的一半N6時(shí)分離到原頻率范圍的1/64。但信號(hào)的有效成分不一定那么低頻。如果采樣率很高信號(hào)的實(shí)際頻率又比較低層數(shù)再多也只是把低頻部分繼續(xù)細(xì)分對(duì)去噪沒有增益。我在處理工業(yè)現(xiàn)場(chǎng)振動(dòng)數(shù)據(jù)時(shí)總結(jié)過一個(gè)大致的參考采樣率在1000Hz左右信號(hào)主頻在50Hz以下5層足夠采樣率5000Hz主頻在100Hz以內(nèi)7層比較合適但如果要提取的特征是高頻沖擊成分反而應(yīng)該保持在4層以內(nèi)避免沖擊特征被閾值化掉。6.2 硬閾值萬能論硬閾值在數(shù)學(xué)上的確保留了系數(shù)的原始幅值但它處理后的系數(shù)會(huì)產(chǎn)生不連續(xù)點(diǎn)重構(gòu)時(shí)容易引發(fā)人為的局部振蕩。很多人發(fā)現(xiàn)去噪后波形在沖擊點(diǎn)附近出現(xiàn)細(xì)小的波紋多半就是硬閾值造成的。實(shí)踐建議先用軟閾值看效果如果具體場(chǎng)景需要保留準(zhǔn)確的信號(hào)幅值比如振動(dòng)幅值分析再考慮硬閾值配合更高一層的平滑處理。6.3 忘了歸一化信號(hào)單位不同、量綱不同直接影響閾值的幅值。比如信號(hào)幅值范圍是0~1噪聲標(biāo)準(zhǔn)差是0.1但另一個(gè)信號(hào)是0~1000那么同樣的噪聲標(biāo)準(zhǔn)差算出來的閾值完全不同。去噪前統(tǒng)一做歸一化是個(gè)好習(xí)慣x_min min(x); x_max max(x); x_norm (x - x_min) / (x_max - x_min); % 去噪 xd_norm wave_denoise(x_norm, ...); % 還原 xd xd_norm * (x_max - x_min) x_min;這樣處理還有一個(gè)好處閾值選擇時(shí)不容易被信號(hào)的絕對(duì)幅值干擾參數(shù)調(diào)一次可以復(fù)用到不同的數(shù)據(jù)集上。7. 為什么用“不同小波基分解層數(shù)”這個(gè)組合做對(duì)比說回項(xiàng)目本身。很多人直接用wden(s, sqtwolog, s, mln, N, sym6)一行代碼就把小波去噪跑完了但完全不理解為什么選擇sym6為什么是5層。這個(gè)項(xiàng)目的價(jià)值恰恰在于把“正確性”的假設(shè)打破讓你親手對(duì)比不同配置下的去噪結(jié)果差異。從影響權(quán)重來看小波基和分解層數(shù)這兩個(gè)因素占去噪效果的80%。閾值規(guī)則也很重要但它更像是微調(diào)。先從這兩個(gè)因素入手建立判斷標(biāo)準(zhǔn)后期遇到任何新的去噪任務(wù)都能快速定位參數(shù)范圍這才是做這個(gè)對(duì)比實(shí)驗(yàn)的最終目的。我在跑這個(gè)實(shí)驗(yàn)時(shí)遇到過一個(gè)小波基選擇的有趣現(xiàn)象對(duì)于同一個(gè)含噪信號(hào)db2的去噪結(jié)果中信號(hào)細(xì)節(jié)損失嚴(yán)重但sym2表現(xiàn)就明顯好很多盡管它們消失矩相同。原因就在于sym2近似對(duì)稱相位偏移更小。這個(gè)細(xì)節(jié)如果你不跑對(duì)比實(shí)驗(yàn)光看理論分析根本想不到。這也解釋了為什么我強(qiáng)烈建議不要只盯著某個(gè)小波基用到底——數(shù)據(jù)會(huì)告訴你哪個(gè)更合適。代碼的話全部函數(shù)可以按這個(gè)結(jié)構(gòu)組織gen_test_signal.m生成含突變分量和周期分量的測(cè)試信號(hào)。wave_denoise.m自定義去噪函數(shù)支持小波基、層數(shù)、閾值模式、軟硬閾值切換。run_basis_compare.m對(duì)比不同小波基。run_level_compare.m對(duì)比不同分解層數(shù)。run_auto_search.m自動(dòng)掃參輸出最優(yōu)參數(shù)組合。如果要把去噪應(yīng)用到自己的數(shù)據(jù)上只需要改加載信號(hào)那一行把仿真信號(hào)換成你的實(shí)際數(shù)據(jù)其余邏輯完全復(fù)用。最后分享一個(gè)非常實(shí)用的小技巧在判斷去噪效果時(shí)不要只看SNR或RMSE的數(shù)值一定要同時(shí)繪制“去噪誤差曲線”s - xd。理想誤差看起來應(yīng)該像隨機(jī)噪聲——沒有明顯的周期性成分也沒有集中的大偏差。如果誤差曲線出現(xiàn)明顯的低頻波動(dòng)說明去噪過程把有效信號(hào)的一部分也濾掉了這比任何指標(biāo)都直觀。我實(shí)測(cè)過很多數(shù)據(jù)集誤差曲線一旦出現(xiàn)“波浪形”哪怕SNR數(shù)值再高這個(gè)參數(shù)組合也不可用。這也是為什么每次做完自動(dòng)掃參后我都不會(huì)直接選最高SNR的那組參數(shù)而是把它和次高SNR參數(shù)的誤差曲線都畫出來肉眼對(duì)比后再做決定。這套方法論希望你也能用在自己的項(xiàng)目里。