频谱校正及数据密度聚类法语音欠定盲识别方法和装置

摘要

本发明属于数字信号处理技术领域，为提高算法精度；显著提高算法的效率；使算法对噪声的鲁棒性更强。本发明采用的技术方案是，频谱校正及数据密度聚类法语音欠定盲识别方法，包括如下步骤：Step1：做加汉宁窗L点50％重叠的短时傅立叶变换STFT，得到观测频谱X

说明书

技术领域

本发明属于数字信号处理技术领域。具体涉及对于语音信号在观测数目小于源数的欠定情形下混合矩阵的盲估计问题。

背景技术

盲源分离(BlindSourceSeparation，简称BSS)[1]是在没有信道先验知识的前提下，仅依靠接收到的观测信号估计信道参数，进而恢复源信号的过程。其应用涉及语音信号处理[2]、图像处理[3]、机械故障诊断[4]、通信信道估计[5]等等领域，是信号处理领域的热点问题之一。

按照源信号数目及观测信号的数目的相对关系，盲源分离可分为超定情形(源信号数目N小于观测信号数M)、正定情形(N＝M)及欠定情形(N>M)。显然，欠定情形是最难处理的情形。

盲源分离(BSS)问题通常可分为两个阶段[6]：1)盲识别(Blindidentification,BI)阶段；2)盲恢复(Blindrecovery,BR)阶段。盲识别阶段是指从观测信号中估计信道参数的过程；盲恢复阶段是指利用盲估计的结果恢复源信号的过程。显然，盲识别阶段性能的好坏直接影响后续阶段性能。因而，高效、精确的盲识别算法是盲源分离问题的基础。

盲识别方法主要分为2类：基于统计量的盲识别方法及基于稀疏表示的盲识别方法[3]。

基于统计量的盲识别方法假定源信号相互独立。该类方法始于独立分量分析(IndependentComponentAnalysis,ICA)[6-8]的提出。然而ICA方法仅适用于超定情形。为了解决欠定情形盲估计问题，De.Lathauwer等人提出了基于四阶累积量的盲识别(Fourth-OrderCumulant-BasedBlindIdentification,FOOBI)方法[9]，引入四阶张量分析及矩阵瞬时对角化解决欠定盲识别问题。

基于稀疏表示的盲识别方法假定源信号本身呈现稀疏分布，或者可以找到一个变换域(如傅立叶变换(FastFourierTransform,FFT),短时傅立叶变换(ShortTimeFourierTransform，STFT),魏格纳分布(Wigner-VileDistribution,WVD)等)，使得信号在该域上呈现稀疏分布。文献[10]分析了信号稀疏程度与可分离性之间的关系。Bofill及Zibulevsky提出了基于快速傅立叶变换的势函数法[11]，利用势函数作为目标函数，通过势函数的极值数目及位置可估计源的数目及混合矩阵。文献[12]提出了LOST算法(LineOrientationSeparationTechnique,LOST)，假定语音信号符合拉普拉斯模型，并利用最大化似然函数方法估计混合矩阵。Zhou[13]等人提出非线性投影及列屏蔽(NonlinearProjectionandColumnMasking,NPCM)方法，利用该思想设计出目标函数，并利用粒子群算法求其目标函数的最优解，作为混合矩阵的估计。文献[14]提出退化解混估计(DegenerateUnmixingEstimationTechnique,DUET)，假设源信号符合W-不相邻正交条件，则可通过时频屏蔽的方法估计混合矩阵。ShengliXie[3]等人提出基于魏格纳分布的混合矩阵估计，并提出一种单源域识别方法。文献[15]提出基于混合时频比(TimeFrequencyRatioOfMixtures,TIFROM)方法，对观测信号进行分段，依靠小段内的方差进行单源域的判断，而小段的均值作为混合矩阵某列的估计。

上述方法主要存在三个缺陷：

应用范围受到限制。ICA[6-8]仅能解决超定问题，当涉及本发明中的欠定情形时，ICA方法失效。势函数法[11]及DUET[14]法仅适用于观测数为2的盲识别问题，而现实中，观测数目往往是多个。文献[12]提出的LOST无法实现对源数的估计，而在欠定盲分离问题中，源数目往往是未知的。

耗费计算量较大。FOOBI(FourthOrderBlindIdentification)[9]引入张量分析及同时对角化方法，需多次进行特征值分解。文献[13]提出的NPCM可以自动估计源数，然而其观测矩阵的向量要一个接一个的估计，不能够同时估计。同时，粒子群算法为了得到较准确的结果，须保证一定的种群数目和进化代数，因此引入庞大的计算量。

估计精度不够，对噪声不鲁棒。由于语音信号不能完全符合拉普拉斯模型，LOST方法[12]的估计精度十分不稳定。TIFROM[15]以某一小段的均值作为最终估计，因而易受噪声的影响。

由以上分析可知，现有算法无法实现高精确度及高效率的兼备。针对源数大于观测数的欠定情形，本发明提出一种基于频谱校正及数据密度聚类的欠定盲辨识方法。该方法利用语音信号的结构特征[16,17]，利用短时傅立叶变换对语音信号进行稀疏表达，引入频谱校正技术[18]来对信号进行谐波提取，并提出一种单源域判别准则筛选出单源域，最终利用最近提出的数据密度聚类法[19]对单源域进行聚类。该方法兼具高精度、高效率及对噪声的强鲁棒性。

参考文献

[1]AguileraP,CrucesS,etal.BlindseparationofdependentsourceswithaBoundedComponentAnalysisdeflationaryalgorithm[J].IEEESignalProcessingLetters,2013,20(7):709-12。

[2]LiyanQ,YinC,XuH,etal.Blindseparationofspeechsourcesinmultichannelcompressedsensing[C].2013IEEEInternationalInstrumentationandMeasurementTechnologyConference(I2MTC),1771-1774

[3]XieS,YangL,YangJM,etal.Time-frequencyapproachtounderdeterminedblindsourceseparation[J].IEEETransactionsonNeuralNetworks&LearningSystems,2012,23(2):306-16.

[4]LiZ,YanX,TianZ,etal.Blindvibrationcomponentseparationandnonlinearfeatureextractionappliedtothenonstationaryvibrationsignalsforthegearboxmulti-faultdiagnosis[J].Measurement,2013,46(1):259-71。

[5]LiK,SunG,XiaoM.DualchanneldelayspeechdenoisingbasedonunderdeterminedBSS[C].2014IEEEInternationalConferenceonMechatronicsandControl(ICMC),689-692

[6]HyvRinenA,OjaE.Afastfixed-pointalgorithmforindependentcomponentanalysis[J].Neuralcomputation,1997,9(7):1483-92。

[7]HyvarinenA.Fastandrobustfixed-pointalgorithmsforindependentcomponentanalysis[J].IEEETransactionsonNeuralNetworks,1999,10(3):626-34。

[8]ShenH,KleinsteuberM,HperK.LocalconvergenceanalysisofFastICAandrelatedalgorithms[J].IEEETransactionsonNeuralNetworks,2008,19(6):1022-32。

[9]DeLathauwerL,CastaingJ,CardosoJF.Fourth-ordercumulant-basedblindidentificationofunderdeterminedmixtures[J].IEEETransactionsonSignalProcessing,2007,55(6):2965-73。

[10]GeorgievP,TheisF,CichockiA.Sparsecomponentanalysisandblindsourceseparationofunderdeterminedmixtures[J].IEEETransactionsonNeuralNetworks,2005,16(4):992-6。

[11]ZibulevskyBM.Underdeterminedblindsourceseparationusingsparserepresentations[J].SignalProcessing,2001,81(11):2353–62。

[12]PaulD,GradyO.TheLOSTalgorithm:findinglinesandseparatingspeechmixtures[J].EURASIPJournalonAdvancesinSignalProcessing,2008,2008(1):17。

[13]ZhouG,YangZ,XieS,etal.Mixingmatrixestimationfromsparsemixtureswithunknownnumberofsources[J].IEEETransactionsonNeuralNetworks,2011,22(2):211-21。

[14]YilmazO,RichardS.Blindseparationofspeechmixturesviatime-frequencymasking[J].IEEEtransactionsonSignalProcessing2004,52(7):1830-47。

[15]AbrardF,DevilleY.Atime–frequencyblindsignalseparationmethodapplicabletounderdeterminedmixturesofdependentsources[J].SignalProcessing,2005,85(7):1389-403。

[16]SiegelLJ,BesseyA.Voiced/unvoicedmixedexcitationclassificationofspeech[J].IEEETransactionsonAcousticsSpeech&SignalProcessing,1982,30(3):451-60。

[17]VaseghiSV.Advanceddigitalsignalprocessingandnoisereduction[M].JohnWiley&Sons,2008。

[18]ZhangF,GengZ,YuanW.ThealgorithmofinterpolatingwindowedFFTforharmonicanalysisofelectricpowersystem[J].IEEETransactionsonPowerDelivery,2001,16(2):160-4。

[19]RodriguezA,LaioA.Clusteringbyfastsearchandfindofdensitypeaks[J].Science,2014,344(6191):1492-6。

[20]Aissa-El-BeyA,Linh-TrungN,Abed-MeraimK,etal.Underdeterminedblindseparationofnondisjointsourcesinthetime-frequencydomain[J].IEEETransactionsonSignalProcessing,2007,55(3):897-907。

[21]HeZ,HeZ,CichockiA,etal.ImprovedFOCUSSmethodwithconjugategradientiterations[J].IEEETransactionsonSignalProcessing,2009,57(1):399-404。

发明内容

为克服现有技术的不足，本发明旨在实现：提高算法精度；显著提高算法的效率；使算法对噪声的鲁棒性更强。本发明采用的技术方案是，频谱校正及数据密度聚类法语音欠定盲识别方法，包括如下步骤：

Step1：对观测信号x_m(t),m＝1,…,M做加汉宁窗L点50％重叠的短时傅立叶变换STFT，得到观测频谱X_m(t,k)，t为时间，k为谱线位置，m为当前观测的下标，M为观测数目；

Step2：逐帧对STFT观测频谱做频谱校正，对于某一特定时间帧t₀，具体操作为：

1)收集时间帧t0的观测频谱X_m(t₀,k)的所有峰值的位置k_p，p为当前谱峰的下标；

计算峰值X_m(t₀,k_p)及其次峰的比值v_p，m为当前观测下标：

$>v_p=\frac{X_m(t_0,k_p)}{max\left\{\right|X_m(t_0,k_p-1)|,|X_m(t_0,k_p+1)\left|\right\}}.---\left(1\right)$>

则可得到中间参数u_p

u_p＝(2-v_p)/(1+v_p)(2)

2)根据参数u_p估计频率偏移

则频率估计为

$>{\hat{\omega }}_{m,p}=(k_p+{\hat{\delta }}_p)2\pi /L---\left(4\right)$>

3)幅值估计及相位则分别为

$>{\hat{d}}_{m,p}=2\pi {\hat{\delta }}_p(1-{\hat{\delta }}_p^2)\left|X_m(t_0,k_p)\right|/sin\left(\pi {\hat{\delta }}_p\right),---\left(5\right)$>

其中ang(·)代表取相角运算，得到较多的谐波参数集合

Step3：对于上述特定时间帧t₀，对所有谐波参数对进行模式提纯，具体操作为

1)将频率集合进行混合，为观测1～M的频率成分；按照升序排列，得到Q个簇，其中第q个簇为为当前簇的下标，Γ_q为该簇的元素个数，为当前簇中p′个成分，则合并后的频率为

$>{\bar{\omega }}_q=\frac{1}{{\Gamma }_q}\underset{p^{\prime }=1}{\overset{{\Gamma }_q}{\Sigma }}{\tilde{\omega }}_{q,p^{\prime }},---\left(7\right)$>

则可得到合并后的频率序列

2)找出中满足给定一个小阈值ε>0，对于每个m都存在且仅存在一个下标满足下式：

$>|{\hat{\omega }}_{m,{\bar{p}}_{m,q}}-{\bar{\omega }}_q|<ϵ,$>

则该频率对应一个候选有效模式

其中为观测1～M相应成分的幅值，观测1～M相应成分的相位；

3)将有效模式作为输入，按照下式识别中的SAS模式，并记为

将所有投影的组合情况考虑后，共有种情况，得到最终的判定准则为

C代表组合操作，r及l代表从1到M遍历的两个变量；

Step4：将上述Step2及Step3逐帧进行，收集所有时间帧得到的SAS模式组成单源域

Ω＝{z_i,i＝1,…,P}，P为单源域的模式数目，z_i为单源域中第i个成分；

Step5：对单源域中的SAS模式进行数据密度聚类，具体操作如下

1)计算任两个模式点的特征距离d_i,j

d_i,j＝||z_i-z_j||,1≤i,j≤P,i≠j；(8)

i,j为从1到P变化的变量；

2)计算每个模式点的局部数据密度ρ_i,i＝1,…,P如下式，

$>{\rho }_i=\underset{j}{\Sigma }\chi (d_{i,j}-d_c),$>

其中，χ(·)为一个阈值函数，d_i,j为当前模式与其他模式的距离其中参数，d_c为邻域截止距离，

$>\chi \left(t\right)=\left(\begin{array}{c}1,t\le 0\\ 0,t>0\end{array}\right);---\left(9\right)$>

3)将ρ_i,i＝1,…,P按照降序排列，得到下表序列{q_i,i＝1,…,P}，使得ρ_q1≥ρ_q2≥…≥ρ_qP；

4)按照下式计算每个模式点的特征距离：

$>{\delta }_i=\left(\begin{array}{cc}\underset{j\ge 2}{\max }\left(d_{q_i,q_j}\right),& i=1;\\ \underset{q_j,j<i}{\min }\left(d_{q_i,q_j}\right),& i\ge 1;\end{array}\right)$>

其中，将所求局部数据密度序列按照降序排列，得到一个下标序列{q_i,i＝1,…,P}；

5)计算每个模式点的局部数据密度及特征距离的乘积γ_i＝ρ_iδ_i,i＝1,…,P，则必有N个γ_i显著大于其他γ_i，即将γ_i降序排列，必存在一个下标集合满足

$>{\gamma }_{{\stackrel{·}{q}}_1}\ge {\gamma }_{{\stackrel{·}{q}}_2}\ge ...\ge {\gamma }_{{\stackrel{·}{q}}_N}>>{\gamma }_j,j\notin \Lambda ---\left(10\right)$>

因而，即为N个聚类中心，对于中心将其周围距离小于的模式进行统计平均，即得到最终的混合矩阵列向量估计对所有聚类中心进行此操作，最终得混合矩阵估计

频谱校正及数据密度聚类法语音欠定盲识别装置，包括模数转化器、DSP，多路观测信号x(t)经过模数转化器采样得到样本序列x(n)，以并行数字输入的形式进入DSP芯片，经过DSP芯片的内部算法处理，得到混合矩阵的估计；最后显示混合矩阵的估计值；

其中，DSP包括如下模块

(1)核心算法模块，用于完成信号的STFT变换，频谱校正，模式提纯，单源模式筛选得到信号的单源域；

(2)聚类模块，用于利用基于数据密度的聚类法对单源域进行聚类，得到源数目及混合矩阵的估计，并将结果输出至驱动和显示模块；

(3)驱动和显示模块。

本发明的特点及有益效果是：

本发明提出的基于频谱校正及数据密度聚类法的语音欠定盲识别方法，在实际应用中可以产生如下有益效果：

第一该分离装置针对语音信号的欠定盲分离系统可，高精度地估计其源信号数目和混合矩阵。

第二利用频谱校正法处理观测信号，使得算法核心得以快速实现，大大降低了后续步骤的计算量，提高时效性。

第三提出了模式提纯过程，大大减少了计算复杂度，且精准的找出了候选有效模式，更进一步提高了性能。

第四提出了单源模式的识别方法，进一步降低计算复杂度，使本方法更加稳健，应用范围更广。

第五该方法引入了最近提出的基于数据密度的聚类法，十分简洁高效，且具备高精度。

第六具备足够的抗噪声性能，复杂度要求低，实效性好。

附图说明：

图1：g(t)频谱图。

图2：源信号波形图。

图3：观测信号波形图。

图4：第71帧观测信号频谱图。

图5：数据参数散点图。

图6：本发明的源数目和混合矩阵估计流程。

图7：模式提纯图。

图8：估计混合矩阵及理想混合矩阵散点图。

图9：不同观测信噪比下抗噪性能。

图10：本发明的硬件实施图。

图11：DSP内部程序流图。

具体实施方式

本发明的提出的基于频谱校正及数据密度聚类法的语音欠定盲识别装置具备以下功能：

1、引入频谱校正，可大大提高算法精度；

2、提出模式逐步筛选过程，并引入先进的数据密度聚类算法，可显著提高算法的效率；

3、经过筛选后的模式为强单源成分，又经过频谱校正后，更加接近理想值，因而使算法对噪声的鲁棒性更强。

本发明采用的技术方案是：

1.欠定盲分离系统的盲辨识问题

按照混合系统的不同模型，盲分离可分为线性瞬时混合、线性卷积混合、非线性混合。线性瞬时混合是盲信号分离中结构形式最简单，也是最经典的混合模型形式，因此本发明采用线性瞬时混合系统模型，其模型为

x(t)＝Αs(t)+n(t)(11)

式中s(t)＝[s₁(t)s₂(t)…s_N(t)]^T为源信号矢量，N为源信号数目；x(t)＝[x₁(t)x₂(t)…x_M(t)]^T为观测信号矢量，M为观测信号数目；n(t)＝[n₁(t)n₂(t)…n_N(t)]^T为噪声矢量，A为M×N维未知混合矩阵。考虑理想无噪情况，将(1)式表示为更易理解的向量形式

$>x\left(t\right)=\underset{n=1}{\overset{N}{\Sigma }}{a}_n{s}_n\left(t\right)---\left(12\right)$>

其中t代表时刻，a_n代表混合矩阵的n个列向量。若M<N，则称上述系统为欠定盲分离系统。

对于欠定情形，盲分离问题可分为两个阶段：1)盲辨识阶段，利用观测信号估计混合矩阵；2)盲恢复阶段，利用盲识别阶段的结果恢复源信号。从上述分析可知，盲辨识阶段的性能直接影响着后续步骤的性能，从而影响整个盲分离过程的精度。本发明旨在提出高效、精确、强噪声鲁棒性的盲辨识方法，再结合文献^[20,21]所提出的盲恢复算法，就可得到完整的盲分离系统。

2.基于频谱校正与数据密度聚类的语音欠定盲识别

2.1.基于频谱校正的谐波提取

2.1.1.短时傅立叶变换谱泄漏及其对盲识别的影响

假定时频变换工具是理想的。若在某时频点(t₀,ω₀)下仅有一个源s_n活跃，则(12)式可变为：

$>\left(\begin{array}{c}{X}_1(t,\omega )\\ X_2(t,\omega )\\ .\\ .\\ .\\ X_M(t,\omega )\end{array}\right)=\left(\begin{array}{c}{a}_{1,n}\\ a_{2,n}\\ .\\ .\\ .\\ a_{M,n}\end{array}\right)c_0{e}^{{\mathrm{j\phi }}_0}\delta (t-t_0,\omega -{\omega }_0),---\left(13\right)$>

其中n＝1,…,N，c₀,φ₀为源信号s_n在时频点(t₀,ω₀)的幅度参数及相位参数，而δ(·)为冲击函数，t,ω为时间和频率的自变量。因此观测向量[X₁(t,ω)X₂(t,ω)…X_M(t,ω)]^T平行于混合矩阵第n个列向量a_n。因而，仅凭借一个时频点(t₀,ω₀)即可得到观测矩阵第n列a_n的估计。

然而在实际当中，常用时频分析工具与(13)中理想时频工具间总是存在偏差。其原因在于：

(1)由于不确定准则的限制，任何时频分析工具的时间分辨率与频率分辨率不能够同时高；

(2)语音混合信号中包含大量成分，这些成分之间会相互干扰，从而降低盲识别方法的性能；

(3)不同时频工具有不同的本质缺陷。例如魏格纳分布会引入交叉项干扰，而短时傅立叶变换由平移加窗的傅立叶变换实现，因而会受到傅立叶变换本质的频谱泄漏的影响。

本发明引入短时傅立叶变换作为稀疏表示，因而有必要对其频谱泄漏现象进行研究。为了便于理解，以例1为例介绍频谱泄漏现象及其影响。

例1：考虑一个包含三个谐波成分的信号g(t)＝cos(2πf₁t)+cos(2πf₂t)+cos(2πf₃t)，其中f₁＝152Hz，f₂＝2f₁＝304Hz，f₃＝3f₁＝456Hz，t为时间自变量。假定采样频率为f_s＝16000Hz，对信号做L＝512点的傅立叶变换可得其频谱|G(k)|，则频率分辨率Δf＝f_s/L＝31.25Hz。由频偏公式

δ＝f/Δf-[f/Δf],(14)

可求得δ₁＝-0.1360，δ₂＝2δ₁＝-0.2720，δ₃＝3δ₁＝-0.4080，其中[·]代表取整运算。其频谱见图1。

从图1中可以看出理想的谱线(如图中虚线所示)演变为一簇近邻谱线(如图中黑色实线所示)。如图，频偏越大，谱峰与理想谱线的偏离越大，从而频谱泄漏越严重。如文献[16,17]所述，语音信号由近似谐波的浊音部分与形状与噪声相似的清音部分组成。浊音部分占据语音信号能量的大部分，且主要集中在低频频段。

因此，若将语音信号中谐波成分的理想谱线位置估计出来，则短时傅立叶变换(STFT)表示将被大大简化，从而使盲识别过程更具高效性。

2.1.2.基于频谱校正的谐波提取

为了实现语音信号中谐波成分的理想谱线的估计，本发明引入频谱校正技术。该技术可以精确的估计公式(13)中列出的谐波的三个参数(频率ω₀，幅值c₀及相位)，用三个参数代替一簇谱线描述信号频率成分。

对第m路观测信号进行频谱校正的方法为：在短时傅立叶变换(STFT)某一特定时间帧t₀的谱图中，第p簇谱线(代表一个频率成分)内选取相邻最大的两根进行比值(即将主谱线的幅值除以旁边幅值最大的一根谱线的幅值)，将该比值记为v_m,p，然后根据v_m,p求取参数u_m,p，如式(15)所示

u_m,p＝(2-v_m,p)/(v_m,p+1).(15)

根据u_m,p求出频率偏差(频偏)估计利用进行频率、幅值和相位的校正，若当前谱线处于谱峰k_p谱线处，则校正后的频率为

$>{\hat{\omega }}_{m,p}=(k_p+{\hat{\delta }}_{m,p})2\pi /L,---\left(16\right)$>

相位和幅值分别为

$>{\hat{d}}_{m,p}=\left|X_m(t_0,k_p)\right|\left(2\pi {\hat{\delta }}_{m,p}\right)\times (1-{\hat{\delta }}_{m,p}^2)/sin\left(\pi {\hat{\delta }}_{m,p}\right),---\left(17\right)$>

其中|X_m(t₀,k_p)|当前FFT谱峰的幅值，为k＝k_p谱线处的观测相位。经过谱校正，我们可以从混合的语音信号中得到各个峰值的频率、幅值、相位的精确估计(p为谱线簇数)，即语音信号谐波成分的频率、幅值、相位3参数的精确估计。以例2为例来说明频谱校正谐波提取及后续步骤的过程。

例2考虑3×4的混合系统(即源信号数目为4，观测信号数目为3)，指定一个归一化的混合矩阵

$>A=\left(\begin{array}{cccc}0.9356& 0.6354& 0.2813& 0.2858\\ 0.2433& 0.2102& 0.4571& 0.9147\\ 0.2557& 0.7430& 0.8438& 0.2858\end{array}\right).---\left(18\right)$>

其源信号s₁(t)～s₄(t)选自TIMIT语音数据库，对应的采样频率为f_s＝16000Hz。其波形如图2所示，3路观测信号x₁(t)～x₃(t)波形如图3所示。对x₁(t)～x₃(t)做STFT变换(FFT点数L＝512，加重叠率50％的汉宁窗，时间窗总数为219)可得其STFT频谱X₁(t,ω)～X₃(t,ω)，t,ω为时间和频率的自变量。

由上述分析可知，对3路观测信号做频谱校正，可以得到语音信号中得到m观测谱峰p的频率、幅值、相位的精确估计。以例2的第71帧观测信号为例，其频谱图如图4。从x₁(t)可以得到20组谐波成分，即共有20组；从x₂(t)可以得到55组谐波成分，因而共有55组；从x₃(t)可以得到14组谐波成分，同理共有14组。且不同观测信号的频率成分并不完全重合。

2.2.基于频谱校正的模式提纯

2.2.1.频率合并

值得注意的是由于噪声干扰及频率成分之间的相互干扰，即使一个单源成分，不同观测的频率校正结果会有轻微的差异，因而有必要将其进行合并。

假定每个观测信号的频率分量的元素数为Υ₁，…，Υ_m。将混合并升序排列，可以得到若干组仅有很小差异的的簇。不失一般性，假定第q(q＝1,…,Q)簇频率为(Γ_q为第q簇元素的数目，p′为1到Γ_q变化的变量)，因而存在关系Υ₁+…+Υ_m＝Γ₁+…+Γ_Q，Γ_Q为第Q簇元素的数目。

对于第q个簇，其频率估计可合并为

$>{\bar{\omega }}_q=\frac{1}{{\Gamma }_q}\underset{p^{\prime }=1}{\overset{{\Gamma }_q}{\Sigma }}{\tilde{\omega }}_{q,p^{\prime }}.---\left(19\right)$>

仍以上述例2为例，对及中总共89组数据进行上述合并操作，得到66个合并后的频率

2.2.2.组建候选有效模式

尽管合并后得到了Q个合并频率，这Q个合并频率却不一定被所有观测包含，因为各个观测的频率成分本身就不重合。

然而对于BSS模型来说，只要观测矩阵A不为零，所有源信号的成分都应出现在所有观测中。因此，合并成分中，只有被所有观测同时包含的成分，才可能是来自源信号的成分，也可称为候选有效模式，通常这类成分的能量都较高。其余成分则很可能来自噪声或者频率间的相互干扰。因而，候选有效模式的判断方式为，给定一个小阈值ε，对于每个观测都存在且唯一存在一个下标满足

$>|{\hat{\omega }}_{m,{\bar{p}}_{m,q}}-{\bar{\omega }}_q|<ϵ,---\left(20\right)$>

则即可被当作候选有效成分。假定总共有个下标集合符合条件，则可组成个候选有效模式为从1到变化的变量。

其中为观测1～M相应成分的幅值，观测1～M相应成分的相位。

对于上述例2，经过上述步骤，得到的候选有效模式为5组。其幅值与相位见表一。

表一71帧候选有效模式及其校正参数

为第个成分在观测1到观测3中对应的幅值及相位。从表一中可以看出，66个频率成分中，仅有5组是候选的有效模式(即被所有观测包含的成分)，这5组成分的幅值均较大。从实验中可观测到，排除的61组成分的幅值均较小，说明候选有效模式的筛选大幅度的减小了模式数目，精确的排除了噪声成分及成分间的干扰而产生的伪峰。

2.2.3.单源成分的判定

当观测中某一成分仅来自一个源信号s_n时，BSS模型可简化为

$>\begin{array}{c}x=A[0,...,s_n,...,0]\\ =s_n{a}_n.\end{array}---\left(22\right)$>

因而，理论上观测信号向量x与混合矩阵第n列a_n平行。为了表示方便，将式(21)中的表示为

其中为观测1～M相应成分的幅值，观测1～M相应成分的相位。

由2.1.1.节分析可知，单源(single-active-source,SAS)模式可用来估计混合矩阵列向量，因为模式方向与激活的源所对应的混合矩阵列向量平行。因而在中，只有单源的模式才是真正的有效模式，而这些有效模式必包含在候选有效模式中。为了选择真正有效的单源模式(SAS模式)，需提出一个SAS模式判别准则。若为SAS模式，由于(22)式中的s_n仅为一个复数标量，在允许幅值不确定性的条件下，则应平行于实数向量a_n。因此的M个相位必须相等，或者相互差π(即其相位对应的单位向量指向同一方向，或部分指向相反的方向)。换句话说，任两个相位构成的单位向量相互投影，其值的绝对值等于，或者接近于1，即

其中，‘〈·〉’代表内积操作。其中，r及l代表从1到M遍历的两个变量。

将所有投影的组合情况考虑后(共有种情况)，可以得到最终的判定准则为

若(25)成立，则说明可被认定为SAS模式。其中，C代表组合操作，M代表观测数目

从式(23)～(25)中可以看出，本发明所提出的SAS模式判定准则十分简易，仅需要少数几次内积与平均的操作即可完成。又因为该判定准则的输入是个由频谱校正精确获得的候选有效模式，本身已得到大幅度的精简，因此有效的提高了本发明的性能及效率。

同时，该SAS模式判定准则对噪声的鲁棒性较好。其原因在于每个校正后的幅值都是理想谱峰(如图1中虚线圆点所示)的近似估计。由于理想的谱峰幅值强于其所属的簇的其他谱线值，显然校正后对噪声的鲁棒性得到提高。

在经过SAS模式判别准则后，我们可以进一步将模式数目减小到个。个模式如下：

为单源模式，及为在观测1～观测M的幅值及相位。由上述分析可知，单源模式为复数向量。为了便于后续步骤，需组建实数向量。由于相位的一致性，将作为参考相位，因而的实数向量可以由如下表达式得出

其中

m为从1到M变化的变量。同时，实数单源模式必须归一化以方便后续步骤的进行。

对于上述例2，从表1中可以看出，仅仅前两个模式，其最大相位差异较小，大约在1°，其余的模式，最大相位差异均较大。对表1中5个候选有效模式实施SAS模式判别准则后，发现仅前两列符合条件，与实际观察一致。与此同时，可以看出选择出的两个模式z₁和z₂正好为幅值较大的两列模式，因而验证了其噪声鲁棒性较高。

$>z_1=\left(\begin{array}{c}84.5548\\ 28.1076\\ 97.8573\end{array}\right),z_2=\left(\begin{array}{c}28.3142\\ 9.4284\\ 33.6466\end{array}\right)---\left(28\right)$>

将这两个模式进行归一化

$>z_1=\left(\begin{array}{c}0.6389\\ 0.2124\\ 0.7394\end{array}\right),z_2=\left(\begin{array}{c}0.6296\\ 0.2096\\ 0.7481\end{array}\right)---\left(29\right)$>

可以看出，这两个模式与混合矩阵A的第二列十分接近。

2.3.SAS域及基于数据密度聚类

2.3.1.组建SAS域

上述步骤仅仅针对一个时间帧。然而某一帧信号不足以恢复混合矩阵A，原因如下：

(1)语音信号为非平稳信号，一个时间帧无法描述整个信号的特点；

(2)不同SAS模式间存在偏差，可通过统计平均降低这些偏差；

(3)同一帧中很有可能不能同时出现所有源的SAS模式，因而需考虑所有时间帧才能得到完整的估计。

通过一帧接一帧的实施模式提纯过程，可以得到一个SAS域Ω＝{z_i,i＝1,2,…,P}。为了从该SAS域中得到混合矩阵A的估计，需引入高效的聚类算法，本发明引入基于数据密度的聚类算法以获得混合矩阵A的估计

2.3.1.基于数据密度的聚类法及混合矩阵估计

基于数据密度的聚类算法^[19]有以下3个优点：

(1)聚类的数目可以自动确定；

(2)同时拥有高性能及高效性；

(3)不对初始值或某个参数敏感。

其特点是更加强调聚类中数据密度的作用。这主要反映在2个假定中：

(1)聚类中心周围的点都具有更低的局部密度；

(2)聚类中心与其他高密度的点距离叫远。

基于数据密度的聚类算法中，数据由两个比较重要的参数来描述，即局部数据密度和特征距离。局部密度ρ_i是指在某数据点半径为d_c的小邻域内其他数据点的个数

$>{\rho }_i=\underset{j}{\Sigma }\chi (d_{i,j}-d_c),---\left(30\right)$>

其中‘χ(·)’为一个阈值函数，d_i,j为当前模式与其他模式的距离，i,j为数据的下标。局部数据密度反映了某个数据周围小邻域内数据密度的分布情况。特征距离δ_i反映该点与大局部密度的点的距离，为了便于特征距离的求解，将所求局部数据密度序列按照降序排列，可得到一个下标序列{q_i,i＝1,…,P}，P为SAS域数据总数

$>{\delta }_i=\left(\begin{array}{cc}\underset{j\le 2}{\max }\left(d_{q_i,q_j}\right),& i=1;\\ \underset{q_j,j<i}{\min }\left(d_{q_i,q_j}\right),& i\ge 1;\end{array}\right).---\left(31\right)$>

得到每个数据的参数对后，我们就可以选择参数都很大的数据作为数据的聚类中心。绘制数据关于这两个参数的散点图，则聚类中心就可以很直观的观察出来。选取明显远离其他点的几个数据点作为聚类中心，并在聚类中心周围以为半径的小邻域内做统计平均，即可得到混合矩阵的最终估计。

对于上述例2，数据的参数对散点图如图5，最终估计结果为

$>\hat{A}=\left(\begin{array}{cccc}0.9350& 0.6350& 0.2827& 0.2837\\ 0.2456& 0.2088& 0.4564& 0.9148\\ 0.2556& 0.7437& 0.8436& 0.2873\end{array}\right)---\left(32\right)$>

本发明利用恢复信噪比来衡量盲识别效果的好坏。恢复信噪比(RecoverySignalNoiseRatio,RecoverySNR)的表达式为：

$>SNR=10\lg [\underset{i=1}{\overset{M}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{a}_{i,j}^2/\underset{i=1}{\overset{M}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{(a_{i,j}-{\hat{a}}_{i,j})}^2]---\left(33\right)$>

将A与元素值代入(33)中，排除次序不确定性后可得SNR＝53.37dB，lg为10为底的对数运算。

2.本发明的基于频谱校正与数据密度盲识别流程总结

本发明基于频谱校正与数据密度盲识别流程总结如图6。

Step1：对观测信号x_m(t),m＝1,…,M做加汉宁窗L点50％重叠的STFT变换，得到观测频谱X_m(t,k)t为时间，k为谱线位置，m为当前观测的下标，M为观测数目；

Step2：逐帧对STFT观测频谱做频谱校正，对于某一特定时间帧t₀，具体操作为：

4)收集时间帧t₀的观测频谱X_m(t₀,k)的所有峰值的位置k_p，计算峰值X_m(t₀,k_p)(m为当前观测下标)及其次峰的比值v_p

$>v_p=\frac{X_m(t_0,k_p)}{max\left\{\right|X_m(t_0,k_p-1)|,|X_m(t_0,k_p+1)\left|\right\}}.---\left(1\right)$>

则可得到参数u_p

u_p＝(2-v_p)/(1+v_p)(2)

5)根据参数u_p估计频率偏移

则频率估计为(其中L为汉宁窗点数)

$>{\hat{\omega }}_{m,p}=(k_p+{\hat{\delta }}_p)2\pi /L.---\left(4\right)$>

6)幅值估计及相位估计则分别为

$>{\hat{d}}_{m,p}=2\pi {\hat{\delta }}_p(1-{\hat{\delta }}_p^2)\left|X_m(t_0,k_p)\right|/sin\left(\pi {\hat{\delta }}_p\right),---\left(5\right)$>

其中ang(·)代表取相角运算，得到较多的谐波参数集合

Step3：对于上述特定时间帧t₀，对所有谐波参数对进行模式提纯，具体操作为

1)将为观测1～M的频率成分)频率集合进行混合，按照升序排列，得到Q个簇，其中第q个簇为为当前簇的下标，Γ_q为该簇的元素个数，为当前簇中p′个成分，则合并后的频率为

$>{\bar{\omega }}_q=\frac{1}{{\Gamma }_q}\underset{p^{\prime }=1}{\overset{{\Gamma }_q}{\Sigma }}{\tilde{\omega }}_{q,p^{\prime }},---\left(7\right)$>

则可得到合并后的频率序列

2)找出中满足给定一个小阈值ε>0，对于每个m都存在且仅存在一个下标满足(20)，则该频率对应一个候选有效模式如式(21)；

3)将作为输入，按照式(25)识别中的SAS模式，并记为

Step4：将上述Step2及Step3逐帧进行，收集所有时间帧得到的SAS模式组成单源域Ω＝{z_i,i＝1,…,P}；

Step5：对单源域中的SAS模式进行数据密度聚类，具体操作如下

6)计算任两个模式点的特征距离d_i,j

d_i,j＝||z_i-z_j||,1≤i,j≤P,i≠j；(8)

i,j为从1到P变化的变量；

7)计算每个模式点的局部数据密度如(30)，其中参数d_c为邻域截至距离，χ(·)为阈值函数

$>\chi \left(t\right)=\left(\begin{array}{c}1,t\le 0\\ 0,t>0\end{array}\right);---\left(9\right)$>

8)将ρ_i,i＝1,…,P按照降序排列，得到下表序列{q_i,i＝1,…,P}，使得ρ_q1≥ρ_q2≥…≥ρ_qP；

9)按照式(31)计算每个模式点的特征距离；

10)计算每个模式点的局部数据密度及特征距离的乘积γ_i＝ρ_iδ_i,i＝1,…,P，则必有N个γ_i显著大于其他γ_i。即将γ_i降序排列，必存在一个下标集合满足

$>{\gamma }_{{\stackrel{·}{q}}_1}\ge {\gamma }_{{\stackrel{·}{q}}_2}\ge ...\ge {\gamma }_{{\stackrel{·}{q}}_N}>>{\gamma }_j,j\notin \Lambda .---\left(10\right)$>

因而，即为N个聚类中心。对于中心将其周围距离小于的模式进行统计平均，即可得到最终的混合矩阵列向量估计对所有聚类中心进行此操作，最终可得混合矩阵估计

3.实验

本章将选自标准语音信号数据库TIMIT(TexasInstrumentsandMassachusettsInstituteofTechnology)的语音信号用于验证本发明的盲识别方法对欠定盲分离问题中源数目估计及混合矩阵估计的有效性。

实验一：本实验为了验证模式提纯步骤的效果。考虑例2中同样的信号，为了探究SAS模式判别准则对负值的鲁棒性，混合矩阵A设为

$>A=\left(\begin{array}{cccc}0.9356& 0.6354& 0.2813& 0.2858\\ -0.2433& -0.2102& 0.4571& 0.9147\\ 0.2557& -0.7430& -0.8438& 0.2858\end{array}\right).---\left(34\right)$>

STFT傅立叶点数L＝512，加汉宁窗，50％重叠率。

其模式提纯过程散点图如图7。图7(a)为所有短时傅立叶变换时频点的归一化后的散点图，图7(b)为候选有效模式的散点图，图7(c)为SAS模式散点图。

从图7(a)中可以看出，直接从所有短时傅立叶变换时频点中无法估计出混合矩阵，因为其分布杂乱无章，几乎占据了整个单位球面；经过频谱校正及候选有效模式的筛选后，图7(b)中模式数大大减小，而且可以清晰的看到4个簇，尽管存在一些明显的野点；经过SAS模式筛选后，模式数更进一步减小，且聚类倾向更为明显，大部分野点被去除。由上可知，模式提纯过程可大大减少模式数目，有效的提取出SAS模式，可大大提升盲识别效率。

最终，经过聚类后，估计出的观测矩阵的为(允许次序不确定性)

$>\hat{A}=\left(\begin{array}{cccc}0.9356& 0.6354& 0.2813& 0.2858\\ -0.2433& -0.2102& 0.4571& 0.9147\\ 0.2557& -0.7430& -0.8438& 0.2858\end{array}\right)---\left(35\right)$>

估计混合矩阵(‘*’)与原始混合矩阵(‘ο’)的散点图见图8

排除次序不确定性后，恢复信噪比为SNR＝55.25dB。

实验二：本实验旨在验证本发明盲识别算法对不同系统的稳定性。本实验将本发明的盲提取方法与另外三种已有的方法TIFROM^[15]、LOST^[12]、NPCM^[13]进行了对比。为了验证对系统的鲁棒性，本实验将盲提取过程进行了100次，每次的归一化混合矩阵A都是随机生成的。但考虑到任何算法都很可能因为混合矩阵的奇异而性能急剧下降，混合矩阵A必须保证条件数良好。为了发挥所有算法最好的性能，对随机生成的归一化混合矩阵A作如下规定：

(1)混合矩阵A中任两个列向量的夹角不得低于10°；

(2)混合矩阵A中所有元素绝对值均需大于0.1。

四种方法均做傅立叶点数L＝512，50％重叠率的汉宁窗的短时傅立叶变换。对于NPCM，其参数设置为衰减参数ρ＝10⁴，种群数P＝30，衍化代数K＝30，终止参数ε_T＝0.4。对于LOST，其源数固定为真实源数，拉普拉斯密度参数β初始化为1。对于TIFROM，STFT窗长设置为10。

本发明引入恢复信噪比来衡量盲识别的性能。因而，100次实验的平均恢复信噪比用来衡量方法的性能，用100次实验平均消耗时间来衡量方法的效率。不同方法的平均恢复信噪比及平均时间消耗见表二。

实验三：为了验证本发明的盲识别算法对噪声的鲁棒性，本实验对不同观测信号信噪条件比下的4种盲识别算法分别进行蒙特卡罗实验，每个信噪比下进行100次，求取其平均恢复信噪比作为对噪声鲁棒性能的衡量。其结果见图9。其中本发明盲识别方法标记为NPCM方法标记为‘×’，LOST标记为‘*’，TIFROM标记为‘ο’。

表二不同混合矩阵的鲁棒性比较

下面对实施本发明的硬件予以简单说明。

硬件实施图如图10，将采集到的多路观测信号x(t)经过A/D(模数转化器)采样得到样本序列x(n)，以并行数字输入的形式进入DSP芯片，经过DSP芯片的内部算法处理，得到混合矩阵的估计；最后借助输出驱动及其显示模块显示混合矩阵的估计值。

其中，图10的DSP(DigitalSignalProcessor，数字信号处理器)为核心器件，在信号参数估计过程中，完成如下主要功能：

(1)调用核心算法，完成信号的STFT变换，频谱校正，模式提纯，单源模式筛选得到信号的单源域；

(2)利用基于数据密度的聚类法对单源域进行聚类，得到源数目及混合矩阵的估计，并将结果输出至驱动和显示模块；

DSP器件的内部程序流程如图9所示。本发明将所提出的“一种稳健的欠定盲分离源数及混合矩阵估计方法”这一核心估计算法植入DSP器件内，基于此完成高精度、低复杂度、高效的源信号数目及混合矩阵的估计。

图11流程分为如下几个步骤：

1)首先需根据具体应用要求(如医学和军事等的具体测量要求)，设置信号的采样点数L；

2)其次，CPU主控器从I/O端口读采样数据，进入内部RAM；

3)最终，按图6本发明的处理过程进行混合矩阵估计，并将其通过外部显示装置进行显示。