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