nara wpe git地址
博客中demo代码下载
参考论文

NARA - WPE: A Python Package for Weighted Prediction Error Dereverberation in Numpy and Tensorflow for Online and Offline ProcessingNAR

记录一下nara wpe去混响的流程，将nara库的主要代码抽出进行分析。首先分析离线去混响的方法，离线去混响是在已经得到全部数据的情况下对整体音频进行去混响。之后分析了在线去混响的流程，其中单独复习了一下卡尔曼滤波的基本流程，在线去混响摒弃了离线去混响的一些复杂计算，使用卡尔曼滤波预测去混响的滤波器，同时在迭代计算滤波器的中间过程中完成了去混响操作。

1.WPE方法去混响的基本流程

1.1.基本流程

估计混响信号的尾部。当前信号减去混响尾部，获得“早期”语音，也就是我们需要的去混响后的语音。

$${\hat{x}}_{t,f,d}^{(early)}=y_{t,f,d}-\sum _{\tau =\Delta }^{\Delta +K-1}\sum _{d^{\prime }}{g}_{\tau ,f,d,d^{\prime }}^{\ast }{y}_{t-\tau ,f,d^{\prime }}$$

这里是时域的表示，主要用于理解需要处理的数据，但是公式中包含了频率f，需要使用不同频率的特征。实际操作时使用stft先变换，这样同时包含时间和频率信息，处理完毕后通过istft反变换到时域。

下面是各个符号的含义：

t时刻，频率为f，通道d上的早期反射信号x，也就是我们最终希望得到的干净的时域信号。

$${\hat{x}}_{t,f,d}^{(early)}$$

t时刻，频率为f，通道d上的带混响的时域信号y。

$${\mathbf{y}}_{\mathrm{t},f,d}$$

t-τ时刻，频率为f，通道间d’上的带混响的时域信号y。实际操作中，我们会使用多个τ的数据点来估计混响信号的尾部。d’代表了所有通道的影响。这里需要的维度是(时间, 频率, 通道)。后续代码处理是将y通过stft变换后经过延迟，变为为y_tiled后，组成一个tensor( 频率, Kx通道数, 通道上每个时间的系数)，我的理解是 Kx通道数 这个维度包含了不同延迟时间内通道间的关系d’，另外两个维度则体现了频率f和时间。

$${\mathbf{y}}_{\mathrm{t}-\tau ,f,d^{\prime }}$$

下面是对两边的两段时域信号做stft运算，得到当前带混响的信号和估计的混响信号的频域表示。这样两者相减，就可以得到早期反射信号X。最后再通过istft反变换到时域。

$${\hat{\mathbf{x}}}_{t,f}^{(\text{early})}={\mathbf{y}}_{t,f}-{\mathbf{G}}_f^{\mathcal{H}}{\tilde{\mathbf{y}}}_{t-\Delta ,f}$$

首先要看一下如何构造这个y_tiled(y波浪)。y_tiled的基本组成上面已经交代，其中包含了频率信息、不同tap（延迟）内通道间信息（K个不同时间）、每个通道的系数。

下面要做的是求Gf中的系数，这个系数没有办法直接球出来。Gf将当前y_tilde转换为估计的混响尾部信号，Gf的形状为(频点数, D*tap, D)。这个Gf目前没有办法直接求解，只能在单通道直达信号是零均值复高斯分布的假设下，使模型似然最大化。

单通道分布：

$$p(x_{t,f,d}^{(\text{early})};{\lambda }_{t,f})=\mathcal{N}(x_{t,f,d}^{(\text{early})};0,{\lambda }_{t,f})$$

等式右边代表0均值，方差λt,f的复高斯分布概率密度函数，左边代表似然函数在给定方差λt,f的情况下观测到早期语言估计x的概率。注意，这里是信号频率的系数均值假设是0，到这里处理的全部一段音频是stft后的结果。也就是是说可以直接求同一频率下同一时间，多个通道的频点上所有系数的平均能量，就是这里所讲的方差λt,f了。

多通道分布:

$$p({\mathbf{x}}_{t,f}^{(\text{early})};{\lambda }_{t,f})=\mathcal{N}({\mathbf{x}}_{t,f}^{(\text{early})};0,{\lambda }_{t,f}\mathbf{I})$$

到这里知道了方差的求法，根据方差来估计早期语音出现的概率。下面就需要看看具体如何求解Gf，求出Gf就知道了混响尾部与tap次观测值之间的关系了，通过观测信号和Gf，得到混响尾部的估计。当前带混响的信号减去这个估计得到干净信号。

2.离线迭代方法

求方差的方法上一小节已经给出，实际就是求同一频率下同一时间，多个通道的频点系数的平均能量值。下面给出公式：

$${\lambda }_{t,f}=\frac{1}{(\delta +1+\delta )D}\sum _{\tau =t-\delta }^{t+\delta }\sum _d|{\hat{x}}_{\tau ,f,d}^{(\text{early})}{|}^2$$

下面就是求Gf的关键部分。

$${\mathbf{R}}_f=\sum _t\frac{{\tilde{\mathbf{y}}}_{t-\Delta ,f}{\tilde{\mathbf{y}}}_{t-\Delta ,f}^{\mathrm{H}}}{{\lambda }_{t,f}}\in {\mathbb{C}}^{DK\times DK}$$

$${\mathbf{P}}_f=\sum _t\frac{{\tilde{\mathbf{y}}}_{t-\Delta ,f}{\mathbf{y}}_{t,f}^{\mathrm{H}}}{{\lambda }_{t,f}}\in {\mathbb{C}}^{DK\times D}$$

$${\mathbf{G}}_f={\mathbf{R}}_f^{-1}{\mathbf{P}}_f\in {\mathbb{C}}^{DK\times D}$$

没有查到明确的Gf的求解方法的解释，下面是根据大模型给出的结果结合自己的理解给出解释

为什么使用观测信号自相关矩阵Rf？

自相关反映了观测信号在不同时间点上的相似程度。在混响环境下，观测信号包含直达声和各种反射声，不同时刻观测信号之间存在一定关联。可以体现信号在不同通道间、不同延迟间的相关性。

同一频率，同一时间的结果除以多mic间方差，如果方差小说明多个mic收到的数据相同，则使结果得到增强，否则削弱结果。

为什么使用观测信号和当前信号的互相关矩阵Pf？

互相关体现了观测信号和当前信号的关联性，用于挖掘当前信号和延迟信号的间的联系。观测到的延迟信号对当前信号的影响。

最终滤波器系数Gf的确定：

Pf是延迟信号与当前信号的关系，Rf是延迟信号自己的自相关，Gf对Rf做线性组合得到Pf，从而体现出延迟信号与观测信号的关系。这样，求解RfxGf=Pf的方程，就能得到Gf。另外这种通过相关获取当前信号和观测信号方法，是否也可以推广到其他场合使用？

最后求解过程可以使用下面的流程，Y是多路mic采集的音频数据stft的结果，将X经过istft后就可以得到去混响后的音频。

这里用相对nara中实现更直观的方法演示生成$\tilde{\mathbf{y}}$

def build_Y_tilde(Y, taps, delay):# stft 有效点数S = Y.shape[:-2]# mic数D = Y.shape[-2]# time步数T = Y.shape[-1]def build_pad(i):if len(S) == 0:pad = np.zeros((D, delay+i))else:pad = np.zeros(S + (D, delay+i))return paddef build_zero_Y_tilde_part(n_height):# Y_tilde tensor底部补0 代表给定taps长度超出了当前序列T得长度 多与数据补0if len(S) == 0:Y_tilde_part = np.zeros((D*n_height, T))else:Y_tilde_part = np.zeros(S + (D*n_height, T))return Y_tilde_part# 构造最终Y_tilde (..., D*tap, T)n_height = 0if T - delay < taps:# 最后Y_tilde的部分需要补0n_height = taps - (T-delay)lst = []# 构造前几个tap数据 for i in range(taps - n_height):# 第i个tap中剩余bufferremain = Y[..., :, :T-delay-i]# 整个tensor左半部分补0pad = build_pad(i)Y_tilde_part = np.concatenate((pad, remain), axis=-1)lst.append(Y_tilde_part)if n_height > 0:# 再往后后续数据还没有叠加到混响 后几个tap补0Y_tilde_part = build_zero_Y_tilde_part(n_height)lst.append(Y_tilde_part)# 所有tensor竖着摞起来ret = np.concatenate(lst, axis=-2)return ret

离线去混响流程

def slow_wpe(Y, iters=5, taps=10, delay=3):X = Y# Y_tilde代表了Y数据随着时间推移对当前Y的影响 # S=1 D=2(mic数) T=6 delay=3 tap=2(2组时间数据) # 第0组时间数据# mic0 [0, 0, 0, 3, 2, 1]# min1 [0, 0，0, 6, 5, 4]# 第1组时间数据 延迟一个时间单位后 先发生的数据被"挤"出时间窗口# mic0 [0, 0, 0, 0, 3, 2]# min1 [0, 0，0, 0, 6, 5]# Y_tilde如下# [0, 0, 0, 3, 2, 1]# [0, 0，0, 6, 5, 4]# [0, 0, 0, 0, 3, 2]# [0, 0，0, 0, 6, 5]# Y_tilde 形状  (频点数, D*tap, T)# D*tap 看作是tap个时间点内 每个时间点下D个mic 每个mic当前频率的系数# T 看作处理信号维持的时间Y_tilde = build_Y_tilde(Y, taps, delay)for i in range(iters):# 更新归一化能量分母inverse_power = build_inverse_power(X)   # 重新归一化Y_tildeY_tilde_inverse_power = Y_tilde * inverse_power[..., None, :]# 延迟信号的自相关  R(频点数, D*tap, D*tap)R = np.matmul(Y_tilde_inverse_power, hermite(Y_tilde))# 延迟信号与当前信号的互相关 (频点数, D*tap, T)*(频点数, T, D) = P(频点数, D*tap, D)P = np.matmul(Y_tilde_inverse_power, hermite(Y))# R*G=P G是延迟信号和当前信号的关系 G(频点数, D*tap, D)G = _stable_solve(R, P)# 当前信号减去回来的延迟信号 (频点数, D, D*tap)*(频点数, D*tap, T) = X(D, T)X = Y - np.matmul(hermite(G), Y_tilde)return X

3.在线求法

在线处理需要不断处理新来的视频帧，无法像离线一样一次获取到最终的${\mathbf{R}}_f$和${\mathbf{P}}_f$后迭代多次求G。针对这种情况，nara中给出的方法是对于先给一个初始矩阵G，当前时间的的R，P使用上个时间的R和P以及当前预测值共同决定。

$${\hat{\mathbf{x}}}_{t,f}^{(\text{early})}={\mathbf{y}}_{t,f}-{\mathbf{G}}_{t-1,f}^{\mathrm{H}}{\tilde{\mathbf{y}}}_{t-\Delta ,f}$$

$${\mathbf{R}}_{t,f}=\alpha {\mathbf{R}}_{t-1,f}+\frac{{\tilde{\mathbf{y}}}_{t-\Delta ,f}{\tilde{\mathbf{y}}}_{t-\Delta ,f}^{\mathrm{H}}}{{\lambda }_{t,f}}$$

$${\mathbf{P}}_{t,f}=\alpha {\mathbf{P}}_{t-1,f}+\frac{{\tilde{\mathbf{y}}}_{t-\Delta ,f}{\mathbf{y}}_{t,f}^{\mathrm{H}}}{{\lambda }_{t,f}}$$

通过迭代的方式获取到${\mathbf{R}}_f$和${\mathbf{P}}_f$后，并没有直接通过求解线性方程组的方式解，而是使用了卡尔曼方法。下面先回顾一下卡尔曼方法的原理和实现，然后再看如何使用卡尔曼方法直接通过${\mathbf{R}}_{t,f}$递推得到${\mathbf{G}}_{t,f}$。

3.1.回顾卡尔曼方法

卡尔曼滤波是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。其核心思想是结合系统的动态模型（预测部分）和实际观测数据（校正部分），以递归的方式不断更新对系统状态的估计，使估计值尽可能接近真实值。具体主要分为三个模块，分别为状态方程、观测方程、预测与校正。

状态方程：描述系统状态随时间的演变规律，也就是状态转移矩阵。

思考1. 状态转移矩阵${\mathbf{A}}_k$是我们认为的当前系统的变化规律，如果一切随我们所愿，A完全准确，那么应该可以直接从${\mathbf{A}}_k$得到下一个时刻的结果。但是由于传感器或者A本身的误差，导致实际情况中下一个时刻输出不一定如我们所愿，这些误差就是下面的${\mathbf{w}}_k$。卡尔曼就是利用${\mathbf{w}}_k$的统计特性，比如下文的协方差矩阵$Q_k$，得到一个相对准确的预测值。

$${\mathbf{x}}_k={\mathbf{A}}_k{\mathbf{x}}_{k-1}+{\mathbf{B}}_k{\mathbf{u}}_k+{\mathbf{w}}_k$$

这里${\mathbf{A}}_k$是状态方程，${\mathbf{B}}_k$是外部控制，可以不用。${\mathbf{w}}_k$是估计噪声，形状应该和$\mathbf{x}$一样是一个向量，每个维度上均值为0。${\mathbf{w}}_k$无法直接计算，是通过协方差矩阵Q来体现的，而Q也无法直接计算，是一个经验值或者通过物理模型统计的值。

思考2. ${\mathbf{w}}_k$在这里用于体现有误差，真实计算过程没有这一项，而是用统计规律构建协方差矩阵$Q_k$来体现误差。

观测方程：建立系统状态与实际观测之间的联系。

$${\mathbf{z}}_k={\mathbf{H}}_k{\mathbf{x}}_k+{\mathbf{v}}_k$$

${\mathbf{z}}_k$是观测向量，是实际的测量结果。${\mathbf{H}}_k$是观测矩阵，述了系统状态如何映射到观测空间，给定输入分离对观测结果的贡献。${\mathbf{v}}_k$是观测噪声向量。

预测与校正：先根据上一时刻状态预测当前状态，再利用当前观测值对预测进行校正，得到更准确的状态估计。

状态预测：根据上一个状态加当前控制输入，得到预测的${\hat{x}}_{k|k-1}$。

$${\hat{x}}_{k|k-1}=A_k{\hat{x}}_{k-1|k-1}+B_k{u}_k$$

协方差预测：$P_k$代表状态估计协方差矩阵，是当前估计值中每个分量自己的方差和分量间的协方差，表征了状态估计的分量间的关系。$Q_k$是过程噪声${\mathbf{w}}_k$的协方差矩阵。

思考3. 为什么不使用当前观测数据直接计算协昌差矩阵$Q_k$而是要加入$A_k$？这里体现了噪声协方差矩阵经过$A_k$后的噪声协方差矩阵。

思考4. 如何确定P的初始值和Q的值？P认为是系统的不确定性，正常初始值是一个每个值都很小的对角矩阵，认为给出的开始状态方程是准确的，也就是开始认为估计值是准的。Q是根据物理规律或者经验值设定。

思考5. 根据思考3、思考4可以看出，这里对表征预测量的分量间关系的协方差矩阵P，是通过经验设定加后续迭代产生的，整个过程没有直接通过统计数据计算协方差矩阵，而是直接递推得到。

思考6. $P_{k-1|k-1}$经过状态矩阵$A_k$后为什么会变为$A_k{P}_{k-1|k-1}{A}_k^{\top }$？$P_{k-1|k-1}$代表了k-1时刻状态间的不确定关系，$A_k$对状态变量坐线性变换后，状态变量间的关系也会发生变化。状态量经过$A_k$后，状态量间的协方差矩阵就是求$\text{cov}(A\mathbf{x},A\mathbf{x})$。向量$A\mathbf{x}$和 $B\mathbf{x}$间协方差遵循 $\text{cov}(A\mathbf{x},B\mathbf{x})=A\text{cov}(\mathbf{x},\mathbf{x})B^{\top }$，这里就是$A_k{P}_{k-1|k-1}{A}_k^{\top }$的来源。

思考7. 到这一步，$A_k$和$Q_k$都是确定的，对于下一轮迭代来说，下面需要得到一个新的$P_{k|k}$，这个值需要和当前观测值有关，否则的话P就成为一个完全确定的序列了。

$$P_{k|k-1}=A_k{P}_{k-1|k-1}{A}_k^{\top }+Q_k$$

校正：

卡尔曼增益计算：

$$K_k=P_{k|k-1}{\mathbf{H}}_k^{\top }{\left({\mathbf{H}}_k{P}_{k|k-1}{\mathbf{H}}_k^{\top }+{\mathbf{R}}_k\right)}^{-1}$$

卡尔曼增益用于计算观测值所占的权重。

思考8. 这里出现观测噪声$v_k$的协方差矩阵${\mathbf{R}}_k$。可以通过经验或者直接统计传感器数据计算出这个值。区别于估计值的协方差矩阵只能使用经验值或者分析物理规律来得到。

思考9. 卡尔曼增益计算公式是怎么来的？推导流程比较复杂，后续补充。

状态更新

$${\hat{\mathbf{x}}}_{k|k}={\hat{\mathbf{x}}}_{k|k-1}+K_k\left({\mathbf{z}}_k-{\mathbf{H}}_k{\hat{\mathbf{x}}}_{k|k-1}\right)$$

本次的预测值 = 上一次的估计值 + K(本次观测值 - 上一次的估计值输入观测)。可以看到卡尔曼最终的预测值是由

上一次的预测值和上一次的观测值共同决定得到的。上一次的估计值可能和观测值不是一个数量级，需要使用观测矩阵H进行变换，将其变到观测的空间中。

协方差更新

$${\mathbf{P}}_{k|k}=\left(\mathbf{I}-K_k{\mathbf{H}}_k\right){\mathbf{P}}_{k|k-1}$$

总结：卡尔曼增益K是做什么用的？

预测和校正中卡尔曼增益K起到关键作用，它决定了预测值和测量值的权重，是测量值的系数。如果测量噪声小，卡尔曼增益会更大，更信任测量值；反之，如果预测更可靠，增益会小，更信任预测值。

测量模型H将状态映射到测量空间，然后计算卡尔曼增益K。卡尔曼增益的公式涉及到预测协方差、测量矩阵H和测量噪声协方差R。通过K，将预测状态和测量值进行加权平均，得到最终的状态估计，并更新协方差矩阵，减少不确定性。

3.2.在线去混响递推滤波器G方法

递推方法不需要直接求y_tilde的自相关矩阵$\mathbf{R}$，而是通过观测值y_tilde和估计值${\mathbf{G}}_{t-1}$来估计当前应该输出的${\mathbf{G}}_t$

卡尔曼滤波的状态估计协方差矩阵$\mathbf{P}$对应$\mathbf{R}$。

卡尔曼滤波的观测值的协方差是$\alpha {\lambda }_{t,f}$。

卡尔曼滤波的观测矩阵$\mathbf{H}$对应y_tilde。

论文中给出下面三个公式

$${\mathbf{K}}_{t,f}=\frac{{\mathbf{R}}_{t-1,f}^{-1}{\tilde{\mathbf{y}}}_{t-\Delta ,f}}{\alpha {\lambda }_{t,f}+{\tilde{\mathbf{y}}}_{t-\Delta ,f}^{\mathrm{H}}{\mathbf{R}}_{t-1,f}^{-1}{\tilde{\mathbf{y}}}_{t-\Delta ,f}}\in {\mathbb{C}}^{DK\times 1}$$

$${\mathbf{R}}_{t,f}^{-1}=\frac{{\mathbf{R}}_{t-1,f}^{-1}-{\mathbf{K}}_{t,f}{\tilde{\mathbf{y}}}_{t-\Delta ,f}^{\mathrm{H}}{\mathbf{R}}_{t-1,f}^{-1}}{\alpha }\in {\mathbb{C}}^{DK\times DK}$$

$${\mathbf{G}}_{t,f}={\mathbf{G}}_{t-1,f}+{\mathbf{K}}_{t,f}{\left({\hat{\mathbf{x}}}_{t,f}^{\text{(early)}}\right)}^{\mathrm{H}}\in {\mathbb{C}}^{DK\times D}.(16)$$

对应公式解释参考下面代码的流程和注释
主调用流程

    channels = 8sampling_rate = 16000delay = 3alpha=0.9999taps = 10#frequency_bins = stft_options['size'] // 2 + 1file_template = 'AMI_WSJ20-Array1-{}_T10c0201.wav'signal_list = [sf.read(str('./third_example/nara_wpe/data/'+file_template.format(d + 1)))[0]for d in range(channels)]y = np.stack(signal_list, axis=0)# 每段的点数nperseg = 512# 重叠点数 每段错开128个点noverlap = 512-128window = np.blackman(nperseg)f, t, Y = stft(y, sampling_rate, window=window, nperseg=nperseg, noverlap=noverlap, axis=-1)# 时间, 频率, mic数    (998, 257, 8)Y = Y.transpose(2,1,0)T, _, _ = Y.shapedef acquire_framebuffer():# buffer大小(13, 257, 8) 作为一个队列 每次返回(14, 257, 8) 14个时间内8个mic的dftbuffer = list(Y[:taps+delay, :, :])for t in range(taps+delay+1, T):# buffer 尾部增加一个元素buffer.append(Y[t, :, :])yield np.array(buffer)# buffer 头部减少一个元素buffer.pop(0)frequency_bins = nperseg//2 + 1Z_list = []# (257, 8*10, 80) 存储当前估计值R^-1 作用同离线方法的y_tilde自相关R矩阵 # 代表了卡尔曼滤波的估计协方差矩阵P 是结果间各变量的关系 # 自相关矩阵和协方差矩阵的关系# 自相关矩阵:自相关主要用于衡量一个信号在不同时刻取值之间的相关性  求法:每列与其他列队对应项相乘后求和 然后平均# 协方差矩阵:协方差矩阵用于衡量多个随机变量之间的协方差关系   求法:每列都减去列平均数 然后其他列队对应项相乘后求和 然后平均# 如果每列均值为0 那么自相关矩阵和协方差矩阵是一样的# 所以这里使用自相关矩阵R 作为卡尔曼滤波器估计的协方差  在线实际求解过程中不需要进行自相关运算Q = np.stack([np.identity(channels * taps) for a in range(frequency_bins)])# (257, 8*10, 8) 最终估计的得到的结果 也就是当前滤波器参数G 初始化为全0矩阵G = np.zeros((frequency_bins, channels * taps, channels))for Y_step in tqdm(acquire_framebuffer()):# Y_step (14, 257, 8) 14个连续时间 先发生的在前Z, Q, G = online_wpe_step(Y_step,get_power_online(Y_step.transpose(1, 2, 0)),Q,G,alpha=alpha,taps=taps,delay=delay)# (257, 8) 同一个时刻8个mic 每个mic在257个频点上的值Z_list.append(Z)Z_stacked = np.stack(Z_list)Z_stacked = np.asarray(Z_stacked).transpose(2, 1, 0)# z = []# for i in range(Z_stacked.shape[0]):#     Zi = Z_stacked[i, :, :]  #     zi = istft(Zi, sampling_rate, window=window, nperseg=nperseg, noverlap=noverlap)#     z.append(zi[1])z = istft(Z_stacked[0], sampling_rate, window=window, nperseg=nperseg, noverlap=noverlap)file_path = './output/wpe_online_output.wav'sf.write(file_path, z[1], sampling_rate)

卡尔曼估计滤波器系数G，同时在计算观测值和估计值的差值时，直接得到了去掉混响的音频数据。
online_wpe_step去混响的具体流程参见下面代码和注释

def online_wpe_step(input_buffer, power_estimate, inv_cov, filter_taps,alpha, taps, delay):"""nara中在线去混响方法这个函数主要是找到卡尔曼方法中各项和当前数据的对应关系 使用卡尔曼方法中 z_k - H * hat{x_{k|k-1} 使用观测到的音频减去估计的反射音 预测当前无混响音频使用卡尔曼方法预测当前滤波器系数G  G被当作卡尔曼中被估计的量系统协方差矩阵inv_cov 初始给定一个随机矩阵 通过卡尔曼方法递推收敛测量协方差矩阵R 当前观测到数据的功率谱密度One step of online dereverberation.Args:input_buffer: Buffer of shape (taps+delay+1, F, D)power_estimate: Estimate for the current PSD 功率谱密度的估计值inv_cov: Current estimate of R^-1  系统协方差矩阵P 直接估计R的逆矩阵而不是使用R 为了避免求逆计算filter_taps: Current estimate of filter taps (F, taps*D, D) 上次得到的滤波器Galpha (float): Smoothing factor 平滑因子taps (int): Number of filter tapsdelay (int): Delay in framesReturns:Dereverberated frame of shape (F, D) 去混响结果 (频点数, mic数)Updated estimate of R^-1    卡尔曼估计协方差矩阵P更新Updated estimate of the filter taps  卡尔曼上次观测结果 也就是滤波器G"""# F频点数 Dmic通道数 (14, 257, 8)F, D = input_buffer.shape[-2:]# input_buffer[:-delay - 1] (10,257,8) 然后颠倒 时间维度上先发生的放在后面# input_buffer是stft的结果 所以window中存储的是复数# window代表了当前已经发生的历史数据 去混响就是预测window对input_buffer[-1]的影响#       window                 input_buffer[-1]     滤波器G与window作用得到估计混响音# [   通道0历史数据   ] - delay - [通道0当前音频]# [   通道1历史数据   ] - delay - [通道1当前音频]# [     ......      ] - delay - [  ......  ]# [   通道n历史数据   ] - delay - [通道n当前音频]# 已知当前音input_buffer[-1]和 估计混响音# 干净音频pred = input_buffer[-1] - 估计混响音window = input_buffer[:-delay - 1][::-1]# (10,257,8) -> (257,8,10) -> (257, 80) 8个mic 10个时间的数据放到一个维度window = window.transpose(1, 2, 0).reshape((F, taps * D))# 通过预测获取当前去混响后的音频pred = (# input_buffer[-1] 是先发生的 形状(257, 8)# filter_taps(257, 80, 8) window(257,80) -> (257,8)  # 每个频点80个值和window的一个维度80个值相乘并求和 最终得到(257,8)的矩阵# 卡尔曼对应: input_buffer[-1] 是当前的观测值z_k 发出的历史音频数据的stft作为观测值#           window 观测矩阵H 本身是历史音频数据的展开矩阵 (257, 80) 也就是(频点数, taps*通道数) #           filter_taps 是上一次的预测值 对应 hat{x}_{k|k-1} # 观测数据与预测的差: 通过观测矩阵window(之前的数据)将 上次预测值filter_taps(滤波器)转换为预测当前带混响的声音#                  结果(257, 8)        input_buffer[-1] -np.einsum('fid,fi->fd', np.conjugate(filter_taps), window))# 卡尔曼增益分子 P_{k|k-1}*H_k^\top  代表了预测的置信度#              inv_cov是当前系统协方差矩阵P (257, 8*10, 80) #              window 观测矩阵H (257, 80)#              得到 (257, 80) nominator = np.einsum('fij,fj->fi', inv_cov, window)# 卡尔曼分母 H_k P_{k|k-1} H_k^\top + R_k  代表了观测不确定性+观测噪声 总体上体现了观测的不确定性# 中的噪声R_k 直接使用当前自相关数据lambda求出denominator = (alpha * power_estimate).astype(window.dtype)# 求出系统协方差inv_cov经过观测矩阵window后的协方差 inv_cov本身是R的逆矩阵 # 这样得到的结果本身就是逆矩阵 直接与nominator相乘得到卡尔曼增益# denominator是一个1维向量 分母本身是协方差矩阵应该是对角的或者为一个数字 这里为什么可以这样计算?denominator += np.einsum('fi,fi->f', np.conjugate(window), nominator)# 求出卡尔曼增益K 此时分母为一组标量 _stable_positive_inverse(denominator)求出这组标量的倒数kalman_gain = nominator * _stable_positive_inverse(denominator)[:, None]# 更新协方差矩阵P_{k|k}inv_cov_k = inv_cov - np.einsum('fj,fjm,fi->fim',np.conjugate(window),inv_cov,kalman_gain,optimize='optimal')# R^-1inv_cov_k /= alpha# 当前预测的G # hat{x}_{k|k} = hat{x}_{k|k-1} + K * (z_k - H * hat{x_{k|k-1}})filter_taps_k = (filter_taps +np.einsum('fi,fm->fim', kalman_gain, np.conjugate(pred)))return pred, inv_cov_k, filter_taps_k 

下面简单画了一下求卡尔曼系数时的爱因斯坦约定的张量变换过程