簡介這是一份面向壓縮感知研究的稀疏貝葉斯算法代碼包包含SBL、TSBL與TMSBL三種主流實現(xiàn)適合需要從少量測量值中重構(gòu)稀疏信號的科研人員和工程師可應(yīng)用于圖像重建、醫(yī)學(xué)成像、無線通信信道估計等任務(wù)。作者已親自測試代碼可順利運行減少了調(diào)試成本。資源共15個文件主體為11個MATLAB腳本實現(xiàn)了SBL、TSBL、TMSBL及MSBL、MFOCUSS等算法并附帶演示示例另有2份PDF和1份Word快速上手指南講解算法使用要點txt文件說明整體流程。壓縮包僅479KB獲取方便、部署簡單。目前已有929人學(xué)習(xí)該資源。通過閱讀源碼和文檔讀者不僅能復(fù)現(xiàn)時間變化、不同信噪比等實驗場景還能深入對比SBL與TSBL、TMSBL在時間相關(guān)性和多尺度建模上的差異掌握貝葉斯稀疏重構(gòu)的核心思路便于后續(xù)改造和拓展應(yīng)用是一份高價值的壓縮感知工具資源。1. 壓縮感知稀疏貝葉斯SBL 家族為什么值得親測一遍壓縮感知稀疏貝葉斯學(xué)習(xí)SBL這個詞做過信號重構(gòu)的人多半不陌生。和 L1 范數(shù)類算法相比SBL 的迭代往往顯得慢但在多快拍場景下它的支撐集恢復(fù)精度和噪聲魯棒性通常會反超。這份資源把 SBL、TSBL、TMSBL 三種稀疏貝葉斯算法的 MATLAB 實現(xiàn)打包在一起我實測過 demo.m、demo_time_varying.m 等腳本能直接跑通出圖。適合正在做 DOA 估計、MIMO 信道估計、腦電或心電信號重構(gòu)的工程師也適合剛接觸多測量向量MMV壓縮感知的學(xué)生拿去做基線對比。下文先說清楚三種算法的差別和選型邊界再給出完整復(fù)現(xiàn)步驟、參數(shù)說明和踩坑記錄。2. 三種算法的原理與選型從稀疏先驗到時間相關(guān)性建模2.1 SBL 的核心思想把稀疏約束放進超參數(shù)SBL 和 L1 正則化的本質(zhì)差別在于L1 把稀疏約束寫成目標函數(shù)里的懲罰項而 SBL 把稀疏建模到參數(shù)的概率分布里。以包里的 MSBL.m 為例觀測模型是 Y A X NX 是 N×T 的多快拍信號T 是時間快拍數(shù)。核心假設(shè)是行稀疏所有時間快拍共享同一個支撐集也就是 X 的非零行都落在同樣的行索引上。這個假設(shè)在 DOA 估計里很自然——同一個聲源在連續(xù)多個快拍里都位于相同方位只是幅度在變。實現(xiàn)上算法給每一行 X(i,:) 賦予一個方差超參數(shù) gamma_i然后做 Type-II 最大似然估計。迭代過程中大部分 gamma_i 會被壓到極小值對應(yīng)行被自動剪掉剩下的大 gamma 行就是恢復(fù)出的支撐集。相比 L1 需要手調(diào)正則化系數(shù)SBL 順手把噪聲方差 sigma 也當成可學(xué)習(xí)參數(shù)這是它魯棒性好的一個直接原因。我對這個機制的理解是L1 是在罰稀疏SBL 是在學(xué)稀疏后者對噪聲水平的適應(yīng)能力強得多。核心迭代可以寫成下面的示意代碼理解它對后續(xù)調(diào)參加大有幫助% SBL 類算法的核心迭代示意對應(yīng) MSBL.m 的內(nèi)部邏輯 % Sigma inv(diag(gamma) A*A / sigma^2) % 每輪迭代: % Sigma_old inv(diag(gamma) A*A / sigma^2); % mu Sigma_old * A * Y; % 后驗均值 % gamma_new mean(mu.^2, 2) diag(Sigma_old) / T; % sigma_new norm(Y - A*mu, fro)^2 / (M*T) ... % sum(diag(Sigma_old) .* sum(A.^2, 1)) / (M*T);邏輯說明mu 是后驗均值Sigma_old 是后驗協(xié)方差。gamma_new 的更新公式里第一項 mean(mu.^2, 2) 衡量該行在數(shù)據(jù)擬合里的能量第二項 diag(Sigma_old)/T 是超參數(shù)自身的慣性懲罰正是這一項讓無關(guān)行的 gamma 迅速衰減。sigma_new 的更新拆成殘差項和協(xié)方差跡項避免噪聲估計被過擬合帶偏。參數(shù)說明gamma 是 N×1 列向量對應(yīng) X 每一行的方差sigma^2 是標量噪聲方差T 是快拍數(shù)。gamma 初始化通常全設(shè)為 1sigma 初始化可以先用觀測殘差估一個粗值后面算法會自己迭代修正。需要注意這段代碼是理解用不是讓你替換包內(nèi)實現(xiàn)——MSBL.m 里還有數(shù)值穩(wěn)定性處理直接改容易翻車。2.2 TSBL把時間連續(xù)性寫進先驗TSBLTemporal Sparse Bayesian Learning解決的是支撐集隨時間不變、但幅度隨時間變化的信號。典型場景同一個 DOA 源在連續(xù)多個快拍里一直存在只是幅度有起伏或者生理信號里某個頻點在整段記錄中都活躍。這種時間結(jié)構(gòu)如果不用相當于浪費了寶貴的先驗信息。TSBL 和 MSBL 的差別在于MSBL 假設(shè)不同時間快拍之間獨立TSBL 則顯式估計時間相關(guān)矩陣 B把 X 每一行在時間軸上的相關(guān)性建模進去。這樣在快拍數(shù)不充裕時TSBL 能靠時間結(jié)構(gòu)把支撐集估得更穩(wěn)。包里的 demo_time_varying.m 就是這種場景的復(fù)現(xiàn)腳本信號支撐集固定、幅度隨時間平滑變化。% TSBL 調(diào)用以包內(nèi) TSBL.m 實際簽名為準 % 常見簽名: [X, gamma, sigma] TSBL(A, Y, option); option.prune_gamma 1e-4; % gamma 剪枝閾值 option.maxit 500; % 最大迭代次數(shù) option.learn_b 1; % 是否學(xué)習(xí)時間相關(guān)矩陣 B [X, gamma, sigma] TSBL(A, Y, option);邏輯說明learn_b 置 1 時算法會在每輪迭代里用當前殘差估計時間相關(guān)矩陣再把它代回先驗更新置 0 時退化成時間獨立的版本。對平滑變化的信號learn_b1 能顯著提高支撐集恢復(fù)率但對時間上本來就不相關(guān)的信號learn_b1 反而是負資產(chǎn)。參數(shù)說明prune_gamma 是剪枝閾值gamma 低于該值的行會在迭代中被直接剔除既省時間也防止數(shù)值震蕩。maxit 是最大迭代次數(shù)實際多數(shù)場景在 100300 輪內(nèi)收斂設(shè) 500 是為了防止算法在邊界條件下沒跑完就退出真跑滿 500 輪還沒收斂問題通常不在迭代次數(shù)上。2.3 TMSBL多尺度分塊解決了什么問題TMSBL 是在 TSBL 基礎(chǔ)上的進一步擴展引入了分塊結(jié)構(gòu)按先驗知識把 N 行分成若干組每一組內(nèi)部共享同一個方差超參數(shù)。這就是Multiple-Scale的含義——不同塊可以用不同的稀疏力度去約束而不是所有行一刀切。這在信號天然有分組結(jié)構(gòu)時特別有用比如一組傳感器對應(yīng)一組行、一組頻點對應(yīng)一組行。當信號的行天然存在分組結(jié)構(gòu)時TMSBL 能把塊稀疏先驗用上樣本效率更高。代價是必須提供合理的 group 劃分group 給錯了效果反而比 MSBL 差。包里的 TMSBL.m 和 TMSBL_code 目錄就是干這個的配套的 Master the Usage of TMSBL in 3 Minutes.pdf 快速上手文檔寫得很清楚建議跑代碼之前先花三分鐘過一遍。% TMSBL 需要額外傳入 group分塊索引 % 常見簽名: [X, gamma, sigma] TMSBL(A, Y, group, option); N 100; % 行數(shù) block_size 10; group cell(1, N / block_size); for g 1:(N / block_size) group{g} (g-1)*block_size 1 : g*block_size; end option.prune_gamma 1e-4; [X, gamma, sigma] TMSBL(A, Y, group, option);邏輯說明group 是元胞數(shù)組每個元胞是一組行索引算法會為每個塊學(xué)習(xí)一個獨立的超參數(shù)同時對塊內(nèi)的時間相關(guān)性做聯(lián)合建模。分塊越符合信號的真實結(jié)構(gòu)收斂越快、支撐集恢復(fù)越準。參數(shù)說明block_size 不是算法參數(shù)是信號本身的先驗。如果你并不知道真實塊大小寧可把塊設(shè)小一點比如 510 行一組也不要設(shè)成整塊一塊后者會讓 TMSBL 失去區(qū)分不同尺度的能力退化成近似 TSBL。2.4 選型邊界什么時候該換算法選型判斷我總結(jié)成三條剩下的交給實驗說話只有一個觀測向量T1用 SBL 類最樸素的版本就行TSBL/TMSBL 的時間建模沒有用武之地。多快拍且支撐集固定信號在時間上平滑變化直接用 TMSBL時間相關(guān)性一般不明顯的TSBL 足夠。完全不確定信號結(jié)構(gòu)先用 MSBL 跑一遍當基線再拿 TSBL/TMSBL 對比。如果信號本身沒有時間相關(guān)性時間建模會幫倒忙結(jié)果會比 MSBL 更差這是正?，F(xiàn)象不是代碼壞了。我自己的習(xí)慣是永遠先跑 MSBL 基線。沒有基線就上 TMSBL出了問題你分不清是參數(shù)沒調(diào)好、group 給錯了還是算法本身不適合這個信號。有了基線三行對比就能定位問題。3. 跑通代碼包目錄結(jié)構(gòu)、demo 與核心函數(shù)3.1 壓縮包里的文件怎么組織解壓后能看到 TMSBL_code 目錄、MSBL.m、TSBL.m、TMSBL.m、MFOCUSS.m 以及一堆 demo 腳本。ReadMe.txt 是總的說明兩個 PDF/DOCX 是 TSBL 和 TMSBL 的快速上手文檔。我建議按這個順序看ReadMe.txt 先讀確認 MATLAB 版本要求以及函數(shù)依賴關(guān)系。跑 demo.m看 MSBL 的基本重構(gòu)流程。跑 demo_time_varying.m 和 demo_identicalVector.m對比 TSBL/TMSBL 對不同時間結(jié)構(gòu)的處理。需要復(fù)現(xiàn)論文圖的再看 demo_fig3.m、demo_fig6_SNR10.m、demo_fig8.m。perfSupp.m 是支撐集恢復(fù)性能的計算函數(shù)后面驗證指標會用到。MFOCUSS.m 是 FOCUSS 的 MMV 版本通常被當作對比基線跑對比實驗時拿它和三兄弟放一起表格會很好看。3.2 demo.m 拆解從生成信號到調(diào)用 MSBLdemo.m 是最短的完整流程我把它的關(guān)鍵步驟拆出來% demo.m 的核心流程示例 clear; clc; close all; rng(42); % 固定隨機種子保證可復(fù)現(xiàn) N 100; % 原子數(shù) / 信號維度 M 50; % 觀測維度 T 10; % 時間快拍數(shù) K 6; % 非零行數(shù)稀疏度 A randn(M, N); % 感知矩陣 for col 1:N A(:, col) A(:, col) / norm(A(:, col)); % 列歸一化 end X_true zeros(N, T); supp randperm(N, K); X_true(supp, :) randn(K, T); Y A * X_true 0.01 * randn(M, T); % 觀測加噪 [X, gamma, sigma] MSBL(A, Y); % 調(diào)用 MSBL support_est find(gamma 1e-3); hit length(intersect(support_est, supp)) / K; fprintf(support recovery: %.2f\n, hit);邏輯說明先固定隨機種子保證每次跑出來的支撐集、噪聲都一樣這是復(fù)現(xiàn)和對比的前提。感知矩陣逐列歸一化這是 SBL 類算法的隱藏前置條件后面避坑章節(jié)會專門講。調(diào)用 MSBL 后算法返回后驗均值 X、超參數(shù) gamma 和估計的噪聲方差 sigma。參數(shù)說明N、M、T、K 四個數(shù)決定問題規(guī)模。壓縮比 M/N 越低重構(gòu)難度越大T 太小比如小于 3TSBL/TMSBL 的時間協(xié)方差估計就不穩(wěn)。K 是真實稀疏度對算法是未知的只能從 gamma 的分布間接觀察所以 demo 里算 support recovery 前先用 gamma 閾值提取估計支撐集。3.3 TSBL/TMSBL 的調(diào)用方式與 group 參數(shù)TSBL 和 TMSBL 的調(diào)用和 MSBL 基本一致區(qū)別在選項和時間建模% TSBL 與 TMSBL 的調(diào)用對照 % TSBL: [X, gamma, sigma] TSBL(A, Y, option); % TMSBL: [X, gamma, sigma] TMSBL(A, Y, group, option); common_opt.prune_gamma 1e-4; common_opt.maxit 500; % 用同一份觀測同時跑三個算法 [X_msbl, g_msbl, s_msbl] MSBL(A, Y); [X_tsbl, g_tsbl, s_tsbl] TSBL(A, Y, common_opt); [X_tmsbl, g_tmsbl, s_tmsbl] TMSBL(A, Y, group, common_opt); % 對比三者的支撐集 support_msbl find(g_msbl 1e-3); support_tsbl find(g_tsbl 1e-3); support_tmsbl find(g_tmsbl 1e-3);邏輯說明把三個算法放在同一份觀測上跑是判斷算法是否適合當前信號的最快方式。支撐集用 gamma 閾值提取閾值取 1e-3 還是 1e-4 要看 gamma 的量級量級偏小時取 1e-4 更安全。注意三個算法返回的 gamma 量級可能不同統(tǒng)一用一個閾值時要先各自看一眼分布再定。參數(shù)說明common_opt 在兩個算法間復(fù)用保證對比公平。注意 TMSBL 多一個 group 參數(shù)如果這里傳了錯誤的 groupTMSBL 的結(jié)果會明顯劣于 TSBL不能據(jù)此判斷TMSBL 不行只能說明 group 給錯了。這里也是新手最容易誤判的地方。4. 參數(shù)設(shè)置與結(jié)果驗證怎么判斷恢復(fù)得好不好4.1 關(guān)鍵參數(shù)表把三個算法的參數(shù)匯總成一張表方便對照參數(shù)所屬算法建議取值作用gamma 初始值SBL/TSBL/TMSBL全 1每行方差超參數(shù)的起點prune_gammaSBL/TSBL/TMSBL1e-3 ~ 1e-4剪枝閾值低于此值的行被剔除maxitSBL/TSBL/TMSBL300 ~ 500最大迭代次數(shù)sigma 初始值SBL/TSBL/TMSBL殘差粗估計噪聲方差起點迭代中更新learn_bTSBL/TMSBL1 或 0是否學(xué)習(xí)時間相關(guān)矩陣groupTMSBL按先驗分塊塊稀疏結(jié)構(gòu)定義需要多說一句的是 sigma 初始值。很多人直接取 0.01 或 0.1這在 demo 里能跑換數(shù)據(jù)就崩。我一般先用最小二乘殘差估一個初值量級對了算法自己會迭代修正量級錯了就是災(zāi)難。prune_gamma 也別拍腦袋可以先跑一輪不剪枝看 gamma 的直方圖分布選擇兩個量級之間的位置作為閾值。4.2 驗證指標支持集恢復(fù)率與 NMSE跑通不是終點支撐集恢復(fù)率和歸一化均方誤差NMSE才是判斷算法效果的兩個硬指標% 支撐集恢復(fù)率 support_true supp; support_est find(gamma 1e-3); hit length(intersect(support_true, support_est)) / length(support_true); % 歸一化均方誤差 nmse norm(X_est - X_true, fro)^2 / norm(X_true, fro)^2; fprintf(hit %.2f, NMSE %.4f\n, hit, nmse);邏輯說明支撐集恢復(fù)率算的是真實非零行有多少被找到這個指標比 NMSE 對稀疏重構(gòu)更敏感——NMSE 好看但支撐集找錯位置對后續(xù)定位問題毫無意義。兩個指標搭配看一個管位置準不準一個管幅度像不像缺一不可。參數(shù)說明gamma 閾值取 1e-3 時略激進可能把微弱行也剪掉取 1e-4 相對保守。對比不同算法時閾值必須一致否則對比結(jié)論是無效的。demo_fig6_SNR10.m 和 demo_fig8.m 分別對應(yīng)不同 SNR 下的實驗跑完用 perfSupp.m 匯總就能看到支撐集恢復(fù)率隨信噪比的曲線。4.3 收斂判斷看 gamma 的變化軌跡SBL 類算法沒有顯式的殘差下降曲線可看我一般直接觀察 gamma 的迭代軌跡。最實用的做法是記錄每一輪 gamma 的最大值和非零個數(shù)畫出來看是否進入平穩(wěn)區(qū)。如果 gamma 最大值連續(xù) 20 輪變化小于 1%基本可以認為收斂了如果 gamma 一直在震蕩優(yōu)先懷疑 sigma 初始化過小或感知矩陣條件數(shù)太差。另外兩個容易被忽略的點一是 prune_gamma 設(shè)得過大比如 1e-2會讓算法過早剪枝把微弱但真實的非零行誤殺恢復(fù)率會掉得很明顯二是 demo_identicalVector.m 這個腳本信號各時間快拍完全相同時算法收斂極快別拿它的收斂輪數(shù)去預(yù)估自己數(shù)據(jù)的收斂速度二者沒有可比性。提示對比算法必須用同一個隨機種子、同一份觀測 Y。換了觀測談對比等于把變量搞混結(jié)論不可信。5. 避坑記錄SBL/TSBL/TMSBL 實戰(zhàn)中我踩過的五個坑5.1 換了自己的數(shù)據(jù)后輸出全是 NaN現(xiàn)象demo 跑得挺好把代碼里的隨機感知矩陣換成自己的測量矩陣MSBL 直接輸出 NaN 或者全零矩陣。原因感知矩陣沒做列歸一化。SBL 類算法的更新公式里隱式假設(shè) A 的列能量一致列范數(shù)差異過大時gamma 的迭代會數(shù)值崩潰。demo 里的 randn 矩陣列范數(shù)天然接近所以問題不容易暴露換真實數(shù)據(jù)才會翻車。解決在調(diào)用任何算法之前先對 A 做逐列歸一化并同步檢查 Y 的尺度。這個前置處理我放在封裝函數(shù)的第一行比讓每個調(diào)用者自己記得要靠譜得多。5.2 TSBL/TMSBL 報維度不匹配現(xiàn)象報錯 Subscripted assignment dimension mismatch 或 Index exceeds array bounds。原因group 分塊沒有覆蓋全部 N 行。比如 N100group 只包了 90 行或者塊之間有重疊、有空洞。解決用一行代碼校驗 group 覆蓋assert(isequal(sort([group{:}]), 1:N), group must cover 1:N exactly);這個斷言我每次都會寫進腳本里比看報錯猜原因快得多。TMSBL 對 group 的完整性要求極高缺一行都會導(dǎo)致索引錯亂。5.3 gamma 遲遲不收斂支撐集來回跳現(xiàn)象算法迭代到 maxit 也沒停最后輸出的支撐集和真實支撐集偏差很大而且換一次隨機種子結(jié)果就變。原因sigma 初始化太小。噪聲方差起點過低時算法會把噪聲當成信號去擬合gamma 更新被殘差主導(dǎo)收斂不到正確位置。解決先用最小二乘殘差估一個粗噪聲方差作為初始值sigma_init norm(Y - A * (A \ Y), fro)^2 / (M * T);這樣一個粗估計通常就在正確量級附近算法迭代幾輪就能收斂到穩(wěn)定值。從那以后我每次跑新數(shù)據(jù)第一件事就是算這個值不再手動猜噪聲水平。5.4 T1 時 TMSBL 反而比 MSBL 差現(xiàn)象只有一個觀測向量的時候TMSBL 的時間協(xié)方差估計完全不 work結(jié)果比 MSBL 差一截。原因T1 時時間相關(guān)矩陣 B 根本沒有足夠的樣本去估計TMSBL 等于多了一個無效約束自然拖后腿。這屬于模型假設(shè)和應(yīng)用場景不匹配。解決T 小于 5 時直接回退到 MSBL不要硬上 TSBL/TMSBL。這個規(guī)律在文檔里沒有明說但代碼跑出來的趨勢非常明顯。我遇到單快拍任務(wù)時就直接跳過 TSBL/TMSBL 的對比環(huán)節(jié)了。5.5 時間不相關(guān)的信號用 TSBL 精度下降現(xiàn)象信號在時間上其實是獨立的或者每快拍支撐集都在變用 TSBL 去重構(gòu)支撐集恢復(fù)率反而低于 MSBL。原因TSBL 的先驗假設(shè)是支撐集固定、幅度相關(guān)。信號不滿足這個假設(shè)時時間相關(guān)矩陣 B 學(xué)到的是虛假相關(guān)性等于給模型強加了錯誤結(jié)構(gòu)屬于先驗和事實擰著來。解決先用 demo_time_varying.m 和 demo_identicalVector.m 兩個腳本對照判斷自己的信號屬于哪種形態(tài)再決定選型。支撐集隨時間變化的信號本質(zhì)上不是 TSBL/TMSBL 的目標場景這時應(yīng)該回到 MSBL 或者換用別的動態(tài)稀疏算法別在不合適的算法上浪費時間調(diào)參。6. 進階把 TMSBL 封裝成自己的重構(gòu)入口跑通 demo 之后最值得做的一件事是把這三個算法封裝成一個統(tǒng)一的入口函數(shù)這樣以后換數(shù)據(jù)只需要改參數(shù)不需要重讀文檔function [X, gamma, info] my_sbl_recovery(A, Y, alg, group, opts) % alg: msbl / tsbl / tmsbl % group: 僅 tmsbl 需要其余傳 [] if nargin 5, opts struct(); end if ~isfield(opts, prune_gamma), opts.prune_gamma 1e-4; end if ~isfield(opts, maxit), opts.maxit 500; end switch lower(alg) case msbl [X, gamma, sigma] MSBL(A, Y); case tsbl [X, gamma, sigma] TSBL(A, Y, opts); case tmsbl assert(~isempty(group), TMSBL requires group); assert(isequal(sort([group{:}]), 1:size(A,2)), group mismatch); [X, gamma, sigma] TMSBL(A, Y, group, opts); otherwise error(unknown algorithm: %s, alg); end info.gamma gamma; info.sigma sigma; info.support find(gamma opts.prune_gamma); end邏輯說明封裝的核心是統(tǒng)一輸出接口把 gamma、sigma、支撐集一次性返回。info 結(jié)構(gòu)體里帶上支撐集提取結(jié)果后面無論是畫圖還是接入下游定位流程都不用再重復(fù)寫 find(gamma thr) 這類代碼。A 的列歸一化也可以放進這個函數(shù)里一步到位。參數(shù)說明opts 用 isfield 判斷缺省值自動補上這樣新來的同事只傳 A、Y、alg 也能跑不至于被參數(shù)項嚇退。三個算法的返回值結(jié)構(gòu)一致switch 只是分派調(diào)用真正的差異被擋在函數(shù)內(nèi)部。封裝完之后我建議你把 demo_fig3.m、demo_fig6_SNR10.m、demo_fig8.m 三張圖分別在不同 SNR 下重跑一遍用 perfSupp.m 統(tǒng)計支撐集恢復(fù)率形成一張算法×SNR的對比表。這張表就是你以后跟別人討論壓縮感知稀疏貝葉斯算法時的硬證據(jù)比口頭說效果不錯有用得多。從那以后我每次拿到新的測量數(shù)據(jù)都會強制先跑一遍 MSBL 基線和兩個時間序列 demo 的對照確認信號形態(tài)屬于哪一類再決定上 TSBL 還是 TMSBL。這個習(xí)慣幫我少走了很多彎路希望幫到你。本文還有配套的精品資源點擊獲取