企业官网建设流程全解析

简介：面向音频信号分离与盲源分离研究者的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）擅长抓取频谱的低秩模式，但本身不能直接估计解混矩阵。ILRMA（Independent 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_{k=1}^{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 部分会把噪声也建模成“源”，分离后每个信道残留很多串扰。经验上，语音用 K=4 到 8，音乐用 K=8 到 20。具体怎么定，我在第 4 章会给一张参数表。

2.3 目标函数和交替更新策略

把独立性约束放进同一个目标函数，得到如下的代价函数（省略常数项）：

L = sum_{f=1}^{J} [ sum_{t=1}^{T} sum_{n=1}^{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 显式建模这个方差。这样得到的分离矩阵在结构上更稳定，短帧数据下不容易过拟合。

从表中能看出，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 必须是 通道×采样 的矩阵，如果是行向量，nCh=1，但 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、窗长、迭代次数怎么定

这里的 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，后续除法就会爆炸。所以我的排查顺序是：

1. 检查输入 X 是否包含 NaN 或直流偏置。STFT 前先x = x - mean(x)。
2. 检查初始化后的 R_n 是否小于 1e-10。是的话说明 scale 没起效。
3. 在乘法更新分母处统一+ 1e-12，而不是用eps（eps 在双精度下是 2.2e-16，太小）。
4. 更新 W 后检查det(W(:,:,f))是否接近 0。如果接近 0，说明该频率点上的混合矩阵估计不可靠，可以把这个频段的 W 重置为单位阵，让迭代再继续走。
5. 如果只有个别频率发散，说明低秩基个数 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); end

5.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 轮精调，整体收敛速度大约能快一倍。

本文还有配套的精品资源，点击获取