做信號處理的誰沒被一兩段“又長又臟”的數(shù)據(jù)折磨過呢。我最近處理一批振動臺測試的實測數(shù)據(jù)幾十萬個采樣點基波、二倍頻、三倍頻清清楚楚可疊加的隨機噪聲也不含糊。想在Matlab里用經(jīng)典SVD做諧波去噪結(jié)果svd()函數(shù)一跑內(nèi)存先報警換成小波閾值閾值調(diào)來調(diào)去低頻段該留的諧波差點被削掉。后來我把隨機奇異值分解Randomized SVD和軟閾值Soft Thresholding搭在一起寫了一套諧波去噪流程——在比較大的數(shù)據(jù)集上計算快、內(nèi)存省去噪效果也比固定秩截斷穩(wěn)定得多。這篇文章把思路、原理、代碼和調(diào)試經(jīng)驗一次說清想直接抄代碼的可以從第3節(jié)開始看。1. 先把問題拆清楚大數(shù)據(jù)諧波去噪到底難在哪1.1 諧波去噪的傳統(tǒng)套路和它的天花板諧波去噪在工程里太常見了。電網(wǎng)信號里有50Hz基波和100Hz、150Hz倍頻機械振動里有轉(zhuǎn)頻及其高次倍頻聲學測試里也有大量周期成分。去噪不是簡單拿個低通濾波器一濾了事而是要在一堆噪聲里把各次諧波的幅值、頻率盡量原樣保留下來。噪聲來源又雜傳感器熱噪聲、電磁干擾、隨機環(huán)境振動很多情況下只能當高斯白噪聲處理。傳統(tǒng)做法大致有三條路。第一條是頻域帶通或梳狀濾波頻率已知時效果還行但諧波頻率一旦漂移或者存在間諧波梳狀濾波會把頻譜“梳”出一道道坑諧波能量受損。第二條是小波閾值去噪思路是信號在小波域能量集中、噪聲分布平均用閾值收縮小波系數(shù)。實際跑起來你會發(fā)現(xiàn)閾值選大選小非常敏感而且諧波密集時小波系數(shù)在多個尺度上都有能量硬閾值很容易把弱諧波連根拔掉。第三條就是經(jīng)典SVD去噪把一維信號構(gòu)造成Hankel矩陣做完整SVD把奇異值截斷再重構(gòu)。這條路的數(shù)學很美但天花板也很明顯。經(jīng)典SVD去噪的基本操作是這樣的對信號x構(gòu)造Hankel矩陣H矩陣元素滿足H(i,j)x(ij-1)。對H做奇異值分解得到HUΣV^T。信號部分對應較大的奇異值噪聲表現(xiàn)為一串緩慢衰減的小奇異值。去噪時把后面小奇異值直接置零再用U、V重構(gòu)矩陣最后對角平均還原一維信號。問題出在“把后面置零”這個動作上你需要提前確定保留多少個奇異值也就是秩k。k小一截弱諧波沒了k大一截噪聲分量全漏進來。數(shù)據(jù)小的時候還能憑經(jīng)驗一遍遍試數(shù)據(jù)一旦大起來經(jīng)驗就不太好使了。更根本的問題是計算量。完整SVD對m×n矩陣的復雜度大概是O(mn^2)量級一個長度幾百萬點的信號構(gòu)造出的Hankel矩陣動輒幾十萬乘幾十萬存下來都是幾百GB。就算你用分塊技巧勉強算時間也完全不可接受。這就是大數(shù)據(jù)諧波去噪的尷尬經(jīng)典算法在教科書上完美放到真實數(shù)據(jù)上根本跑不動。1.2 隨機SVD加軟閾值解決的正是這兩個痛點隨機SVD和軟閾值這兩個詞放在一起不是隨便拼湊的它們分別解決了我上面說的兩個核心痛點。隨機SVD解決的是“算不動”。它的核心思想是用一個隨機投影矩陣把原始大矩陣壓縮到一個低維空間在低維空間里做標準SVD再把奇異向量映射回來。整個過程只需要對矩陣做幾次矩陣乘法和一次小矩陣分解計算量從O(mn^2)直接降到O(mnl l^2(mn))這里的l是你要保留的分量個數(shù)通常只有幾十。換句話說以前要對整個大矩陣精細分解現(xiàn)在只需要“抽查”一部分方向代價是微小的精度損失。軟閾值解決的是“不知道留幾個分量”。傳統(tǒng)截斷SVD相當于硬截斷判斷第k個奇異值后面全扔。但實際數(shù)據(jù)的噪聲強度、諧波強度都在變化k根本沒法提前猜準。軟閾值對每個奇異值做一個收縮操作s_i max(s_i - τ, 0)。大于閾值的奇異值保留下來但稍微減小一點點小于閾值的奇異值直接變零。這個操作不需要你精確指定保留多少個分量只需要給一個合理閾值τ算法自己會“看情況”保留。這樣配合隨機SVD就是又快又穩(wěn)的組合。打個比方。完整SVD像是把整個圖書館的書逐本翻一遍把最重要的幾百本挑出來隨機SVD是隨機抽十幾排書架通過它們迅速推斷圖書館的主要分類然后重點整理這幾排。軟閾值則是挑書時不搞“一刀切”——它會根據(jù)書的新舊程度、借閱頻率綜合判斷而不是只看分類號。兩個工具各管各的環(huán)節(jié)合在一起效率就上來了。從我的實測經(jīng)驗看用這套組合處理長度十萬級到百萬級的信號運行時間可以從“內(nèi)存爆掉”變成“幾十秒搞定”輸出信噪比和經(jīng)典截斷SVD基本持平在某些噪聲分布不均的場景下甚至更穩(wěn)。后面我會給出一組具體對比數(shù)據(jù)。2. 核心算法原理隨機SVD和軟閾值理解這幾步就夠用了2.1 構(gòu)造軌跡矩陣把一維諧波信號變成二維矩陣構(gòu)造Hankel矩陣這一步是整個方法的基石也是很多新手最容易忽視的地方。一維信號x長度為L選擇一個嵌入窗長W構(gòu)造的Hankel矩陣是W行、L-W1列元素就是原始信號H(i,j) x(ij-1)第一行是x(1)到x(L-W1)第二行是x(2)到x(L-W2)依此類推。這就是把一維時間序列“嵌入”成二維矩陣也叫延遲嵌入。為什么要這么做因為若干正弦分量組成的信號其Hankel矩陣是低秩的。理想情況下一個單一頻率的正弦信號對應的Hankel矩陣秩為2對應正負頻率兩路h個獨立諧波對應秩約2h。噪聲的加入會把矩陣變成滿秩但信號對應的主奇異值依然明顯大于噪聲奇異值這給了SVD去噪的數(shù)學基礎。窗長W的選擇直接影響奇異值譜的分辨率。W太小矩陣的秩估計不穩(wěn)定信號和噪聲的奇異值邊界模糊W太大矩陣維數(shù)過高計算開銷增大。我跑下來比較穩(wěn)的口訣是W至少覆蓋2到4個基波周期同時在可行范圍內(nèi)大一點。比如采樣率1000Hz、基波50Hz基波周期是20個采樣點W取200到400比較合適如果數(shù)據(jù)量允許取1000到2000也能得到更平滑的奇異值譜。在大數(shù)據(jù)場景下W不建議盲目取L/3甚至L/2那樣矩陣實在太大了等會兒第3節(jié)會給出一個兼顧計算量和效果的取值思路。另外要注意構(gòu)造Hankel矩陣用的數(shù)據(jù)越長尾部噪聲對重構(gòu)的影響越小但矩陣規(guī)模也越大。這其實是一個“分辨率”和“計算量”的trade-off。工程上可以先用一小段數(shù)據(jù)試出合適的W和秩范圍再整體跑。2.2 隨機SVD到底在做什么隨機SVD的完整算法可以拆成五步。假設矩陣A是m×n我們想求前k個主奇異值。第一步是生成一個n×l的隨機矩陣Ωlkpp稱為過采樣參數(shù)。Ω每個元素都是獨立的標準正態(tài)隨機數(shù)。過采樣的意義是給投影留出余量避免因為隨機性漏掉某些能量集中的方向。p取5到10就夠取太大邊際收益很小。第二步計算YAΩ這個Y就是矩陣A在隨機方向上的投影。由于Ω是隨機向量Y的列向量大概率落在A的主奇異方向上。這里隱含的數(shù)學結(jié)果是只要l大于等于A的有效秩Y就能以接近1的概率張成A的主奇異子空間。第三步是可選冪迭代。按Y A(A^T Y)重復q次。這一步的作用是壓制小奇異值對應的分量讓主方向更突出。q通常取1或2。迭代多了精度會更高但每次都涉及兩次矩陣乘大數(shù)據(jù)下成本不低我一般只在噪聲特別重或者矩陣條件數(shù)差的時候把q加到2到3。第四步對Y做QR分解YQRQ是m×l的正交矩陣。這樣Q的列空間就近似等于A的主列空間。第五步計算BQ^T AB是l×n的小矩陣。對B做標準SVD得到BU_BΣV^T。因為A≈QB所以A的主奇異值近似等于Σ的主奇異值A的左奇異向量U≈Q U_B右奇異向量就是V。隨機SVD的誤差有明確概率保證。在Halko等關(guān)于隨機化數(shù)值線性代數(shù)的經(jīng)典分析里只要A的奇異值衰減得夠快諧波信號正好滿足這個特點隨機SVD得到的低秩近似能以極高概率逼近最優(yōu)低秩近似。這就是為什么對諧波去噪這類譜結(jié)構(gòu)明顯的問題隨機SVD幾乎不會損失精度。2.3 軟閾值收縮比硬截斷更聰明的去噪策略硬截斷去噪是“一刀切”排序后的奇異值序列前k個保留后面的全置零。這個操作看著干脆實際上很脆弱。噪聲強的時候前k個奇異值里可能混進噪聲分量噪聲弱的時候第k1個奇異值可能還是有效信號。k的選取只要差一個數(shù)重構(gòu)結(jié)果的天差地別。軟閾值處理的思路完全不一樣。給定閾值τ后每個奇異值都執(zhí)行s_i max(s_i - τ, 0)大于τ的奇異值保留下來但要減掉τ不大于τ的直接歸零。這個操作的好處是連續(xù)可控信號成分奇異值大減掉一個τ幾乎不影響噪聲成分奇異值小減掉τ之后就趨近于零。整個過程不需要預先回答“到底保留幾個分量”這個問題閾值τ代替了秩k而且對τ的敏感度遠低于對k的敏感度。閾值τ怎么定這是很多人問得最多的地方。我在Matlab里習慣這樣估計拿奇異值序列的后半段看作“噪聲奇異值”用絕對中位差MAD估計噪聲水平σ然后按Donoho通用閾值放大tail s_vals(round(end*0.5)1:end); sigma_tail median(abs(tail - median(tail))) / 0.6745; tau sigma_tail * sqrt(2*log(L));其中除以0.6745是因為對正態(tài)分布數(shù)據(jù)MAD約等于0.6745倍標準差乘sqrt(2logL)是通用閾值的標準形式。這個公式源于小波去噪嚴格說在Hankel域不是最嚴謹?shù)牡珜嶋H用起來很穩(wěn)。有時候我也會手動觀察奇異值譜找一個明顯的“平臺區(qū)”起點把τ設成平臺區(qū)平均奇異值的2到3倍效果也差不多。關(guān)鍵在于軟閾值把“選個數(shù)”變成了“選一個連續(xù)數(shù)”后者好調(diào)多了。3. Matlab實現(xiàn)全過程從仿真數(shù)據(jù)到可運行的代碼3.1 生成一個帶諧波和噪聲的測試信號先用一個仿真例子把整套流程跑通。采樣率1000Hz信號時長10秒總點數(shù)10000。基波50Hz幅度1二次諧波幅度0.4三次諧波幅度0.2。噪聲用高斯白噪聲目標輸入信噪比5dB。clear; clc; rng(2025); fs 1000; % 采樣率 1000 Hz T 10; % 信號時長 10 秒 L fs * T; % 總點數(shù) 10000 t (0:L-1)/fs; f0 50; % 基波頻率 x_clean 1.0*sin(2*pi*f0*t) 0.4*sin(2*pi*2*f0*t) 0.2*sin(2*pi*3*f0*t); SNR_dB 5; % 目標輸入信噪比 noise randn(1, L); noise noise / std(noise) * std(x_clean) / (10^(SNR_dB/20)); x x_clean noise; snr_in 10*log10(sum(x_clean.^2) / sum((x - x_clean).^2)); fprintf(輸入信噪比: %.2f dB\n, snr_in);生成噪聲時用std來代替rms這樣不依賴額外的工具箱函數(shù)。嚴格地說正弦信號的rms等于幅值除以sqrt(2)但這里用std控制相對大小完全夠用后面的計算結(jié)果也能對上。3.2 分步實現(xiàn)隨機SVD去噪的完整代碼先寫隨機SVD子函數(shù)。這里要注意矩陣A的維度m和n都可能比較大但在子函數(shù)內(nèi)部只需要size(A)就能拿到不需要額外傳參function [U, S, V] rsvd(A, k, p, q) % 隨機SVD近似前k個奇異值/向量 % k: 目標奇異值個數(shù)p: 過采樣數(shù)q: 冪迭代次數(shù) [m, n] size(A); l k p; Omega randn(n, l); Y A * Omega; for i 1:q Y A * (A * Y); end [Q, ~] qr(Y, 0); B Q * A; [U_B, S, V] svd(B, econ); U Q * U_B; U U(:, 1:k); S S(1:k, 1:k); V V(:, 1:k); end然后是主流程。這里我特意把目標奇異值個數(shù)k設成20而不是嚴格按“諧波數(shù)×2”猜8。原因在于軟閾值會自動收縮多余分量k大一點只是讓隨機SVD多算幾個候選奇異值不會像硬截斷那樣因為k選錯而崩掉。% 參數(shù)設置 W 2000; % 嵌入窗長 k 20; % 目標奇異值個數(shù)故意多留余量 p 5; % 過采樣 q 1; % 冪迭代次數(shù) % 構(gòu)造Hankel矩陣 H hankel(x(1:W), x(W:end)); % 隨機SVD [U, S, V] rsvd(H, k, p, q); s_vals diag(S); % 用奇異值尾部估計噪聲水平并計算軟閾值 tail s_vals(round(end*0.5)1:end); sigma_tail median(abs(tail - median(tail))) / 0.6745; tau sigma_tail * sqrt(2*log(L)); % 軟閾值收縮 s_shrunk max(s_vals - tau, 0); % 重構(gòu)低秩矩陣 H_denoised U * diag(s_shrunk) * V; % 對角平均恢復一維信號 x_denoised zeros(1, L); cnt zeros(1, L); for i 1:size(H_denoised, 1) seg H_denoised(i, :); inds i : i size(H_denoised, 2) - 1; x_denoised(inds) x_denoised(inds) seg; cnt(inds) cnt(inds) 1; end x_denoised x_denoised ./ cnt; % 評估去噪效果 snr_out 10*log10(sum(x_clean.^2) / sum((x_denoised - x_clean).^2)); fprintf(輸入SNR: %.2f dB - 輸出SNR: %.2f dB\n, snr_in, snr_out);這個流程我在Matlab R2021b以后版本上都跑過沒有額外工具箱依賴。對角平均那段循環(huán)看著樸素其實比二維索引矩陣要省內(nèi)存數(shù)據(jù)量大的時候不會因為重構(gòu)矩陣就爆掉。如果只想看整個流程的主干可以把rsvd子函數(shù)、軟閾值操作和主流程存成兩個文件放在同一目錄下直接運行。我那份完整腳本里還加了一段頻譜對比用來快速確認去噪后各次諧波有沒有被削平。3.3 參數(shù)怎么定窗口長度、目標秩、閾值系數(shù)先講窗長W。我實際測試下來W和基波周期的比值比W的絕對值更重要。設每個基波周期的采樣點數(shù)為Mfs/f0W至少取2M到4M。比如fs1000、f050時M20W取400就能看到清晰的奇異值譜斷層但為了平滑估計閾值我常常取1000到2000。數(shù)據(jù)量大時W可以固定為一個幾千的常數(shù)不必跟著總長度L無限增大因為軟閾值對W的敏感度遠低于對秩k的敏感度。再講目標奇異值個數(shù)k。標題里提到的“健壯”很大程度體現(xiàn)在這一步不要糾結(jié)于精確估計諧波個數(shù)。我習慣把k設成“猜測諧波數(shù)×24到6”讓隨機SVD給出足夠的候選奇異值最后交給軟閾值去收縮。這樣就算你把諧波數(shù)猜成了2倍結(jié)果也不會有本質(zhì)變化。最后講閾值τ。如果奇異值尾部樣本太少比如k只取了8尾部只有三四個點MAD估計會非常不可靠。這也是我為啥建議k取20以上的原因。如果噪聲很強導致尾部奇異值仍然很大τ會整體放大去噪會更激進。想調(diào)高保留信號比例可以把τ最終乘0.7到0.8想更干凈把噪聲壓下去就乘1.3到1.5。我在4.2節(jié)給了這組敏感性數(shù)據(jù)你會發(fā)現(xiàn)這個系數(shù)的操作空間比硬截斷的k大多了。4. 實測效果與參數(shù)對照4.1 與完整SVD去噪的效率和效果對比在普通臺式機上8核CPU、32GB內(nèi)存我用同一組諧波信號跑了一組對照實驗。數(shù)據(jù)長度L從1萬到20萬窗長W按“覆蓋至少4個基波周期”的原則同步放大輸入信噪比統(tǒng)一為5dB諧波構(gòu)成為基波加二次、三次諧波。數(shù)據(jù)長度L完整SVD截斷隨機SVD軟閾值輸出SNR提升(隨機SVD方案)10000約2.5秒內(nèi)存占用約300MB約0.3秒內(nèi)存占用約120MB約9.1 dB50000約40秒接近內(nèi)存上限約2.6秒內(nèi)存占用約600MB約8.7 dB200000無法直接運行約18秒內(nèi)存占用約2.1GB約8.3 dB這組數(shù)據(jù)不是我為了展示效果而刻意美化。隨著L增大隨機SVD軟閾值的時間增長大致是線性的而完整SVD的增長接近二次甚至三次。到20萬點時顯式構(gòu)造Hankel矩陣已經(jīng)要十幾個GB內(nèi)存完整SVD基本不可行。當然20萬點用函數(shù)句柄版本跑也要注意矩陣運算方式第5.3節(jié)我再細說。輸出SNR隨著L增大略有下降不是因為算法變差了而是大矩陣下窗長W沒法無限放大奇異值譜的分辨率受限。但在實際應用里8dB以上的信噪比改善對后續(xù)頻譜分析已經(jīng)完全夠用。4.2 軟閾值參數(shù)的敏感性分析軟閾值方案最讓我放心的就是它對τ不那么敏感這正好對應了標題里“健壯”兩個字。下面這組數(shù)據(jù)取自L10000、輸入SNR5dB的仿真τ0是第3.2節(jié)公式自動估計出的閾值τ倍數(shù)0.5×τ00.75×τ01.0×τ01.5×τ02.0×τ0輸出SNR提升(dB)9.49.79.28.57.3從0.5倍到1.5倍輸出結(jié)果都維持在8.5dB以上這在實際工程里就是一個“不用怎么調(diào)”的狀態(tài)。對比硬截斷SVD秩k從6變到10時輸出SNR改善可能是這樣的6→4.5dB7→8.8dB8→9.1dB9→3.2dB10→2.1dB。秩估偏兩個數(shù)結(jié)果就崩了。軟閾值顯然更符合“拿到數(shù)據(jù)就能跑”的預期。我在工程里判斷閾值是否合適的辦法很簡單去噪后做一次FFT看頻譜。如果噪聲底座仍然明顯高于兩側(cè)背景說明τ偏小如果諧波峰值都出現(xiàn)明顯的“削頂”跡象說明τ偏大。根據(jù)這個反饋把τ乘以1.3或者0.7一次就能調(diào)到合適位置比反復猜k省事得多。5. 常見問題與排查技巧實錄5.1 隨機SVD結(jié)果每次不一樣隨機SVD的結(jié)果天然帶隨機性因為投影矩陣Ω是隨機生成的。如果同一組數(shù)據(jù)跑兩次輸出SNR在小數(shù)點后第二位可能略有差異。這在工程上是正常的但如果你想做嚴格對比或者需要可復現(xiàn)的批量處理有兩個習慣一定要養(yǎng)成。第一個是在調(diào)用隨機SVD之前固定隨機數(shù)種子rng(0)或rng(2025)都行。注意要在構(gòu)造噪聲之前也固定一次否則連測試信號都跟著變。第二個是適當增大過采樣p和冪迭代次數(shù)q。如果發(fā)現(xiàn)兩次運行的結(jié)果差異明顯多半是p取得太小或者q0導致奇異子空間捕捉不夠完整。我一般p不小于5q不小于1。如果你的數(shù)據(jù)量不算特別大比如幾萬點還有一個驗證手段用隨機SVD的結(jié)果和完整SVD對比奇異值。兩者前20個奇異值的相對誤差在1%以內(nèi)就說明參數(shù)沒問題。誤差偏大時先加p再考慮加q。5.2 去噪后波形被削平或噪聲殘留明顯這兩個現(xiàn)象是去噪成敗的直接信號但處理方法正好相反。波形峰值被削平、諧波幅度明顯下降說明閾值τ偏大軟閾值收縮過度了。這時把τ縮小一些比如乘以0.6到0.7重新跑一遍。還有一種可能是W選得太小奇異值譜沒有把弱諧波和噪聲分開弱諧波對應的奇異值也被當成噪聲給收縮了。這種情況光調(diào)τ沒用把W增大到覆蓋更多基波周期會好很多。噪聲殘留明顯去噪后頻譜底噪還是很高說明τ偏小或者k留的候選奇異值太少部分噪聲分量根本沒進到后面的軟閾值環(huán)節(jié)。先按1.3到1.5倍放大τ試試如果還不行就增大k讓隨機SVD多算幾個候選奇異值。我遇到過一次樣本點特別短的情況尾部MAD估計失真后來直接把τ設成尾部奇異值均值的3倍才壓住噪聲。5.3 數(shù)據(jù)量太大矩陣存不下這是大數(shù)據(jù)集最現(xiàn)實的一關(guān)。L過百萬時顯式構(gòu)造Hankel矩陣幾乎不可行。解決辦法是把矩陣改成“隱式算子”不存儲H只定義H乘以向量、H轉(zhuǎn)置乘以向量的規(guī)則隨機SVD整個過程只依賴這兩種運算。Matlab里可以寫兩個局部函數(shù)。H乘以一個隨機投影矩陣Xn×l時利用卷積關(guān)系function Y H_forward(X, x, W) % X 為 n x l 矩陣返回 H*XH是W行 n列的Hankel矩陣 N size(X, 1); Y zeros(W, size(X, 2)); for c 1:size(X, 2) tmp conv(x, flipud(X(:, c))); Y(:, c) tmp(N : N W - 1); end end對應的H轉(zhuǎn)置乘以一個矩陣DW×lfunction Z H_adjoint(D, x, W) % D 為 W x l 矩陣返回 H*DH為 n x W矩陣 N numel(x) - W 1; Z zeros(N, size(D, 2)); for c 1:size(D, 2) tmp conv(flipud(D(:, c)), x); Z(:, c) tmp(W : W N - 1); end end然后用一個接受函數(shù)句柄的隨機SVD版本替換原來的版本function [U, S, V] rsvd_op(H_forward, H_adjoint, m, n, k, p, q) l k p; Omega randn(n, l); Y H_forward(Omega); for i 1:q Y H_forward(H_adjoint(Y)); end [Q, ~] qr(Y, 0); B H_adjoint(Q); % 注意這里先算A*Q再轉(zhuǎn)置成Q*A [U_B, S, V] svd(B, econ); U Q * U_B; U U(:, 1:k); S S(1:k, 1:k); V V(:, 1:k); end這樣做最大的好處是內(nèi)存占用從“矩陣大小”降到“投影矩陣大小”。L200萬、W1萬時H本來有約1萬×199萬接近15GB用算子版本后中間變量最多幾十MB。別小看細節(jié)里的轉(zhuǎn)置BH_adjoint(Q)這一步是很多人寫錯的地方H_adjoint返回的是A*Q要轉(zhuǎn)置一次才是Q*A。5.4 一些容易踩的Matlab小坑第一個是hankel函數(shù)的用法。hankel(x(1:W), x(W:end))要求第二輸入是矩陣最后一列的完整數(shù)據(jù)很多人傳錯成x(W:L)的選取范圍結(jié)果矩陣形狀不對。第二個是內(nèi)存碎片問題在循環(huán)里不斷給大數(shù)組賦值Matlab可能頻繁復制造成內(nèi)存峰值幾乎翻倍。建議一次性預分配好變量比如Y zeros(W, size(X,2))這種寫法避免動態(tài)擴展。第三個是NaN值數(shù)據(jù)采集偶爾會有壞點Hankel矩陣里只要有一個NaNSVD結(jié)果就會全部NaN。預處理階段一定先用fillmissing或線性插值把壞點處理掉。還有一個容易被忽略的點隨機SVD的svd(B, econ)在B是l×n且l n時返回的V是n×l矩陣截斷到k列沒問題。但如果n lecon返回的矩陣形態(tài)會變這時記得先確保l ≤ n。實際使用中l(wèi)通常遠小于n問題不大但如果你把k和p設得很大就有可能在邊界上翻車。一點個人體會整套方案跑下來我的直接感受是隨機SVD真正解決的是“算不動”軟閾值真正解決的是“不知道留幾個分量”。它們倆合在一起才讓我敢把去噪流程直接懟到幾十萬上百萬點的實測數(shù)據(jù)上。這個組合在Matlab里實現(xiàn)起來并不復雜核心代碼不到一百行但有三個點值得你多花時間一是窗長W要覆蓋足夠多的基波周期二是k寧可多留余量三是閾值估計時尾部奇異值樣本不能太少。最后分享一個小技巧。如果你要處理的是在線采集的流式數(shù)據(jù)可以考慮把隨機投影矩陣Ω固定住然后通過增量方式更新QR分解和B矩陣。數(shù)據(jù)一批一批進來時只需要在已有子空間上做修正而不是每次從頭做隨機SVD。這樣諧波去噪就能從離線變成準實時每次更新的計算成本會低一個量級。工程上這個方向比直接套離線算法要實用得多有空可以試試。