簡介面向信號處理、語音分析與算法復(fù)現(xiàn)的學(xué)習(xí)者這份資源圍繞JADE盲源分離展開適合已有一定MATLAB基礎(chǔ)、希望理解盲信號分離原理并動手實(shí)現(xiàn)的研究生與工程技術(shù)人員。內(nèi)容從盲源分離的基本問題切入在源信號與混合方式均未知的條件下利用觀測信號的統(tǒng)計(jì)獨(dú)立性恢復(fù)原始信號并重點(diǎn)說明JADE依賴高階累積量常用四階也有三階方案與源信號非高斯性的思路。資源還提示了方法邊界除統(tǒng)計(jì)獨(dú)立外源信號中最多只能有一個高斯信號算法主要適用于統(tǒng)計(jì)獨(dú)立的非高斯信號對非白、非平穩(wěn)特性未作專門考慮。壓縮包約426KB文件總數(shù)與類型明細(xì)暫未提供從標(biāo)題看應(yīng)包含原理說明及MATLAB程序便于對照公式理解累積量計(jì)算、分離矩陣求解與結(jié)果驗(yàn)證。已有953人學(xué)習(xí)適合作為盲源分離入門與課程實(shí)驗(yàn)的參考材料。1. 三路信號混在一起時這套 JADE 的 matlab 算法能拆到什么程度實(shí)驗(yàn)室里那臺三通道采集設(shè)備采回來的信號永遠(yuǎn)是混在一起的三路傳感器同時錄同一個沖擊事件在三路里都留下痕跡誰的貢獻(xiàn)占多少說不清。想拆開深度學(xué)習(xí)的路子要先有標(biāo)注、要調(diào) GPU、要反復(fù)試網(wǎng)絡(luò)結(jié)構(gòu)而一套完整的 matlab 算法——JADEJoint Approximate Diagonalization of Eigenmatrices特征矩陣聯(lián)合近似對角化——核心代碼不到兩百行跑一次幾秒鐘就能把線性瞬時混合的多路信號按統(tǒng)計(jì)獨(dú)立性拆回源信號。它押注的是四階累積量高斯噪聲的四階累積量為零所以現(xiàn)場噪聲大一點(diǎn)它反而不慌。適合做課程設(shè)計(jì)、做傳感器陣列預(yù)處理、想給 ICA 找一個確定性基準(zhǔn)的人。下面按能跑通的順序拆原理、代碼結(jié)構(gòu)、參數(shù)、踩坑。2. 四階累積量憑什么能定位那個旋轉(zhuǎn)從白化到聯(lián)合對角化2.1 線性瞬時混合與兩個模糊為什么二階統(tǒng)計(jì)量不夠先把模型釘死后面的代碼才有依據(jù)。JADE 面對的是線性瞬時混合x(t) A s(t) n(t)x 是 n 維觀測s 是 m 維源這里假設(shè) m n源數(shù)靠主成分譜自己判斷A 是 n×m 的未知混合矩陣n(t) 是觀測噪聲。要求只有三條源之間統(tǒng)計(jì)獨(dú)立、源里最多一個高斯分布、A 列滿秩。這三條不滿足后面無論怎么調(diào)參數(shù)都是白費(fèi)。只用二階統(tǒng)計(jì)量能做到什么程度對 x 求協(xié)方差 R E[xx^T] A E[ss^T] A^T特征分解后能得到一個白化矩陣 W使 z Wx 滿足 E[zz^T] I。但白化只消除了二階相關(guān)性剩下的自由度是一個正交旋轉(zhuǎn)z U sU 是任意正交矩陣你怎么轉(zhuǎn)二階統(tǒng)計(jì)量都看不出差別。這就是 ICA 里常說的「白化之后只剩一個旋轉(zhuǎn)」這個旋轉(zhuǎn)必須靠更高階的統(tǒng)計(jì)信息來定JADE 用的是四階累積量。另一個繞不開的點(diǎn)是兩個模糊分離出來的信號幅度不確定源乘任意標(biāo)量混合矩陣除以該標(biāo)量觀測不變順序也不確定源的排列可以任意換。所以判斷分離好壞不能拿分離信號和原始信號逐點(diǎn)相減看誤差必須用帶排列和縮放不變性的指標(biāo)這一點(diǎn)在第 3 章的 PI 函數(shù)里會落實(shí)到代碼。2.2 白化把混合矩陣壓成一個正交旋轉(zhuǎn)白化的實(shí)現(xiàn)有講究不是一句zscore能替代的。標(biāo)準(zhǔn)做法是對協(xié)方差矩陣做特征分解取特征值的負(fù)平方根做縮放。這里有兩個工程細(xì)節(jié)一是協(xié)方差必須強(qiáng)制對稱eig在病態(tài)矩陣上算出來的 R 可能有微小的非對稱導(dǎo)致特征向量出現(xiàn)復(fù)數(shù)二是特征值要設(shè)下限當(dāng)某一路傳感器幾乎沒信號時最小特征值接近零取負(fù)平方根會直接炸掉分離結(jié)果全是 NaN。function [Z, W] whiten(X, r) % X: n x N 觀測矩陣每行一個通道 % r: 保留的主分量個數(shù)缺省等于通道數(shù) n [n, N] size(X); if nargin 2 || isempty(r), r n; end Xc X - mean(X, 2); % 逐通道去均值累積量對均值敏感 R (Xc * Xc) / N; % 協(xié)方差 [E, D] eig((R R) / 2); % 強(qiáng)制對稱避免復(fù)數(shù)特征向量 [d, idx] sort(diag(D), descend); E E(:, idx); d d(1:r); floor_val max(d) * 1e-10; % 特征值下限防止病態(tài)放大 d(d floor_val) floor_val; W diag(d.^(-1/2)) * E; % 白化矩陣 Z W * Xc; % 白化后 E[zz^T] ≈ Ir這個參數(shù)決定保留幾個主分量也就間接決定了后面認(rèn)為有幾個源。常見錯誤是直接取r n明明只有兩路有效信號硬塞三路進(jìn)去多出來的那一維全是噪聲累積量矩陣束會被噪聲主導(dǎo)分離結(jié)果看起來「拆開了但每一路都很臟」。判斷方法在 4.2 節(jié)給看特征值譜有沒有斷崖。白化還有一個容易忽略的前提——各通道的量綱要大致可比。三路傳感器一個是毫伏級、一個是伏級協(xié)方差矩陣的條件數(shù)會非常大白化相當(dāng)于把弱通道放大一千倍同時把它的噪聲也放大一千倍。常見做法是白化前對每個通道做一次標(biāo)準(zhǔn)差歸一化讓協(xié)方差矩陣的對角線都在 1 附近這一步在實(shí)測數(shù)據(jù)上經(jīng)常是決定性的。2.3 四階累積量矩陣束為什么是 Q(M) 而不是 n 的四次方張量四階累積量完整形式是一個 n×n×n×n 的張量元素個數(shù)是 n 的四次方。三通道是 81 個八通道就是 4096 個直接對張量做聯(lián)合對角化計(jì)算量和存儲都不劃算。JADE 的做法是把張量沿著若干方向投影成一組 n×n 的矩陣這叫矩陣束matrix pencil。對白化后的數(shù)據(jù) z定義累積量矩陣Q(M) E[ (z^T M z) · z z^T ] ? 2M ? tr(M)·IM 是任意 n×n 的權(quán)重矩陣。這個式子里最關(guān)鍵的性質(zhì)是如果 z 服從高斯分布Q(M) 恒等于零矩陣。推導(dǎo)不復(fù)雜高斯分布的四階累積量為零而白化后的二階項(xiàng)剛好被 ?2M ? tr(M)I 抵消掉。這意味著觀測噪聲只要是高斯的無論多大方差對 Q(M) 的貢獻(xiàn)都是零JADE 的抗噪能力就來自這里而不是來自什么濾波預(yù)處理。權(quán)重矩陣 M 取多少組理論上取 n2 個線性無關(guān)的基就夠了實(shí)際實(shí)現(xiàn)中利用 Q(M) 的對稱性取 n(n1)/2 個對稱基即可覆蓋。三通道是 6 個矩陣八通道是 36 個規(guī)?？煽?。function Q cum4mat(Z, M) % Z: n x N已經(jīng)白化且去均值 % M: n x n 權(quán)重矩陣 [n, N] size(Z); MZ M * Z; % n x N u sum(Z .* MZ, 1); % 1 x N即每個時刻的 z^T M z Q (Z .* u) * Z / N - 2*M - trace(M)*eye(n);這段代碼里Z .* u是隱式擴(kuò)展把第 j 個時刻的標(biāo)量 u(j) 乘到 z(:,j) 上再乘 Z 求平均得到的就是 E[(z^T M z) z z^T]。老版本 MATLAB 不支持隱式擴(kuò)展要寫成bsxfun(times, Z, u)這個兼容性坑在第 5 章還會提一次。2.4 聯(lián)合對角化的旋轉(zhuǎn)角閉式解現(xiàn)在手上有一組矩陣 Q?…Q_K性質(zhì)是當(dāng)分離矩陣 U 取對了U^T Q_k U 對每個 k 都是對角陣。原因很直接源之間獨(dú)立四階互累積量為零只有自累積量留在對角線上。于是問題變成找一個正交矩陣 U讓這組矩陣同時盡可能對角目標(biāo)函數(shù)寫成非對角元素平方和J(U) Σ_k off(U^T Q_k U)2求解用 Jacobi 旋轉(zhuǎn)每次挑兩個坐標(biāo)軸 (p, q)做一個 2×2 的 Givens 旋轉(zhuǎn)把所有矩陣的 (p,q) 和 (q,p) 元素往零推逐對掃完一遍算一次迭代反復(fù)掃到收斂。關(guān)鍵是這個旋轉(zhuǎn)角有閉式解不用做一維搜索。記 u_k Q_k(p,p) ? Q_k(q,q)v_k Q_k(p,q) Q_k(q,p)累加 C Σu_k2、D Σu_k v_k、E Σv_k2則最優(yōu)旋轉(zhuǎn)角滿足θ 0.25 · atan2(2D, C ? E)這個式子來自「最大化旋轉(zhuǎn)后對角元素平方和」的推導(dǎo)和最小化非對角平方和是等價(jià)的。驗(yàn)證一下單個矩陣 [[1, 0.5],[0.5, 1]]u 0、v 1C 0、D 0、E 1atan2(0, ?1) πθ π/4正好把它轉(zhuǎn)成對角陣 [[0.5, 0],[0, 1.5]]。旋轉(zhuǎn)周期是 π/2θ 落在 (?π/4, π/4] 恰好覆蓋一個周期不會漏解。提示這里的 θ 定義和教科書上某些 Givens 旋轉(zhuǎn)差一個符號取決于旋轉(zhuǎn)矩陣寫成 G^T A G 還是 G A G^T。抄公式的時候必須和代碼里的更新方式對齊否則迭代會發(fā)散而不是收斂。3. 代碼怎么組織五個環(huán)節(jié)的接口與回代順序3.1 主腳本先造一組自己知道答案的混合信號驗(yàn)證任何盲分離程序第一步都是造一個「知道真值」的數(shù)據(jù)集。不要上來就拿實(shí)測數(shù)據(jù)調(diào)那樣分離錯了你根本不知道是算法問題還是數(shù)據(jù)問題。下面的主腳本生成三路特性不同的源純正弦二階統(tǒng)計(jì)量強(qiáng)、四階弱、方波四階很強(qiáng)、調(diào)幅信號時變包絡(luò)考驗(yàn)算法對非平穩(wěn)性的魯棒性再用一個條件數(shù)適中的矩陣混合最后疊一點(diǎn)高斯噪聲。rng(7); % 固定隨機(jī)種子保證結(jié)果可復(fù)現(xiàn) fs 500; T 4; N fs * T; t (0:N-1) / fs; S_true [ sin(2*pi*7*t); % 源1純正弦 sign(sin(2*pi*13*t)); % 源2方波四階累積量強(qiáng) (1 0.6*sin(2*pi*2.5*t)) .* sin(2*pi*41*t) ]; % 源3調(diào)幅 A_true [1.0 0.6 0.3; 0.4 1.0 0.5; 0.2 0.5 1.0]; % 對角線占優(yōu)條件數(shù)約 2數(shù)值友好 X A_true * S_true; X X 0.05 * randn(size(X)); % 高斯觀測噪聲JADE 對它不敏感 [S_est, W, A_est] jade(X, 3); % r3認(rèn)為有三個源 PI perf_index(W, A_true); fprintf(分離矩陣與真實(shí)混合矩陣的 PI %.4f\n, PI);這里的rng(7)不是可有可無的裝飾。四階累積量的估計(jì)方差比二階大得多換一個隨機(jī)種子PI 可能從 0.02 跳到 0.15尤其是樣本量在兩三千點(diǎn)的時候。做對比實(shí)驗(yàn)必須固定種子否則你調(diào)的「參數(shù)優(yōu)化」很可能只是隨機(jī)波動。A_true我特意選成對角占優(yōu)條件數(shù)約 2。如果混合矩陣接近奇異比如兩列幾乎平行那本質(zhì)上就是兩路信號在觀測空間里分不開任何盲分離算法都會翻車這屬于數(shù)據(jù)問題而不是算法問題。3.2 jade 主函數(shù)白化、累積量、回代三步的接口約定主函數(shù)只做三件事白化、構(gòu)造矩陣束、聯(lián)合對角化然后把分離矩陣回代到原始觀測上。接口設(shè)計(jì)上有個細(xì)節(jié)值得說返回分離矩陣 W而不是只返回分離信號 S。因?yàn)樽鲂阅茉u估、做后續(xù)的通道匹配、做在線分離都需要 W 這個線性算子只給 S 的話新來一批數(shù)據(jù)你還得重新跑一遍整個算法。function [S, W, A_est] jade(X, r) % X: n x N 觀測 % r: 源數(shù)目保留的主分量個數(shù) [n, N] size(X); if nargin 2 || isempty(r), r n; end Xc X - mean(X, 2); % 去均值 [Z, Wz] whiten(Xc, r); % 白化Wz 是 r x n Qset cum4_set(Z); % n x n x K 的矩陣束 [U, ~] joint_diag(Qset, 1e-8, 200); % U 是 r x r 正交陣 W U * Wz; % 總分離矩陣 r x n S W * Xc; % 分離信號 r x N A_est pinv(W); % 估計(jì)的混合矩陣 n x r矩陣束的生成函數(shù)把 n(n1)/2 個基矩陣遍歷一遍注意非對角基要對稱化寫成 E(p,q) E(q,p) 1這樣 Q(M) 才是實(shí)對稱矩陣后面的 Jacobi 旋轉(zhuǎn)才能用實(shí)數(shù)的閉式解。如果只寫 E(p,q) 1得到的 Q 是非對稱的聯(lián)合對角化就得處理復(fù)數(shù)代碼復(fù)雜度翻倍不說收斂也慢。function Qset cum4_set(Z) n size(Z, 1); K n * (n 1) / 2; Qset zeros(n, n, K); k 0; for p 1:n for q p:n E zeros(n); if p q E(p,p) 1; % 對角基 else E(p,q) 1; E(q,p) 1; % 對稱非對角基 end k k 1; Qset(:,:,k) cum4mat(Z, E); end end接口一覽函數(shù)輸入輸出說明whitenX (n×N), rZ (r×N), W (r×n)特征值裁剪輸出白化數(shù)據(jù)與白化矩陣cum4matZ (n×N), M (n×n)Q (n×n)單個權(quán)重矩陣對應(yīng)的累積量矩陣cum4_setZ (n×N)Qset (n×n×K)生成完整矩陣束K n(n1)/2joint_diagQset, tol, maxiterU (n×n)正交聯(lián)合對角化返回旋轉(zhuǎn)矩陣perf_indexW, A_true標(biāo)量帶排列與縮放不變性的性能指數(shù)本文還有配套的精品資源點(diǎn)擊獲取