CN116052702A - 一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法 - Google Patents

一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法 Download PDF

Info

Publication number
CN116052702A
CN116052702A CN202211647281.0A CN202211647281A CN116052702A CN 116052702 A CN116052702 A CN 116052702A CN 202211647281 A CN202211647281 A CN 202211647281A CN 116052702 A CN116052702 A CN 116052702A
Authority
CN
China
Prior art keywords
signal
noise
covariance matrix
calculating
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202211647281.0A
Other languages
English (en)
Inventor
薛建清
张捷阳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fujian Xingwang Wisdom Software Co ltd
Original Assignee
Fujian Xingwang Wisdom Software Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Fujian Xingwang Wisdom Software Co ltd filed Critical Fujian Xingwang Wisdom Software Co ltd
Priority to CN202211647281.0A priority Critical patent/CN116052702A/zh
Publication of CN116052702A publication Critical patent/CN116052702A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0208—Noise filtering
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0208—Noise filtering
    • G10L21/0216—Noise filtering characterised by the method used for estimating noise
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0208—Noise filtering
    • G10L21/0216—Noise filtering characterised by the method used for estimating noise
    • G10L21/0232—Processing in the frequency domain
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0208—Noise filtering
    • G10L2021/02082—Noise filtering the noise being echo, reverberation of the speech
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00—Reducing energy consumption in communication networks
    • Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Computational Linguistics (AREA)
  • Quality & Reliability (AREA)
  • Signal Processing (AREA)
  • Health & Medical Sciences (AREA)
  • Audiology, Speech & Language Pathology (AREA)
  • Human Computer Interaction (AREA)
  • Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Multimedia (AREA)
  • Circuit For Audible Band Transducer (AREA)

Abstract

本发明提供一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法，包括：采集信号，并将采集信号经过预处理得到短时傅立叶域上的信号；计算多通道噪声协方差矩阵；利用上一帧估计的经过延迟的带混响无噪声信号和当前帧的采集信号估计多通道自回归参数，并基于上一帧的声源变化检测结果决定卡尔曼状态噪声的方差值；利用估计的自回归参数、当前帧的采集信号以及估计的多通道噪声协方差矩阵估计无噪声带混响信号；将预估的无噪声带混响信号经过延迟，并结合自回归系数计算得到预估的无噪晚期混响信号；将无噪声带混响信号减去无噪声晚期混响信号得到期望的直达声以及早期混响信号。本发明降低了运算复杂度，实现在嵌入式产品的实时应用。

Description

一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法
技术领域
本发明涉及音频处理技术领域，尤其涉及一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法。
背景技术
在使用传声器在房间内采集到的说话人的语音信号时,会同时采集到房间壁面的反射声,这些反射声被称为混响。当混响时间较长时,混响会影响语音通信中语音的清晰度,也会降低语音识别系统的识别率。
谱减法可以用来实现语音去混响，例如现有文献“Lebart K,Boucher JM,DenbighP N.A New Method Based on Spectral Subtraction for Speech Dereverberation[J].Acta Acustica United with Acustica,2001,87(3):359-366.”用短时傅里叶变换将单通道语音信号变换到时-频域,然后使用谱减法将当前帧的语音信号功率谱与估计的晚期混响功率谱相减,得到去混响信号的功率谱,最后通过短时傅里叶逆变换得到时域的去混响语音信号。然而,这种基于谱减法的去混响方法对语音音质有较大程度的损伤。
卡尔曼滤波是一种自适应滤波方法，将卡尔曼滤波与多通道预测模型相结合,可以用来作为自适应去混响的方法。例如文献“BraunS,HabetsEAP.Online Dereverberationfor Dynamic Scenarios Using a Kalman Filter With an Autoregressive Model[J].IEEE Signal Processing Letters,2016,23(12):1741-1745.”指出卡尔曼滤波具有较好的去混响性能。
上述多通道去混响算法的模型建立中都没有对环境噪声的存在进行假设。在实际应用中，环境噪声的存在影响了多通道去混响算法的性能。在低频段的语音信号通常由于混杂着环境噪声而出现可预测性，导致低频的语音信号受到过抑制。在论文“MasahitoTogami，MULTICHANNEL ONLINE SPEECH DEREVERBERATION UNDERNOISY ENVIRONMENTS”中提出不使用带噪声的混响信号来估计晚期混响信号。而是通过无噪声的麦克风输入信号来优化去混响滤波器，从而在有噪声的环境下来获得良好的去混响滤波器。对于麦克风采集信号，先使用了多通道维纳滤波器对无噪声的输入信号进行估计；其次再更新多通道去混响滤波器来预测晚期混响信号。该算法中的维纳滤波器依赖去混响滤波器的自回归参数。但是混响路径是时变的，上一帧的去混响滤波参数已经不适用于当前帧的环境。因此该算法存在因果错误，导致降噪能力不足。
在现有的论文“Sebastian Braun,Linear prediction based onlinedereverberation and noise reduction using alternating Kalman filters”中考虑到混响路径的时变性，解决了多通道降噪去混响算法的因果问题，提出了一种先对混响自回归参数估计再进行噪声抑制的顺序结构。能够较好地在低信噪比环境下进行去混响降噪。然而，该算法利用了两个交替的卡尔曼滤波，计算量较大难以实现嵌入式设备的实时应用，其也没有给出多通道噪声协方差矩阵的估计方法，且该算法利用了位矩阵作为状态转移矩阵,在声源位置发生突变之后,卡尔曼滤波的性能明显下降。
现有文献“T.Dietzen,S.Doclo,A.Spriet,W.Tirry,M.Moonen,and T.vanWaterschoot,“Low-Complexity Kalman filter for multi-channel linear-prediction-basedblind speechdereverberation,”in 2017IEEE Workshop onApplications of Signal Processing to Audio and Acoustics(WASPAA),2017,pp.284–288.”提到使用了一个小于1的常数乘单位矩阵作为状态转移矩阵,没有出现声源位置改变前后算法性能相差显著的现象,但是算法总体的性能较差。在实际的应用中,说话者在语音间隙的移动或者说话者的改变都会导致声源位置发生突变的场景的产生。因此,使用卡尔曼滤波去混响时,需要解决好声源位置发生突变带来的问题。
总之，现有技术中多通道去混响算法计算复杂度较高，难以在嵌入式产品中实时应用，在低信噪比环境下表现较差，且在混响环境改变时，现有技术跟踪混响的能力较差，需要较长时间才能收敛。
发明内容
本发明要解决的技术问题，在于提供一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法，实现了混响路径的快速收敛，解决了噪声环境下的去混响问题，降低了运算复杂度，可满足在嵌入式产品的实时应用。
本发明要解决的技术问题是这样实现的：一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法，包括如下步骤：
步骤S1、采集信号，并将采集信号经过预处理得到短时傅立叶域上的信号；
步骤S2、利用短时傅立叶域的采集信号，计算多通道噪声协方差矩阵；
步骤S3、利用上一帧估计的经过延迟的带混响无噪声信号以及当前帧的采集信号，基于卡尔曼滤波算法估计时变的多通道自回归参数；所述多通道自回归参数的估计过程中所用到的卡尔曼状态噪声的方差值根据上一帧的声源变化检测结果调整；
步骤S4、利用估计的自回归参数、当前帧的采集信号以及估计的多通道噪声协方差矩阵估计无噪声带混响信号；
步骤S5、将预估的无噪声带混响信号经过延迟，并结合自回归系数计算得到预估的无噪晚期混响信号；从预估的无噪声带混响信号中减去预估的无噪声晚期混响信号得到期望的直达声以及早期混响信号。
进一步的，所述步骤S1具体为：
假设混响环境下有未知数量的声源，且使用固定在任意位置的M个子麦进行采集，则给出采集信号的stft域的表达：
y(k,n)＝[Y1(k,n),…,YM(k,n)]T
其中，Ym(k,n)是第m路信号第k个子频带，第n帧的频域表达；
假设多通道麦克风信号有两部分组成：
y(k,n)＝x(k,n)+v(k,n)
其中，向量x(k,n)以及v(k,n)分别表示阵列上麦克风采集的混响语音信号以及加性噪声；
所述混响语音信号x(k,n)表达式如下：
Figure BDA0004010242620000041
其中，向量s(k,n)＝[S1(k,n),…,SM(k,n)]T表示采集信号中希望获取的直达声和早期混响信号的stft域，Sm(k,n)表示第m个子麦的第n帧第k个子频带频域表达，矩阵Cl(k,n)∈CM×M，表示对于第n-l,l∈[D,D+1,…,L]帧的采集信号stft域x(k,n-l)∈CM,1的滤波参数；D为延迟参数，L表示滤波器长度，r(k,n)为晚期混响信号。
进一步的，所述步骤S2中多通道噪声协方差矩阵的计算过程具体如下：
步骤a1、预设瞬时后验信噪比阈值φ0以及长时后验信噪比阈值
Figure BDA0004010242620000042
步骤a2、初始化采集信号协方差矩阵
Figure BDA0004010242620000043
以及噪声协方差矩阵
Figure BDA0004010242620000044
步骤a3、在算法的前Linit帧内，假设初期采集信号只有噪声信号，所述Linit为音频初期纯噪声的帧数；
迭代计算采集信号协方差矩阵以及噪声协方差矩阵:
Figure BDA0004010242620000045
Figure BDA0004010242620000046
其中，αv为噪声信号的迭代系数；αy为采集信号的迭代系数，H代表矩阵共厄转置操作；
步骤a4、在Linit帧之后，进行如下计算：
步骤a41、迭代估计采集信号协方差矩阵：
Figure BDA0004010242620000047
步骤a42、考虑到语音信号和噪声信号的非相关性，估计语音信号协方差矩阵：
Figure BDA0004010242620000048
步骤a43、计算瞬时后验信噪比：
Figure BDA0004010242620000049
步骤a44、计算长时后验信噪比：
Figure BDA0004010242620000051
其中，tr{.}代表矩阵求迹操作；
步骤a45、计算先验信噪比：
Figure BDA0004010242620000052
其中，M表示通道数即子麦的个数；
Figure BDA0004010242620000053
步骤a46、计算平滑迭代语音存在概率:
计算local尺度的语音不存在概率:
Figure BDA0004010242620000054
计算加窗平滑的后验信噪比并计算平滑后的语音不存在概率：
Figure BDA0004010242620000055
其中wglobal表示汉宁窗函数，窗长定义为2K1+1；
Figure BDA0004010242620000056
计算第n帧各频点的后验信噪比均值，并计算帧尺度的语音不存在概率：
Figure BDA0004010242620000057
Figure BDA0004010242620000058
结合三个尺度计算语音不存在概率
Figure BDA0004010242620000059
Figure BDA00040102426200000510
基于估计的语音不存在概率，计算多通道先验语音存在概率
Figure BDA00040102426200000511
Figure BDA00040102426200000512
步骤a47、计算平滑迭代语音存在概率：
Figure BDA0004010242620000061
其中，αp表示语音存在概率的平滑系数；
步骤a48、基于语音存在概率决定噪声协方差矩阵估计的平滑系数，并更新多通道噪声协方差矩阵：
Figure BDA0004010242620000062
Figure BDA0004010242620000063
其中，
Figure BDA0004010242620000064
为噪声协方差矩阵，
Figure BDA0004010242620000065
为上一帧估计的噪声协方差矩阵；
至此，噪声协方差矩阵估计完成。
进一步的，所述步骤3中的自回归参数的估计具体如下：
步骤31、建立第一个卡尔曼模型：
构造卡尔曼观测矩阵为：
Figure BDA0004010242620000066
其中
Figure BDA0004010242620000067
表示Kronecker乘积，IM表示M维的单位矩阵,上标T表示向量转置的操作，x(n)表示第n帧无噪声混响信号；同时定义自回归参数作为卡尔曼模块的状态向量为：
c(n)＝Vec{[CL(n)…CD(n)]T}；
CL(n)是状态向量中的一部分，下标L表示该部分是对第(n-L)帧的自回归参数，中间省略号表示省略了对应(n-L)到(n-D)帧的自回归参数，Vec{.}是矩阵拉直操作，表示将大括号的矩阵的列按照从左到右的顺序首尾拼接起来，得到一个新的向量c(n)，c(n)的长度Lc＝M×M×(L-D)；X(n)是一个形状为M×Lc的稀疏矩阵；
步骤32、第一个卡尔曼滤波模块计算步骤：
步骤321、计算先验状态误差协方差：
Figure BDA0004010242620000068
其中，φw(n)表示状态噪声协方差；
步骤322、计算状态误差e(n)：
e(n)＝y(n)-X(n-D)c(n-1)；
其中，y(n)表示麦克风采集信号，X(n-D)表示观测矩阵，c(n-1)表示上一帧计算得到的自回归参数；
步骤323、计算卡尔曼增益K(n)：
Figure BDA0004010242620000071
其中，
Figure BDA0004010242620000072
表示观测噪声协方差；
步骤324、计算后验状态误差协方差
Figure BDA0004010242620000073
Figure BDA0004010242620000074
步骤325、计算自回归参数c(n)：
c(n)＝c(n-1)+K(n)e(n)；
步骤326、计算观测噪声
Figure BDA0004010242620000075
Figure BDA0004010242620000076
进一步的，所述状态噪声协方差φw(n)通过如下方式获得：
状态噪声协方差φw(n)的大小根据相邻两帧的自回归参数的变化量来决定，同时增加一个极小的正数来模拟当相邻两帧估计自回归参数没有变化时，其实际的连续变化性；
Figure BDA0004010242620000077
φw为状态噪声方差,Lc为自回归参数c(n)的长度，
Figure BDA0004010242620000078
为Lc阶的单位矩阵；
所述观测噪声协方差φu(n)通过如下方式获得：
计算先验观测噪声协方差矩阵：
Figure BDA0004010242620000079
结合上一帧计算的后验观测噪声协方差矩阵
Figure BDA00040102426200000710
计算当前帧的观测噪声协方差矩阵：
Figure BDA00040102426200000711
更新后验观测噪声协方差矩阵：
Figure BDA00040102426200000712
其中，α表示迭代系数，后验观测噪声协方差矩阵的初始值定义为全为0的矩阵，
Figure BDA0004010242620000081
为观测噪声；
计算观测噪声协方差φu(n)：
Figure BDA0004010242620000082
进一步的，所述步骤S3中“所述多通道自回归参数的估计过程中所用到的卡尔曼状态噪声的方差值根据上一帧的声源变化检测结果调整”具体如下：
计算采集信号能量以及估计的直达声以及早期信号的能量比值，当能量比值发生突变时，判定为当前帧声源发生变化(因为帧与帧之间时间差大概是32ms，如果32ms之前判定到声源位置变化，那么大概率当前帧的声源位置也发生变化)；并当检测到声源发生变化时，短暂提升去混响模块卡尔曼推导过程中的状态噪声方差值为原先的十倍，直至两者比值恢复到阈值之上，从而加强状态变化的追踪能力；
其中，所述采集信号的能量计算公式如下：
Figure BDA0004010242620000083
所述估计的去混响后的信号的能量计算公式如下:
Figure BDA0004010242620000084
αpy,αps分别为采集信号能量以及去混响后的信号能量的平滑系数，Py(n)表示第n帧采集信号的能量，Ps(n)表示第n帧的去混响后信号(也就是直达声和早期混响)的能量，K表示stft域的频点个数，y(k,n)指的是频点为k，第n帧的采集信号，
Figure BDA0004010242620000085
表示计算得到的直达声信号的第k个频点，第n帧信号；当两者比值Py(n)/Ps(n)<threshold(threshold根据需要设定，比如取0.65)过低时，说明混响泄漏，判定为当前帧声源发生变化。
进一步的，所述步骤S4中无噪声带混响信号通过创建第二个卡尔曼模型进行估计，具体如下：
步骤S41、建立第二个卡尔曼模型：
根据无噪声带混响信号x(n)，构造第二个卡尔曼滤波的状态向量x(n)，表征L帧内的采集信号中去除噪声后的信号，是一个一维向量，长度为L×M，其中M指的是通道数；
x(n)＝[xT(n-L+1),…,xT(n)]T；
其中，x(n)中的每一个x(l)表示第l帧的无噪声混响信号，是一个长度为M的向量；
同时，定义第二个卡尔曼模块的观测噪声为s(n)，其构造方式为：
s(n)＝[01×M(L-1) sT(n)]T；
基于上一阶段卡尔曼滤波估计的自回归参数构c(n)造状态转移矩阵：
Figure BDA0004010242620000091
(表示C_L到C_D之间还有符号，分别是C_{L-1},C_{L-1}.....直到C_D)
构造观测方程中的观测矩阵H:
H＝[0M×M(L-1) IM]；
其中，IM表示阶数为M的单位矩阵；0M×M(L-1)表示尺寸为M行、M(L-1)列的全0矩阵；
至此，构造第二个卡尔曼滤波的状态转移方程以及观测方程如下：
x(n)＝F(n)x(n-1)+s(n)；
y(n)＝Hx(n)+v(n)；
其中，v(n)表示麦克风采集的噪声信号；y(n)表示采集信号；
步骤S42、计算第二个卡尔曼滤波：
计算先验状态误差协方差矩阵
Figure BDA0004010242620000092
Figure BDA0004010242620000093
其中，Φ s (n)表示状态噪声的协方差矩阵；
Figure BDA0004010242620000094
表示第二个卡尔曼模块的先验状态误差协方差矩阵，
Figure BDA0004010242620000095
表示第二个卡尔曼模块上一帧的后验状态误差协方差矩阵；
计算先验状态向量x(n|n-1):
x(n|n-1)＝F(n)x(n-1)；
其中，F(n)表示状态转移矩阵，x(n-1)表示上一帧的状态向量估计值，x(n|n-1)表示当前帧的先验状态向量；
计算卡尔曼增益Kx(n)：
Figure BDA0004010242620000101
计算状态误差ex(n)：
ex(n)＝y(n)-Hx(n|n-1)；
其中，Kx(n)为第二个卡尔曼模块的卡尔曼增益，
Figure BDA0004010242620000102
表示噪声协方差矩阵，y(n)表示麦克风采集信号，ex(n)表示第二个卡尔曼模块的状态误差；
计算后验状态误差协方差矩阵
Figure BDA0004010242620000103
Figure BDA0004010242620000104
计算状态向量x(n)：
x(n)＝x(n|n-1)+Kx(n)ex(n)；
从x(n)中获取预估的无噪声混响信号
Figure BDA0004010242620000105
x(n)是一个长度为L×M的一维向量，取该向量的最后M长度，即可得到预估的无噪声混响信号
Figure BDA0004010242620000106
进一步的，所述第二个卡尔曼滤波计算中的状态噪声的协方差矩阵Φ s (n)通过如下方式获得：
对每一帧信号估计先验协方差矩阵；接着和上一帧计算的后验协方差矩阵结合来对当前帧的协方差矩阵进行估计：
Figure BDA0004010242620000107
其中，γ参数是权衡先验和后验配比的权重参数，Φ s (n)表示状态噪声的协方差矩阵，
Figure BDA0004010242620000108
表示上一帧的后验状态噪声协方差矩阵，
Figure BDA0004010242620000109
表示当前帧估计的先验状态噪声协方差矩阵；
所述后验协方差矩阵通过在时间帧之间进行平滑迭代获得：
Figure BDA00040102426200001010
其中，α表示平滑系数，
Figure BDA00040102426200001011
表示估计的早期混响信号和直达声信号；
所述先验协方差矩阵通过多通道维纳滤波算法来获得：
Figure BDA0004010242620000111
其中
Figure BDA0004010242620000112
为[M×M]尺寸的维纳滤波权重矩阵，其计算方式为:
Figure BDA0004010242620000113
其中，
Figure BDA0004010242620000114
为噪声信号的协方差矩阵；
Φr(n)为无噪声的晚期混响信号的协方差矩阵，通过第一个卡尔曼滤波的产物可以迭代计算得到；Φy(n)为麦克风采集信号的协方差矩阵，同样的也是通过时间帧之间的平滑迭代得到。
进一步的，所述步骤5具体为：
将预估的无噪声带混响信号
Figure BDA0004010242620000115
(即无噪声的采集信号)经过D帧延迟得到
Figure BDA0004010242620000116
将延迟后的无噪声带混响信号和混响自回归参数
Figure BDA0004010242620000117
乘积得到预估的晚期混响信号：
Figure BDA0004010242620000118
从预估的无噪声带混响信号
Figure BDA0004010242620000119
中减去预估的晚期混响信号
Figure BDA00040102426200001110
得到估计的早期混响信号和直达声信号
Figure BDA00040102426200001111
Figure BDA00040102426200001112
本发明具有如下优点：
1、本发明提出了简化算法，通过将第一个卡尔曼滤波状态向量误差协方差矩阵以及观测噪声协方差矩阵近似完全对角化，避免了大尺度的矩阵求逆步骤，极大降低了运算复杂度，从而可实现在嵌入式产品中实时进行去混响降噪处理；
2、本发明提出了一种检测混响环境改变的技术，基于去混响模块前后的能量比值实现声源变化检测模块，当位置变化时增大状态转换噪声方差来加快收敛速度，加强算法混响跟踪能力，使得算法在100ms内完成混响路径的收敛；
3、本发明结合了噪声与混响路径的时变性，构造优先进行混响路径的估计再进行降噪处理的级联顺序，解决在低信噪比环境下的去混响问题。
附图说明
下面参照附图结合实施例对本发明作进一步的说明。
图1为本发明一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法的执行流程图。
图2为本发明整体信号流示意图。
图3为本发明双卡尔曼去混响降噪算法原理示意图。
具体实施方式
如图1至图3所示，本发明提供的一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法，包括如下步骤：
步骤S1、采集信号，并将采集信号经过预处理得到短时傅立叶域上的信号；
步骤S2、利用短时傅立叶域的采集信号，计算多通道噪声协方差矩阵；
步骤S3、利用上一帧估计的经过延迟的带混响无噪声信号以及当前帧的采集信号，基于卡尔曼滤波算法估计时变的多通道自回归参数；所述多通道自回归参数的估计过程中所用到的卡尔曼状态噪声的方差值根据上一帧的声源变化检测结果调整；
步骤S4、利用估计的自回归参数、当前帧的采集信号以及估计的多通道噪声协方差矩阵估计无噪声带混响信号；
步骤S5、将预估的无噪声带混响信号经过延迟，并结合自回归系数计算得到预估的无噪晚期混响信号；从预估的无噪声带混响信号中减去预估的无噪声晚期混响信号得到期望的直达声以及早期混响信号。
其中，所述步骤S1具体为：
假设混响环境下有未知数量的声源，且使用固定在任意位置的M个子麦进行采集，则给出采集信号的stft域的表达：
y(k,n)＝[Y1(k,n),…,YM(k,n)]T
其中，Ym(k,n)是第m路信号第k个子频带，第n帧的频域表达；
假设多通道麦克风信号有两部分组成：
y(k,n)＝x(k,n)+v(k,n)
其中，向量x(k,n)以及v(k,n)分别表示阵列上麦克风采集的混响语音信号以及加性噪声；
所述混响语音信号x(k,n)表达式如下：
Figure BDA0004010242620000131
其中，向量s(k,n)＝[S1(k,n),…,SM(k,n)]T表示采集信号中希望获取的直达声和早期混响信号的stft域，Sm(k,n)表示第m个子麦的第n帧第k个子频带频域表达，矩阵Cl(k,n)∈CM×M，表示对于第n-l,l∈[D,D+1,…,L]帧的采集信号stft域x(k,n-l)∈CM,1的滤波参数；D为延迟参数，L表示滤波器长度，r(k,n)为晚期混响信号。
较佳的，所述步骤S2中多通道噪声协方差矩阵的计算过程具体如下：
步骤a1、预设瞬时后验信噪比阈值φ0以及长时后验信噪比阈值
Figure BDA0004010242620000132
步骤a2、初始化采集信号协方差矩阵
Figure BDA0004010242620000133
以及噪声协方差矩阵
Figure BDA0004010242620000134
步骤a3、在算法的前Linit帧内，假设初期采集信号只有噪声信号，所述Linit为音频初期纯噪声的帧数；
迭代计算采集信号协方差矩阵以及噪声协方差矩阵:
Figure BDA0004010242620000135
Figure BDA0004010242620000136
其中，αv为噪声信号的迭代系数；αy为采集信号的迭代系数，H代表矩阵共厄转置操作；
步骤a4、在Linit帧之后，进行如下计算：
步骤a41、迭代估计采集信号协方差矩阵：
Figure BDA0004010242620000137
步骤a42、考虑到语音信号和噪声信号的非相关性，估计语音信号协方差矩阵：
Figure BDA0004010242620000138
步骤a43、计算瞬时后验信噪比：
Figure BDA0004010242620000141
步骤a44、计算长时后验信噪比：
Figure BDA0004010242620000142
其中，tr{.}代表矩阵求迹操作；
步骤a45、计算先验信噪比：
Figure BDA0004010242620000143
其中，M表示通道数即子麦的个数；
Figure BDA0004010242620000144
步骤a46、计算平滑迭代语音存在概率:
计算local尺度的语音不存在概率:
Figure BDA0004010242620000145
计算加窗平滑的后验信噪比并计算平滑后的语音不存在概率：
Figure BDA0004010242620000146
其中wglobal表示汉宁窗函数，窗长定义为2K1+1；
Figure BDA0004010242620000147
计算第n帧各频点的后验信噪比均值，并计算帧尺度的语音不存在概率：
Figure BDA0004010242620000148
Figure BDA0004010242620000149
结合三个尺度计算语音不存在概率
Figure BDA00040102426200001410
Figure BDA00040102426200001411
基于估计的语音不存在概率，计算多通道先验语音存在概率
Figure BDA00040102426200001412
Figure BDA0004010242620000151
步骤a47、计算平滑迭代语音存在概率：
Figure BDA0004010242620000152
其中，αp表示语音存在概率的平滑系数；
步骤a48、基于语音存在概率决定噪声协方差矩阵估计的平滑系数，并更新多通道噪声协方差矩阵：
Figure BDA0004010242620000153
Figure BDA0004010242620000154
其中，
Figure BDA0004010242620000155
为噪声协方差矩阵，
Figure BDA0004010242620000156
为上一帧估计的噪声协方差矩阵；
至此，噪声协方差矩阵(即是第二个卡尔曼模型中观测噪声协方差矩阵)的估计完成。
较佳的，所述步骤3中的自回归参数的估计具体如下：
步骤31、建立第一个卡尔曼模型：
构造卡尔曼观测矩阵为：
Figure BDA0004010242620000157
其中
Figure BDA0004010242620000158
表示Kronecker乘积，IM表示M维的单位矩阵,上标T表示向量转置的操作，x(n)表示第n帧无噪声混响信号；同时定义自回归参数作为卡尔曼模块的状态向量为：
c(n)＝Vec{[CL(n)…CD(n)]T}；
CL(n)是状态向量中的一部分，下标L表示该部分是对第(n-L)帧的自回归参数，中间省略号表示省略了对应(n-L)到(n-D)帧的自回归参数，Vec{.}是矩阵拉直操作，表示将大括号的矩阵的列按照从左到右的顺序首尾拼接起来，得到一个新的向量c(n)，c(n)的长度Lc＝M×M×(L-D)；X(n)是一个形状为M×Lc的稀疏矩阵；
步骤32、第一个卡尔曼滤波模块计算步骤：
步骤321、计算先验状态误差协方差：
Figure BDA0004010242620000161
其中，φw(n)表示状态噪声协方差；
步骤322、计算状态误差e(n)：
e(n)＝y(n)-X(n-D)c(n-1)；
其中，y(n)表示麦克风采集信号，X(n-D)表示观测矩阵，对于卡尔曼模块是观测矩阵，对于实际意义来说是无噪混响，c(n-1)表示上一帧计算得到的自回归参数；
步骤323、计算卡尔曼增益K(n)：
Figure BDA0004010242620000162
其中，
Figure BDA0004010242620000163
表示观测噪声协方差；
步骤324、计算后验状态误差协方差
Figure BDA0004010242620000164
Figure BDA0004010242620000165
步骤325、计算自回归参数c(n)：
c(n)＝c(n-1)+K(n)e(n)；
步骤326、计算观测噪声
Figure BDA0004010242620000166
Figure BDA0004010242620000167
较佳的，所述状态噪声协方差φw(n)通过如下方式获得：
状态噪声协方差φw(n)的大小根据相邻两帧的自回归参数的变化量来决定，同时增加一个极小的正数来模拟当相邻两帧估计自回归参数没有变化时，其实际的连续变化性；
Figure BDA0004010242620000168
φw为状态噪声方差,Lc为自回归参数c(n)的长度，
Figure BDA0004010242620000169
为Lc阶的单位矩阵；
所述观测噪声协方差φu(n)通过如下方式获得：
计算先验观测噪声协方差矩阵：
Figure BDA00040102426200001610
结合上一帧计算的后验观测噪声协方差矩阵
Figure BDA0004010242620000171
计算当前帧的观测噪声协方差矩阵：
Figure BDA0004010242620000172
更新后验观测噪声协方差矩阵：
Figure BDA0004010242620000173
其中，α表示迭代系数，后验观测噪声协方差矩阵的初始值定义为全为0的矩阵，
Figure BDA0004010242620000174
为观测噪声；
计算观测噪声协方差φu(n)：
Figure BDA0004010242620000175
由于本申请中对上述第一个卡尔曼模型的计算上进行了简化，从而有效提高计算速率，并保证了最终计算准确度。本申请上述第一个卡尔曼模型简化计算包括将后验状态误差协方差矩阵、先验状态误差协方差矩阵对角化以及语音信号协方差矩阵对角化来进行运算量简化操作，其简化计算的中间推导过程具体如下：
首先，将先验状态误差协方差矩阵近似对角化；将先验状态误差协方差矩阵近似为单位阵乘上一个系数
Figure BDA0004010242620000176
定义状态预测协方差矩阵为:
Figure BDA0004010242620000177
其中，
Figure BDA0004010242620000178
表示矩阵阶数为Lc的单位矩阵，
Figure BDA0004010242620000179
表示先验状态误差协方差；
同样的，将后验状态误差协方差矩阵近似为：
Figure BDA00040102426200001710
其中，
Figure BDA00040102426200001711
表示矩阵阶数为Lc的单位矩阵，
Figure BDA00040102426200001712
表示后验状态协方差；
将现有的卡尔曼滤波黄金五条中的更新矩阵
Figure BDA00040102426200001713
进行近似：
Figure BDA0004010242620000181
其中tr{}为求迹符号；
然后，将观测噪声协方差矩阵近似对角化：即第一个卡尔曼模型中观测噪声协方差矩阵进行近似得到：
Figure BDA0004010242620000182
其中，
Figure BDA0004010242620000183
表示观测噪声协方差，IM表示阶数为M的单位矩阵，e(n)表示状态误差，其中
Figure BDA0004010242620000184
表示取二范数操作后再平方；
其中，后验观测噪声协方差矩阵的更新也进行了近似处理，近似为：
Figure BDA0004010242620000185
其中，
Figure BDA0004010242620000186
表示后验观测噪声协方差，IM为阶数为M的单位矩阵；
进而，得到本发明的上述简化后的卡尔曼滤波计算流程为：
计算先验状态误差协方差
Figure BDA0004010242620000187
Figure BDA0004010242620000188
计算状态误差e(n)：
e(n)＝y(n)-X(n-D)c(n-1)；
计算卡尔曼增益K(n)：
Figure BDA0004010242620000189
计算后验状态误差协方差
Figure BDA00040102426200001810
Figure BDA0004010242620000191
计算自回归参数c(n)：c(n)＝c(n|n-1)+K(n)e(n)；
并计算观测噪声
Figure BDA0004010242620000192
Figure BDA0004010242620000193
上述即为步骤S3中第一个卡尔曼滤波的简化流程的详细推导过程。
较佳的，所述步骤S3中“所述多通道自回归参数的估计过程中所用到的卡尔曼状态噪声的方差值根据上一帧的声源变化检测结果调整”具体如下：
计算采集信号能量以及估计的直达声以及早期信号的能量比值，当能量比值发生突变时，判定为当前帧声源发生变化(因为帧与帧之间时间差大概是32ms，如果32ms之前判定到声源位置变化，那么大概率当前帧的声源位置也发生变化)；并当检测到声源发生变化时，短暂提升去混响模块卡尔曼推导过程中的状态噪声方差值为原先的十倍，直至两者比值恢复到阈值之上，从而加强状态变化的追踪能力；
其中，所述采集信号的能量计算公式如下：
Figure BDA0004010242620000194
所述估计的去混响后的信号的能量计算公式如下:
Figure BDA0004010242620000195
αpy,αps分别为采集信号能量以及去混响后的信号能量的平滑系数，Py(n)表示第n帧采集信号的能量，Ps(n)表示第n帧的去混响后信号(也就是直达声和早期混响)的能量，K表示stft域的频点个数，Py(n-1)指的是n-1帧采集信号的能量，y(k,n)指的是频点为k，第n帧的采集信号，
Figure BDA0004010242620000196
表示计算得到的直达声信号的第k个频点，第n帧信号，这里由于需要将所有的频点k进行累加，因此，引入频点k进行计算，若不引入频点k，则表示所有的频点k都按照同一方法计算；当两者比值Py(n)/Ps(n)<threshold(threshold根据需要设定，比如取0.65)过低时，说明混响泄漏，判定为当前帧声源发生变化。
本发明的第一个卡尔曼模块参数说明表如下表1所示：
表1
Figure BDA0004010242620000201
较佳的，所述步骤S4中无噪声带混响信号通过创建第二个卡尔曼模型进行估计，具体如下：
步骤S41、建立第二个卡尔曼模型：
根据无噪声带混响信号x(n)，构造第二个卡尔曼滤波的状态向量x(n)，表征L帧内的采集信号中去除噪声后的信号，是一个一维向量，长度为L×M，其中M指的是通道数；
x(n)＝[xT(n-L+1),…,xT(n)]T；
其中，x(n)中的每一个x(l)表示第l帧的无噪声混响信号，是一个长度为M的向量；
同时，定义第二个卡尔曼模块的观测噪声为s(n)，其构造方式为：
s(n)＝[01×M(L-1)sT(n)]T；
基于上一阶段卡尔曼滤波估计的自回归参数构c(n)造状态转移矩阵：
Figure BDA0004010242620000211
(表示C_L到C_D之间还有符号，分别是C_{L-1},C_{L-1}.....直到C_D)
构造观测方程中的观测矩阵H:
H＝[0M×M(L-1)IM]；
其中，IM表示阶数为M的单位矩阵；0M×M(L-1)表示尺寸为M行、M(L-1)列的全0矩阵；
至此，构造第二个卡尔曼滤波的状态转移方程以及观测方程如下：
x(n)＝F(n)x(n-1)+s(n)；
y(n)＝Hx(n)+v(n)；
其中，v(n)表示麦克风采集的噪声信号；y(n)表示采集信号；
步骤S42、计算第二个卡尔曼滤波：
计算先验状态误差协方差矩阵
Figure BDA0004010242620000212
Figure BDA0004010242620000213
其中，Φ s (n)表示状态噪声的协方差矩阵；
Figure BDA0004010242620000214
表示第二个卡尔曼模块的先验状态误差协方差矩阵，
Figure BDA0004010242620000215
表示第二个卡尔曼模块上一帧的后验状态误差协方差矩阵；
计算先验状态向量x(n|n-1):
x(n|n-1)＝F(n)x(n-1)；
其中，F(n)表示状态转移矩阵，x(n-1)表示上一帧的状态向量估计值，x(n|n-1)表示当前帧的先验状态向量；
计算卡尔曼增益Kx(n)：
Figure BDA0004010242620000216
计算状态误差ex(n)：
ex(n)＝y(n)-Hx(n|n-1)；
其中，Kx(n)为第二个卡尔曼模块的卡尔曼增益，
Figure BDA0004010242620000217
表示噪声协方差矩阵，y(n)表示麦克风采集信号，ex(n)表示第二个卡尔曼模块的状态误差；
计算后验状态误差协方差矩阵
Figure BDA0004010242620000221
Figure BDA0004010242620000222
计算状态向量x(n)：
x(n)＝x(n|n-1)+Kx(n)ex(n)；
从x(n)中获取预估的无噪声混响信号
Figure BDA0004010242620000223
x(n)是一个长度为L×M的一维向量，取该向量的最后M长度，即可得到预估的无噪声混响信号
Figure BDA0004010242620000224
较佳的，所述第二个卡尔曼滤波计算中的状态噪声的协方差矩阵Φ s (n)通过如下方式获得：
对每一帧信号估计先验协方差矩阵；接着和上一帧计算的后验协方差矩阵结合来对当前帧的协方差矩阵进行估计：
Figure BDA0004010242620000225
其中，γ参数是权衡先验和后验配比的权重参数，Φ s (n)表示状态噪声的协方差矩阵，
Figure BDA0004010242620000226
表示上一帧的后验状态噪声协方差矩阵，
Figure BDA0004010242620000227
表示当前帧估计的先验状态噪声协方差矩阵；
所述后验协方差矩阵通过在时间帧之间进行平滑迭代获得：
Figure BDA0004010242620000228
其中，α表示平滑系数，
Figure BDA0004010242620000229
表示估计的早期混响信号和直达声信号；
所述先验协方差矩阵通过多通道维纳滤波算法来获得：
Figure BDA00040102426200002210
其中
Figure BDA00040102426200002211
为[M×M]尺寸的维纳滤波权重矩阵，其计算方式为:
Figure BDA00040102426200002212
其中，
Figure BDA00040102426200002213
为噪声信号的协方差矩阵；
Φr(n)为无噪声的晚期混响信号的协方差矩阵，通过第一个卡尔曼滤波的产物可以迭代计算得到；Φy(n)为麦克风采集信号的协方差矩阵，同样的也是通过时间帧之间的平滑迭代得到。
较佳的，所述步骤5具体为：
将预估的无噪声带混响信号
Figure BDA0004010242620000231
(即无噪声的采集信号)经过D帧延迟得到
Figure BDA0004010242620000232
将延迟后的无噪声带混响信号和混响自回归参数
Figure BDA0004010242620000233
乘积得到预估的晚期混响信号：
Figure BDA0004010242620000234
从预估的无噪声带混响信号
Figure BDA0004010242620000235
中减去预估的晚期混响信号
Figure BDA0004010242620000236
得到估计的早期混响信号和直达声信号
Figure BDA0004010242620000237
Figure BDA0004010242620000238
本发明上述实施例提供的技术方案至少具有如下优点：
简化算法，通过将第一个卡尔曼滤波状态向量误差协方差矩阵以及观测噪声协方差矩阵近似完全对角化，避免了大尺度的矩阵求逆步骤，极大降低了运算复杂度，从而可实现在嵌入式产品中实时进行去混响降噪处理；提出了一种检测混响环境改变的技术，基于去混响模块前后的能量比值实现声源变化检测模块，当位置变化时增大状态转换噪声方差来加快收敛速度，加强算法混响跟踪能力，使得算法在100ms内完成混响路径的收敛；结合了噪声与混响路径的时变性，构造优先进行混响路径的估计再进行降噪处理的级联顺序，解决在低信噪比环境下的去混响问题。
虽然以上描述了本发明的具体实施方式，但是熟悉本技术领域的技术人员应当理解，我们所描述的具体的实施例只是说明性的，而不是用于对本发明的范围的限定，熟悉本领域的技术人员在依照本发明的精神所作的等效的修饰以及变化，都应当涵盖在本发明的权利要求所保护的范围内。

Claims (9)

1.一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法，其特征在于：包括如下步骤：
步骤S1、采集信号，并将采集信号经过预处理得到短时傅立叶域上的信号；
步骤S2、利用短时傅立叶域的采集信号，计算多通道噪声协方差矩阵；
步骤S3、利用上一帧估计的经过延迟的带混响无噪声信号以及当前帧的采集信号，基于卡尔曼滤波算法估计时变的多通道自回归参数；所述多通道自回归参数的估计过程中所用到的卡尔曼状态噪声的方差值根据上一帧的声源变化检测结果调整；
步骤S4、利用估计的自回归参数、当前帧的采集信号以及估计的多通道噪声协方差矩阵估计无噪声带混响信号；
步骤S5、将预估的无噪声带混响信号经过延迟，并结合自回归系数计算得到预估的无噪晚期混响信号；从预估的无噪声带混响信号中减去预估的无噪声晚期混响信号得到期望的直达声以及早期混响信号。
2.根据权利要求1所述的方法，其特征在于：所述步骤S1具体为：
假设混响环境下有未知数量的声源，且使用固定在任意位置的M个子麦进行采集，则给出采集信号的stft域的表达：
y(k,n)＝[Y1(k,n),…,YM(k,n)]T
其中，Ym(k,n)是第m路信号第k个子频带，第n帧的频域表达；
假设多通道麦克风信号有两部分组成：
y(k,n)＝x(k,n)+v(k,n)
其中，向量x(k,n)以及v(k,n)分别表示阵列上麦克风采集的混响语音信号以及加性噪声；
所述混响语音信号x(k,n)表达式如下：
Figure FDA0004010242610000011
其中，向量s(k,n)＝[S1(k,n),…,SM(k,n)]T表示采集信号中希望获取的直达声和早期混响信号的stft域，Sm(k,n)表示第m个子麦的第n帧第k个子频带频域表达，矩阵Cl(k,n)∈CM×M，表示对于第n-l,l∈[D,D+1,…,L]帧的采集信号stft域x(k,n-l)∈CM,1的滤波参数；D为延迟参数，L表示滤波器长度，r(k,n)为晚期混响信号。
3.根据权利要求1所述的方法，其特征在于：所述步骤S2中多通道噪声协方差矩阵的计算过程具体如下：
步骤a1、预设瞬时后验信噪比阈值φ0以及长时后验信噪比阈值
Figure FDA0004010242610000021
步骤a2、初始化采集信号协方差矩阵
Figure FDA0004010242610000022
以及噪声协方差矩阵
Figure FDA0004010242610000023
步骤a3、在算法的前Linit帧内，假设初期采集信号只有噪声信号，所述Linit为音频初期纯噪声的帧数；
迭代计算采集信号协方差矩阵以及噪声协方差矩阵:
Figure FDA0004010242610000024
Figure FDA0004010242610000025
其中，αv为噪声信号的迭代系数；αy为采集信号的迭代系数，H代表矩阵共厄转置操作；
步骤a4、在Linit帧之后，进行如下计算：
步骤a41、迭代估计采集信号协方差矩阵：
Figure FDA0004010242610000026
步骤a42、考虑到语音信号和噪声信号的非相关性，估计语音信号协方差矩阵：
Figure FDA0004010242610000027
步骤a43、计算瞬时后验信噪比：
Figure FDA0004010242610000028
步骤a44、计算长时后验信噪比：
Figure FDA0004010242610000029
其中，tr{.}代表矩阵求迹操作；
步骤a45、计算先验信噪比：
Figure FDA0004010242610000031
其中，M表示通道数即子麦的个数；
Figure FDA0004010242610000032
步骤a46、计算平滑迭代语音存在概率:
计算local尺度的语音不存在概率:
Figure FDA0004010242610000033
计算加窗平滑的后验信噪比并计算平滑后的语音不存在概率：
Figure FDA0004010242610000034
其中wglobal表示汉宁窗函数，窗长定义为2K1+1；
Figure FDA0004010242610000035
计算第n帧各频点的后验信噪比均值，并计算帧尺度的语音不存在概率：
Figure FDA0004010242610000036
Figure FDA0004010242610000037
结合三个尺度计算语音不存在概率
Figure FDA0004010242610000038
Figure FDA0004010242610000039
基于估计的语音不存在概率，计算多通道先验语音存在概率
Figure FDA00040102426100000310
Figure FDA00040102426100000311
步骤a47、计算平滑迭代语音存在概率：
Figure FDA00040102426100000312
其中，αp表示语音存在概率的平滑系数；
步骤a48、基于语音存在概率决定噪声协方差矩阵估计的平滑系数，并更新多通道噪声协方差矩阵：
Figure FDA0004010242610000041
Figure FDA0004010242610000042
其中，
Figure FDA0004010242610000043
为噪声协方差矩阵，
Figure FDA0004010242610000044
为上一帧估计的噪声协方差矩阵；
至此，噪声协方差矩阵估计完成。
4.根据权利要求1所述的方法，其特征在于：所述步骤3中的自回归参数的估计具体如下：
步骤31、建立第一个卡尔曼模型：
构造卡尔曼观测矩阵为：
Figure FDA0004010242610000045
其中
Figure FDA0004010242610000046
表示Kronecker乘积，IM表示M维的单位矩阵,上标T表示向量转置的操作，x(n)表示第n帧无噪声混响信号；同时定义自回归参数作为卡尔曼模块的状态向量为：
c(n)＝Vec{[CL(n)…CD(n)]T}；
CL(n)是状态向量中的一部分，下标L表示该部分是对第(n-L)帧的自回归参数，中间省略号表示省略了对应(n-L)到(n-D)帧的自回归参数，Vec{.}是矩阵拉直操作，表示将大括号的矩阵的列按照从左到右的顺序首尾拼接起来，得到一个新的向量c(n)，c(n)的长度Lc＝M×M×(L-D)；X(n)是一个形状为M×Lc的稀疏矩阵；
步骤32、第一个卡尔曼滤波模块计算步骤：
步骤321、计算先验状态误差协方差：
Figure FDA0004010242610000047
其中，φw(n)表示状态噪声协方差；
步骤322、计算状态误差e(n)：
e(n)＝y(n)-X(n-D)c(n-1)；
其中，y(n)表示麦克风采集信号，X(n-D)表示观测矩阵，c(n-1)表示上一帧计算得到的自回归参数；
步骤323、计算卡尔曼增益K(n)：
Figure FDA0004010242610000051
其中，
Figure FDA0004010242610000052
表示观测噪声协方差；
步骤324、计算后验状态误差协方差
Figure FDA0004010242610000053
Figure FDA0004010242610000054
步骤325、计算自回归参数c(n)：
c(n)＝c(n-1)+K(n)e(n)；
步骤326、计算观测噪声
Figure FDA0004010242610000055
Figure FDA0004010242610000056
5.根据权利要求4所述的方法，其特征在于：所述状态噪声协方差φw(n)通过如下方式获得：
状态噪声协方差φw(n)的大小根据相邻两帧的自回归参数的变化量来决定，同时增加一个极小的正数来模拟当相邻两帧估计自回归参数没有变化时，其实际的连续变化性；
Figure FDA0004010242610000057
φw为状态噪声方差,Lc为自回归参数c(n)的长度，
Figure FDA0004010242610000058
为Lc阶的单位矩阵；
所述观测噪声协方差φu(n)通过如下方式获得：
计算先验观测噪声协方差矩阵：
Figure FDA0004010242610000059
结合上一帧计算的后验观测噪声协方差矩阵
Figure FDA00040102426100000510
计算当前帧的观测噪声协方差矩阵：
Figure FDA00040102426100000511
更新后验观测噪声协方差矩阵：
Figure FDA00040102426100000512
其中，α表示迭代系数，后验观测噪声协方差矩阵的初始值定义为全为0的矩阵，
Figure FDA0004010242610000061
为观测噪声；
计算观测噪声协方差φu(n)：
Figure FDA0004010242610000062
6.根据权利要求1所述的方法，其特征在于：所述步骤S3中“所述多通道自回归参数的估计过程中所用到的卡尔曼状态噪声的方差值根据上一帧的声源变化检测结果调整”具体如下：
计算采集信号能量以及估计的直达声以及早期信号的能量比值，当能量比值发生突变时，判定为当前帧声源发生变化；并当检测到声源发生变化时，短暂提升去混响模块卡尔曼推导过程中的状态噪声方差值为原先的十倍，直至两者比值恢复到阈值之上，从而加强状态变化的追踪能力；
其中，所述采集信号的能量计算公式如下：
Figure FDA0004010242610000063
所述估计的去混响后的信号的能量计算公式如下:
Figure FDA0004010242610000064
αpy,αps分别为采集信号能量以及去混响后的信号能量的平滑系数，Py(n)表示第n帧采集信号的能量，Ps(n)表示第n帧的去混响后信号的能量，K表示stft域的频点个数，y(k,n)指的是频点为k，第n帧的采集信号，
Figure FDA0004010242610000065
表示计算得到的直达声信号的第k个频点，第n帧信号；当两者比值Py(n)/Ps(n)<threshold过低时，说明混响泄漏，判定为当前帧声源发生变化。
7.根据权利要求1所述的方法，其特征在于：所述步骤S4中无噪声带混响信号通过创建第二个卡尔曼模型进行估计，具体如下：
步骤S41、建立第二个卡尔曼模型：
根据无噪声带混响信号x(n)，构造第二个卡尔曼滤波的状态向量x(n)，表征L帧内的采集信号中去除噪声后的信号，是一个一维向量，长度为L×M，其中M指的是通道数；
x(n)＝[xT(n-L+1),…,xT(n)]T；
其中，x(n)中的每一个x(l)表示第l帧的无噪声混响信号，是一个长度为M的向量；
同时，定义第二个卡尔曼模块的观测噪声为s(n)，其构造方式为：
s(n)＝[01×M(L-1)sT(n)]T；
基于上一阶段卡尔曼滤波估计的自回归参数构c(n)造状态转移矩阵：
Figure FDA0004010242610000071
构造观测方程中的观测矩阵H:
H＝[0M×M(L-1)IM]；
其中，IM表示阶数为M的单位矩阵；0M×M(L-1)表示尺寸为M行、M(L-1)列的全0矩阵；
至此，构造第二个卡尔曼滤波的状态转移方程以及观测方程如下：
x(n)＝F(n)x(n-1)+s(n)；
y(n)＝Hx(n)+v(n)；
其中，v(n)表示麦克风采集的噪声信号；y(n)表示采集信号；
步骤S42、计算第二个卡尔曼滤波：
计算先验状态误差协方差矩阵
Figure FDA0004010242610000072
Figure FDA0004010242610000073
其中，Φ s (n)表示状态噪声的协方差矩阵；
Figure FDA0004010242610000074
表示第二个卡尔曼模块的先验状态误差协方差矩阵，
Figure FDA0004010242610000075
表示第二个卡尔曼模块上一帧的后验状态误差协方差矩阵；
计算先验状态向量x(n|n-1):
x(n|n-1)＝F(n)x(n-1)；
其中，F(n)表示状态转移矩阵，x(n-1)表示上一帧的状态向量估计值，x(n|n-1)表示当前帧的先验状态向量；
计算卡尔曼增益Kx(n)：
Figure FDA0004010242610000081
计算状态误差ex(n)：
ex(n)＝y(n)-Hx(n|n-1)；
其中，Kx(n)为第二个卡尔曼模块的卡尔曼增益，
Figure FDA0004010242610000082
表示噪声协方差矩阵，y(n)表示麦克风采集信号，ex(n)表示第二个卡尔曼模块的状态误差；
计算后验状态误差协方差矩阵
Figure FDA0004010242610000083
Figure FDA0004010242610000084
计算状态向量x(n)：
x(n)＝x(n|n-1)+Kx(n)ex(n)；
从x(n)中获取预估的无噪声混响信号
Figure FDA0004010242610000085
x(n)是一个长度为L×M的一维向量，取该向量的最后M长度，即可得到预估的无噪声混响信号
Figure FDA0004010242610000086
8.根据权利要求7所述的方法，其特征在于：所述第二个卡尔曼滤波计算中的状态噪声的协方差矩阵
Figure FDA00040102426100000812
通过如下方式获得：
对每一帧信号估计先验协方差矩阵；接着和上一帧计算的后验协方差矩阵结合来对当前帧的协方差矩阵进行估计：
Figure FDA0004010242610000087
其中，γ参数是权衡先验和后验配比的权重参数，Φ s (n)表示状态噪声的协方差矩阵，
Figure FDA0004010242610000088
表示上一帧的后验状态噪声协方差矩阵，
Figure FDA0004010242610000089
表示当前帧估计的先验状态噪声协方差矩阵；
所述后验协方差矩阵通过在时间帧之间进行平滑迭代获得：
Figure FDA00040102426100000810
其中，α表示平滑系数，
Figure FDA00040102426100000811
表示估计的早期混响信号和直达声信号；
所述先验协方差矩阵通过多通道维纳滤波算法来获得：
Figure FDA0004010242610000091
其中
Figure FDA0004010242610000092
为[M×M]尺寸的维纳滤波权重矩阵，其计算方式为:
Figure FDA0004010242610000093
其中，
Figure FDA0004010242610000094
为噪声信号的协方差矩阵；
Φr(n)为无噪声的晚期混响信号的协方差矩阵，通过第一个卡尔曼滤波的产物可以迭代计算得到；Φy(n)为麦克风采集信号的协方差矩阵，同样的也是通过时间帧之间的平滑迭代得到。
9.根据权利要求1所述的方法，其特征在于：所述步骤5具体为：
将预估的无噪声带混响信号
Figure FDA0004010242610000095
经过D帧延迟得到
Figure FDA0004010242610000096
将延迟后的无噪声带混响信号和混响自回归参数
Figure FDA0004010242610000097
乘积得到预估的晚期混响信号：
Figure FDA0004010242610000098
从预估的无噪声带混响信号
Figure FDA0004010242610000099
中减去预估的晚期混响信号
Figure FDA00040102426100000910
得到估计的早期混响信号和直达声信号
Figure FDA00040102426100000911
Figure FDA00040102426100000912
CN202211647281.0A 2022-12-21 2022-12-21 一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法 Pending CN116052702A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211647281.0A CN116052702A (zh) 2022-12-21 2022-12-21 一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211647281.0A CN116052702A (zh) 2022-12-21 2022-12-21 一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法

Publications (1)

Publication Number Publication Date
CN116052702A true CN116052702A (zh) 2023-05-02

Family

ID=86119200

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211647281.0A Pending CN116052702A (zh) 2022-12-21 2022-12-21 一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法

Country Status (1)

Country Link
CN (1) CN116052702A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117275528A (zh) * 2023-11-17 2023-12-22 浙江华创视讯科技有限公司 语音存在概率的估计方法及装置

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117275528A (zh) * 2023-11-17 2023-12-22 浙江华创视讯科技有限公司 语音存在概率的估计方法及装置
CN117275528B (zh) * 2023-11-17 2024-03-01 浙江华创视讯科技有限公司 语音存在概率的估计方法及装置

Similar Documents

Publication Publication Date Title
CN108172231B (zh) 一种基于卡尔曼滤波的去混响方法及系统
Schwartz et al. Online speech dereverberation using Kalman filter and EM algorithm
Yoshioka et al. Generalization of multi-channel linear prediction methods for blind MIMO impulse response shortening
US8849657B2 (en) Apparatus and method for isolating multi-channel sound source
CN111418012B (zh) 用于处理音频信号的方法和音频处理设备
Mohammadiha et al. Single channel speech enhancement using Bayesian NMF with recursive temporal updates of prior distributions
EP3685378B1 (en) Signal processor and method for providing a processed audio signal reducing noise and reverberation
US11483651B2 (en) Processing audio signals
JP6225245B2 (ja) 信号処理装置、方法及びプログラム
Cord-Landwehr et al. Monaural source separation: From anechoic to reverberant environments
CN110111802B (zh) 基于卡尔曼滤波的自适应去混响方法
JP6748304B2 (ja) ニューラルネットワークを用いた信号処理装置、ニューラルネットワークを用いた信号処理方法及び信号処理プログラム
CN115424627A (zh) 基于卷积循环网络和wpe算法的语音增强混合处理方法
KR20220022286A (ko) 잔향 제거 오토 인코더를 이용한 잔향 환경 임베딩 추출 방법 및 장치
CN116052702A (zh) 一种基于卡尔曼滤波的低复杂度多通道去混响降噪方法
Nesta et al. Robust Automatic Speech Recognition through On-line Semi Blind Signal Extraction
Schwartz et al. Maximum likelihood estimation of the late reverberant power spectral density in noisy environments
Schwartz et al. Multi-microphone speech dereverberation using expectation-maximization and kalman smoothing
Yu et al. Multi-channel $ l_ {1} $ regularized convex speech enhancement model and fast computation by the split bregman method
Yoshioka et al. Dereverberation by using time-variant nature of speech production system
Parchami et al. Speech reverberation suppression for time-varying environments using weighted prediction error method with time-varying autoregressive model
CN113160842B (zh) 一种基于mclp的语音去混响方法及系统
Dionelis On single-channel speech enhancement and on non-linear modulation-domain Kalman filtering
Gao et al. A Physical Model-Based Self-Supervised Learning Method for Signal Enhancement Under Reverberant Environment
Ephraim et al. A brief survey of speech enhancement 1

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination