ILRMA盲源分离算法原理详解与MATLAB实现

简介面向音频信号分离与盲源分离研究者的MATLAB实现包聚焦独立低秩矩阵分析ILRMA算法。ILRMA结合独立成分分析与低秩矩阵假设常用于多通道信号分解、声源提取等场景。压缩包共15个文件体积18.06MB以9个.m脚本为主体涵盖whitening预处理、STFT/ISTFT变换、ILRMA主算法及一致ILRMA、ISS变体等另含piano.wav与drums.wav两段测试音频便于直接运行main.m查看分离效果3篇PDF为Kitamura等作者的期刊文献帮助理解算法理论README.md提供使用说明。目前已有163人学习下载。通过研读脚本与对比实验结果可掌握ILRMA从数据预处理、模型迭代到结果解析的完整流程代码目录将主程序与函数模块分离便于按需调用、二次开发并可在示例音频或自有数据上快速验证适合具备MATLAB基础、希望系统学习盲源分离算法的学生与工程师。1. 为什么ILRMA值得在MATLAB里重新写一遍盲源分离BSS里有两派常用的方法独立成分分析ICA强调源信号的统计独立性却对语音、音乐这类有强时频结构的信号视而不见非负矩阵分解NMF擅长抓取频谱的低秩模式但本身不能直接估计解混矩阵。ILRMAIndependent Low-Rank Matrix Analysis用概率模型把这两件事放在同一个目标函数里每个源在时频域变成一个非负低秩张量源与源之间再做独立性约束。这个模型尤其适合双通道或更多通道的语音/音乐分离效果往往好过单独跑IVA或NMF。MATLAB写ILRMA主循环不算难难点在预处理、初始化、更新规则的细节和收敛判断。下面从数学讲到可运行脚本按我实际调试的路径来。2. ILRMA的概率模型先懂低秩假设再写循环2.1 从复时频观测到低秩源模型假设麦克风观测信号经过短时傅里叶变换STFT后得到一个复数张量 X ∈ C^{I×J×T}其中 I 是通道数J 是频率点数T 是帧数。瞬时混合模型在频域写为X(f,t) A(f) S(f,t)A(f) 是 I×I 的复混合矩阵S(f,t) 是源信号的复频谱向量。ILRMA 的出发点是直接对 S 做独立假设不够因为语音和音乐的频谱包络随时间变化有很强相关性。于是它假设每个源 n 的功率谱密度满足r_n(f,t) sum_{k1}^{K} u_{n,k}(f) v_{n,k}(t)这里 u_{n,k}(f) 是非负的谱基v_{n,k}(t) 是对应的时间激活。也就是说第 n 个源在时频平面上的方差可以分解成 K 个秩一矩阵的和。这个低秩假设比“平稳高斯”更接近语音的实际情况也让 NMF 的乘法更新可以直接派上用场。需要特别强调这里的 f 索引的是 STFT 的频率点而不是物理频率。混合矩阵 A(f) 之所以随频率变化是因为相位差随频率变化在近场或混响环境下A(f) 甚至是复数矩阵。ILRMA 对 A(f) 不做过强的结构化约束所以它天然支持多通道并且不要求通道数必须等于源数。不过脚本里一般假设两者相等这样 W(f) 是方阵迭代投影更新会更平稳解混矩阵的可逆性也更容易保持。2.2 NMF式方差先验为什么是K个秩一基如果用最小化负对数似然的视角看ILRMA 在更新 W 时把时频方差 r_n(f,t) 当作已知量源信号被建模为复高斯分布y_n(f,t) ~ N_c(0, r_n(f,t))这里的 y_n(f,t) 是分离后的第 n 个源的频谱。等价地|y_n(f,t)|^2 的期望是 r_n(f,t)。由于 r_n 被分解成非负的 u 和 v整个模型对时间方向要求所有帧共享相同的谱基对频率方向要求所有频点共享相同的时间激活。这比 ICA 里常用的非高斯分布更灵活因为低秩结构覆盖了频谱包络的动态变化。语音和音乐的谐波结构在时频图上恰恰表现为能量集中在少数秩一成分上所以用 K 个基去逼近功率谱比单纯假设源“超高斯”要合理得多。K 的取值直接决定模型容量K 太小r_n 无法刻画频谱细节分离会像滤波而不是选择K 太大NMF 部分会把噪声也建模成“源”分离后每个信道残留很多串扰。经验上语音用 K4 到 8音乐用 K8 到 20。具体怎么定我在第 4 章会给一张参数表。2.3 目标函数和交替更新策略把独立性约束放进同一个目标函数得到如下的代价函数省略常数项L sum_{f1}^{J} [ sum_{t1}^{T} sum_{n1}^{I} ( |y_n(f,t)|^2 / r_n(f,t) log r_n(f,t) ) - 2T log |det W(f)| ]其中第一项是数据拟合项第二项是 W(f) 的雅可比项它保证了 y 和 x 之间的概率密度变换关系。优化上面这个 L 的思路就是交替固定当前 W更新 U 和 V。这本质上是 Itakura-Saito 散度下的 NMF 分解乘法更新规则保证非负性。固定 U、V更新 W。这一步对每个频率 f 独立进行用迭代投影iterative projection保持 W(f) 的可逆性。这两种更新都保证目标函数非增。在 MATLAB 里每次迭代可以先做一次 NMF 更新再做一次迭代投影也可以反过来差别不大但通常 W 更新更耗时所以我一般把 NMF 更新放在内层。ILRMA 与独立向量分析IVA的区别在于IVA 假定 y_n 的时间方差是随机的标量而 ILRMA 用 NMF 显式建模这个方差。这样得到的分离矩阵在结构上更稳定短帧数据下不容易过拟合。方法源模型是否显式解混矩阵适合场景ICA非高斯 i.i.d.是瞬时混合、超高斯信号IVA各源多变量独立是多帧语音NMF低秩谱模型否需额外估计矩阵单通道源分离ILRMA低秩谱 独立源是多通道语音/音乐分离从表中能看出ILRMA 的目标不是替代所有方法而是把 NMF 的低秩表达能力和 IVA 的多通道独立性约束拼到一个框架里。后面所有 MATLAB 代码都是围绕这套交替更新展开的。3. MATLAB实现ILRMA主循环预处理、初始化、迭代更新3.1 用spectrogram把多通道音频变成复数时频张量在 MATLAB 里做 STFT 最直接的是用 spectrogram 函数但 spectrogram 默认返回单边频谱而且一次只处理一个通道。我一般先写一个小的转换函数把多通道音频读成统一的复数张量 X 三维频率 × 帧 × 通道。注意如果你想做可逆分离还需要保存窗函数、跳数、FFT 点数等参数。function X audios_to_stft(x, win, hop, nfft) [nCh, nSamp] size(x); if nSamp hop, error(音频太短); end x x ./ (max(abs(x(:))) eps); % 简单幅值归一化 nFrames floor((nSamp - nfft) / hop) 1; nFreq nfft / 2 1; % 单边频谱 X zeros(nFreq, nFrames, nCh); % 复数张量 window hamming(nfft, periodic); for c 1:nCh for t 1:nFrames seg x(c, (t-1)*hop (1:nfft)) .* window; X(:, t, c) fft(seg, nfft); % 复数谱 end end X X(1:nFreq, :, :); % 保留正频率 end这段代码里nFreq是单边谱的维度hop不要取得比窗长小太多否则帧之间重复计算会让收敛变慢。我在 16kHz 下常用 1024 点窗、256 点 hop也就是 75% 重叠。注意输入 x 必须是 通道×采样 的矩阵如果是行向量nCh1但 ILRMA 至少要 2 个通道所以先检查维度。此外FFT 结果去掉负频率段因为正频率段已经包含完整幅度和相位信息逆变换时再补回共轭对称部分。3.2 初始化解混矩阵和低秩因子W 的初始化不能全零也不能随机复数因为迭代投影需要 W(f) 可逆。常见做法是每个频率点都用单位矩阵然后加很小的噪声或者用随机正交矩阵。NMF 的 U、V 初始化成正的随机数并做一次简单缩放让 r_n 的初始均值接近 1会减少前几次迭代的震荡。function [W, U, V] init_ilrma(nFreq, nFrames, nSrc, K) W zeros(nSrc, nSrc, nFreq); for f 1:nFreq A randn(nSrc, nSrc) 1i*randn(nSrc, nSrc); [Q, ~] qr(A); % 正交基 W(:,:,f) Q; % 保持可逆 end U cell(nSrc, 1); V cell(nSrc, 1); for n 1:nSrc U{n} rand(nFreq, K) 0.1; V{n} rand(K, nFrames) 0.1; scale sqrt(mean(U{n}(:)) * mean(V{n}(:))); U{n} U{n} / scale; V{n} V{n} / scale; end end这里的qr确保 W 的每一列正交scale让 U、V 的乘积不至于一开始就很大。初始化 W 用复数随机矩阵的原因是混合矩阵的相位随频率变化实数初始化会让某些频率点接近奇异影响收敛。你可以保持 W 为单位矩阵但那样分离能力会弱一些因为初始点离最优解更远。3.3 交替更新乘法规则与迭代投影NMF 更新部分对每个源独立。对于源 n先计算当前功率谱 P_n |y_n|^2再计算由 U_n、V_n 重构出的方差 R_n U_n * V_n。乘法更新规则如下U_n ← U_n .* sqrt( (P_n ./ R_n.^2) * V_n^T ./ ( (1./R_n) * V_n^T ) ) V_n ← V_n .* sqrt( U_n^T * (P_n ./ R_n.^2) ./ ( U_n^T * (1./R_n) ) )这些运算在 MATLAB 里要用./和.*小心处理分母加 eps 防止除零。W 的更新用迭代投影下面是单个频率 f 的处理逻辑function Wf ip_update(Xf, Rf, Wf) [nSrc, T] size(Xf); for n 1:nSrc Vn (Xf ./ Rf(n,:)) * Xf / T; % 加权协方差 w (Wf * Vn) \ eye(nSrc, n); w w / sqrt(real(w * Vn * w) 1e-12); Wf(n,:) w.; end end这里Xf ./ Rf(n,:)表示把矩阵 Xf 的每一列除以标量 r_n(f,t)MATLAB 的广播机制会自动完成。注意eye(nSrc, n)是第 n 列单位向量也就是 e_n。迭代投影的本质是把当前 W 的一行替换为在加权协方差意义下的最小方差响应然后归一化。这个更新既保持了 W 的可逆性又让分离后的 y_n 的功率和目标方差 r_n 对齐。3.4 主循环的完整伪代码把上面的片段串起来主循环结构为for iter 1:n_iter % 1) 分离当前源 Y zeros(J, T, nSrc); for f 1:J Y(f,:,:) (W(:,:,f) * squeeze(X(f,:,:))).; end % 2) NMF更新 for n 1:nSrc Yn squeeze(Y(:,:,n)); Pn abs(Yn).^2; Rn U{n} * V{n}; U{n} U{n} .* sqrt( (Pn ./ (Rn.^2 reg)) * V{n} ./ ((1 ./ (Rn reg)) * V{n}) ); V{n} V{n} .* sqrt( U{n} * (Pn ./ (Rn.^2 reg)) ./ (U{n} * (1 ./ (Rn reg))) ); end % 3) 迭代投影更新W for f 1:J Xf squeeze(X(f,:,:)); Yf W(:,:,f) * Xf; Rf zeros(nSrc, T); for n 1:nSrc Rn U{n} * V{n}; Rf(n,:) Rn(f,:); end W(:,:,f) ip_update(Xf, Rf, W(:,:,f)); end end这段代码已经接近可以运行但缺少边界条件和对 U、V 的尺度保护。我在第 5 章给一个更完整的版本并说明怎么验证分离结果。注意每次迭代结束后最好把 U、V 稍微做一次归一化否则 NMF 的尺度会漂移影响 W 更新的数值稳定性。由于 NMF 更新只做一次所以每一步其实是坐标下降的一小步而不是内层循环完全收敛这在 ILRMA 里是正常的外层迭代会逐步逼近最优。4. ILRMA的参数设置与收敛性诊断别等跑完再后悔4.1 关键参数表K、窗长、迭代次数怎么定参数推荐范围对结果的影响我的常用值低秩基个数 K2~20K 过小不能刻画谐波结构过大会把噪声建模成源语音 6音乐 12STFT 窗长 nfft256~2048窗长决定频率分辨率短窗时间分辨率高但频谱粗糙51216kHz帧移 hopnfft/4 ~ nfft/2帧移小则帧间冗余多收敛更平稳但更慢nfft/4最大迭代次数30~200ILRMA 收敛较慢但 30 轮后分离比提升不明显80W 初始化单位矩阵/随机正交随机正交可避免奇异但重复实验差异大随机正交U/V 初始化正随机方差量级最好接近 1见 3.2 节代码这里的 K 是最需要调的参数。如果你发现分离出的源里有明显“音乐噪声”或嗡嗡声往往是 K 太大把噪声低秩化了如果分离不彻底源里还混着另一路声音则 K 太小模型无法表达源内的时间变化。另外窗长选择要匹配采样率16kHz 下 512 点窗对应 31.25ms刚好能分辨语速较快的辅音8kHz 下用 256 点窗比较合理。4.2 收敛性监测算目标函数和分离指标不要只靠听结果判断是否收敛我习惯每 5 轮算一次目标函数 L看它是不是在单调下降。实现起来不高只需要在更新前存一份 W 和 U/V 的旧值更新后按公式计算。也可以用一个更简单的代理指标所有源之间在时域上的相关系数绝对值。如果两路输出 y_1 和 y_2 的相关系数接近于 0说明独立性已建立。function cost ilrma_cost(Y, U, V, W) [J, T, nSrc] size(Y); cost 0; for f 1:J Yf squeeze(Y(f,:,:)); % nSrc x T detW abs(det(W(:,:,f))); for n 1:nSrc Rn U{n} * V{n}; r Rn(f,:); % 1 x T p abs(Yf(n,:)).^2; cost cost sum(p ./ r log(r 1e-12)); end cost cost - 2 * T * log(detW 1e-12); end end这段代码里log(r)要加 1e-12 防止 log(0)。cost是一个实数标量每轮更新完后重新计算如果出现连续三次不降反升就要考虑是不是步长或者初始化出了问题。注意这里Y是用当前 W 和 X 计算得到的如果你在更新过程中复用了旧 Y算出来的 cost 会不准确。我一般在主循环里每隔mod(iter,5)0调用一次这个函数把结果打印出来。4.3 常见发散的排查顺序我调试 ILRMA 时遇到最多的问题就是 NaN 和无穷大。出现 NaN 的第一反应不是加 eps而是先看 R_n 是不是有零元素。R_n 是由 U_n 和 V_n 相乘得到的只要 U、V 中有元素在乘法更新中变成 0后续除法就会爆炸。所以我的排查顺序是检查输入 X 是否包含 NaN 或直流偏置。STFT 前先x x - mean(x)。检查初始化后的 R_n 是否小于 1e-10。是的话说明 scale 没起效。在乘法更新分母处统一 1e-12而不是用epseps 在双精度下是 2.2e-16太小。更新 W 后检查det(W(:,:,f))是否接近 0。如果接近 0说明该频率点上的混合矩阵估计不可靠可以把这个频段的 W 重置为单位阵让迭代再继续走。如果只有个别频率发散说明低秩基个数 K 对该频段不匹配试着增大或减小 K 再看。这套排查顺序能解决 90% 的 NaN 问题。注意迭代投影自身不发散发散的根因几乎都在 NMF 更新产生的 0 元素上。另一个常见问题是如果你看到 cost 在下降但分离结果很差那大概率是通道顺序发生置换或者频谱排列有问题需要检查 STFT 和 ISTFT 的对称性。5. 一个可运行的完整ILRMA分离脚本5.1 整体脚本与输入输出约定下面这个脚本demo_ilrma.m接收一个双通道 WAV 文件作为混合信号输出两个分离后的 WAV 文件。它把前面几节的函数串在一起并加了必要的保护。为了可以直接复制运行我把audios_to_stft、init_ilrma、ip_update、stft_back都放在同一个函数文件里这样不依赖额外附件。核心参数放在文件头部修改起来方便。5.2 脚本代码function demo_ilrma(infile, outfile1, outfile2) [x, fs] audioread(infile); if size(x,2) 2, error(需要至少两个麦克风通道); end x x(:,1:2); nSrc 2; K 6; nfft 512; hop 128; n_iter 80; reg 1e-12; X audios_to_stft(x, hamming(nfft,periodic), hop, nfft); [J, T, ~] size(X); [W, U, V] init_ilrma(J, T, nSrc, K); for iter 1:n_iter Y zeros(J, T, nSrc); for f 1:J Y(f,:,:) (W(:,:,f) * squeeze(X(f,:,:))).; end for n 1:nSrc Yn squeeze(Y(:,:,n)); Pn abs(Yn).^2; Rn U{n} * V{n}; U{n} U{n} .* sqrt( (Pn ./ (Rn.^2 reg)) * V{n} ./ ((1 ./ (Rn reg)) * V{n}) ); V{n} V{n} .* sqrt( U{n} * (Pn ./ (Rn.^2 reg)) ./ (U{n} * (1 ./ (Rn reg))) ); end for f 1:J Xf squeeze(X(f,:,:)); Yf W(:,:,f) * Xf; Rf zeros(nSrc, T); for n 1:nSrc Rn U{n} * V{n}; Rf(n,:) Rn(f,:); end W(:,:,f) ip_update(Xf, Rf, W(:,:,f)); end end y_est zeros(nSrc, size(x,2)); for n 1:nSrc Yn squeeze(Y(:,:,n)); y_est(n,:) stft_back(Yn, hamming(nfft,periodic), hop, nfft, size(x,2)); end audiowrite(outfile1, y_est(1,:), fs); audiowrite(outfile2, y_est(2,:), fs); end function X audios_to_stft(x, win, hop, nfft) [nCh, nSamp] size(x); nFrames floor((nSamp - nfft) / hop) 1; nFreq nfft / 2 1; X zeros(nFreq, nFrames, nCh); for c 1:nCh for t 1:nFrames seg x(c, (t-1)*hop (1:nfft)) .* win(:); X(:, t, c) fft(seg, nfft); end end X X(1:nFreq, :, :); end function [W, U, V] init_ilrma(nFreq, nFrames, nSrc, K) W zeros(nSrc, nSrc, nFreq); for f 1:nFreq A randn(nSrc, nSrc) 1i*randn(nSrc, nSrc); [Q, ~] qr(A); W(:,:,f) Q; end U cell(nSrc, 1); V cell(nSrc, 1); for n 1:nSrc U{n} rand(nFreq, K) 0.1; V{n} rand(K, nFrames) 0.1; scale sqrt(mean(U{n}(:)) * mean(V{n}(:))); U{n} U{n} / scale; V{n} V{n} / scale; end end function Wf ip_update(Xf, Rf, Wf) [nSrc, T] size(Xf); for n 1:nSrc Vn (Xf ./ Rf(n,:)) * Xf / T; w (Wf * Vn) \ eye(nSrc, n); w w / sqrt(real(w * Vn * w) 1e-12); Wf(n,:) w.; end end function y stft_back(Y, win, hop, nfft, nSamp) [J, T] size(Y); y zeros(1, (T-1)*hop nfft); win win(:).; for t 1:T spec Y(:,t); spec [spec; conj(spec(J-1:-1:2))]; seg real(ifft(spec, nfft)) .* win; idx (t-1)*hop (1:nfft); y(idx) y(idx) seg; end y y(1:nSamp); end5.3 运行方式和结果验证运行前用audiowrite准备一个混合文件或者直接用mix cat(2, speech1, speech2)在内存里混合。调用时只需要一句demo_ilrma(mix.wav, sep1.wav, sep2.wav);然后用audioread检查输出。这里要说一下ILRMA 的输出会有幅度和顺序的置换不确定性这是盲源分离的固有问题不能要求输出顺序和输入一致。验证分离效果时我会先做幅度归一化再计算每路输出与真实源的相关性取绝对值最大的一对作为匹配。如果你没有真实源就听分离结果里有没有明显的串扰声。另外这段脚本没有做帧间的重叠相加归一化所以在窗函数边界处会有轻微的幅值抖动。对分离任务来说听感差别不大但如果要拿去做后续特征提取建议在stft_back里加一个正常的窗和重合补偿或者直接用istft函数。6. 把ILRMA脚本改得快一点的三个工程技巧6.1 对频率轴用parfor并行ILRMA 主循环里 W 的更新是对每个频率点独立执行的天然适合并行。在 MATLAB 里把第 5 章脚本中的频率循环替换成parfor只需要注意循环里不能修改共享变量。常见做法是先准备一个临时数组Wf_list每个频率计算后写入循环结束后再拼回 W。parfor f 1:J Xf squeeze(X(f,:,:)); Yf W(:,:,f) * Xf; Rf zeros(nSrc, T); for n 1:nSrc Rn U{n} * V{n}; Rf(n,:) Rn(f,:); end W_f ip_update(Xf, Rf, W(:,:,f)); Wf_list{f} W_f; end for f 1:J W(:,:,f) Wf_list{f}; end注意parfor的循环体里不能直接用squeeze(X(f,:,:))来更新 W否则 MATLAB 会让你报错。把 W 的赋值放到循环外合并这是最稳的写法。6.2 早期用短窗、后期用长窗窗长决定了频率分辨率。短窗 STFT 的计算量小且时间帧数多NMF 的时间基能更快更新。我经常在迭代前 20 轮用 256 点窗分离出一个粗糙结果然后再切换到 512 点窗继续精调。两种窗对应的频率轴长度不同需要把 W 在频率上插值U 的谱基也要插值。这个技巧在 16kHz 下很实用能省下三分之一的时间而且最终分离质量不会明显下降。6.3 对U和V做快速重新初始化如果发现分离结果陷入局部最优比如两个源各分到一半没必要重新跑全部迭代。常见做法是把 U 和 V 重新随机初始化但保留已经收敛的 W。这样做能让 NMF 部分跳出局部极小而 W 的解混方向大致保持。实现上只需要在当前的 U、V 上乘一个随机的波动系数然后重新运行主循环 20 轮。我会在调试阶段把这个过程封装成一个函数方便反复试。在 16kHz 下我会先把窗长降到 256 跑 20 轮粗分离再用 512 跑 80 轮精调整体收敛速度大约能快一倍。本文还有配套的精品资源点击获取