CN105355212A - 一种稳健的欠定盲分离源数及混合矩阵估计方法及装置 - Google Patents

Abstract

本发明公开了一种稳健的欠定盲分离源数及混合矩阵估计方法及装置，方法包括：对观测频谱进行频谱校正；利用谱校正结果构造M×Q维的谐波参数矩阵，并对构造的矩阵进行归一化；对归一化后的矩阵建立距离矩阵，判断距离矩阵的每一列中元素是否小于阈值ξ，元素都不小于阈值的列对应的频率为重叠频率，将该频率对应的列从归一化后的矩阵中剔除；对处理后的矩阵计算势函数，并搜索势函数的峰值，确定源数目估计和混合矩阵估计。装置包括：输入模块、处理模块和输出模块。本发明针对周期平稳信号的欠定盲分离系统，可高精度地估计源信号数目和混合矩阵；利用频谱校正法处理观测信号，使得势函数求取得以快速实现，降低了后续计算量，提高了时效性。

Description

一种稳健的欠定盲分离源数及混合矩阵估计方法及装置

技术领域

本发明涉及数字信号处理技术领域，尤其涉及一种稳健的欠定盲分离源数及混合矩阵估计方法及装置，具体涉及对于在源信号间存在频率重叠时，周期性观测信号数目小于源信号数目(即欠定情况)的盲信号分离场合，仅依据观测信号进行处理而估计出源的数目以及混合矩阵。

背景技术

盲源分离(BlindSourceSeparation，BSS)是指在源信号及混合信道均未知的情况下，只利用传感器的观测信号来恢复源信号的过程^[1]。该问题在信号处理领域应用广泛。依据观测数目M与源信号的数目N的关系，可分为超定(M>N)、正定(M＝N)及欠定(M<N)三种情况。独立分量分析(IndependentComponentAnalysis，ICA)^[2]及改进的FastICA^[3]可处理超定情形，直接求出解混矩阵，即混合矩阵的伪逆。

为了实现欠定情形的源数估计和混合矩阵估计，现有的方法基本分3步：1)选择某个变换域对观测信号做稀疏表示；2)对观测信号的稀疏表示进行某种映射；3)对映射结果做聚类，估计源数和混合矩阵；例如2001年PauBofill及MichaelZibulevsky提出一种基于稀疏表达的两步法^[4]，第一步利用FFT作为稀疏表达变换，利用势函数法估计混合矩阵；2005年PandoGeorgiev等人研究了欠定盲分离的欠定程度与稀疏程度的关系，同时提出了一种超平面上的聚类法估计混合矩阵^[5]；2007年李宁，史铁林在传感器数与源数关系不明确源信号相互独立情况下，通过混合信号功率谱密度函数比值求解混合阵，并通过混合阵判定观测信号是完备混合、超定混合还是欠定混合^[6]；2008年谭北海与谢胜利提出一种基于源数估计的欠定盲分离算法，通过FFT来对信号做稀疏表示，并将观测信号的稀疏表示结果归一化到单位圆周上，利用弧线距离的统计分布的极值估计源数目及混合矩阵^[7]。2011年GuoxuZhou等人提出一种非线性投影列屏蔽方法，可估计源数未知的欠定盲分离问题的混合矩阵，同时定义了方向稀疏性，来衡量观测信号在某方向的稀疏程度^[8]。2012年ShengliXie等人利用时频分析方法魏格纳分布作为稀疏表达工具，并用K-均值方法估计混合矩阵的列向量^[9]。2012年毕晓君，宫汝江提出基于人工蜂群算法和K-均值的混合聚类方法对信号数据进行聚类，利用网格密度法修正每一类的聚类中心的混合矩阵估计法^[10]。2014年付卫红等人利用STFT(ShorttimeFourierTransform,短时傅立叶变换)对信号进行稀疏表达，并将K-均值法进行改进，使混合矩阵的估计更加精确^[11]。

以上各种估计方法，存在两个缺陷：

1)耗费的计算量很大。这是因为，以上方法在对观测信号的稀疏表示结果做映射和聚类时，需遍历所有的稀疏表示结果；例如文献^[4-8]均利用FFT作为稀疏表达的方式，显然FFT固有的频谱泄漏及栅栏效应等会导致混合结果的误差。同时，根据FFT的特性，样本数目足够大时才能精确估计出混合矩阵，样本数目愈多，需遍历的FFT谱线数目越多，因此计算量较大。文献^[9,11]利用时频工具(魏格纳分布和STFT)来做稀疏表达和K-均值聚类，时频工具是二维的稀疏表示手段，遍历二维数据会比一维情况引入更大的数据量，计算用时长。

2)对源信号成分的分布要求过于苛刻，应用范围窄。以上方法中有一个共同的前提，就是要求各路源信号的频率在所选定的稀疏表示域上，不能有重叠。当存在重叠频率时，会导致源数估计错误或混合矩阵列向量的方向偏移。最终导致盲分离失败。

参考文献

[1]AguileraP,CrucesS,Duran-DiazI,etal.BlindSeparationofDependentSourcesWithaBoundedComponentAnalysisDeflationaryAlgorithm[J].SignalProcessingLetters,IEEE,2013,(7):709-12.

[2]ComonP.Independentcomponentanalysis,Anewconcept？[J].SIGNALPROCESSING,1994,36(3):287–314.

[3]贾银洁,许鹏飞.基于FastICA的混合音频信号盲分离[J].信息与电子工程,2009,7(4):321-5.

[4]BofillP,ZibulevskyM.Underdeterminedblindsourceseparationusingsparserepresentations[J].SIGNALPROCESSING,2001,81(11):2353–62.

[5]PG,FT,AC.Sparsecomponentanalysisandblindsourceseparationofunderdeterminedmixtures[J].NeuralNetworks,IEEETransactionson,2005,16(4):992-6.

[6]李宁,史铁林.基于功率谱密度的盲信号分离[J].振动工程学报,2007,20(3):255-9.

[7]谭北海,谢胜利.基于源信号数目估计的欠定盲分离[J].电子与信息学报,2008,30(863-7).

[8]ZhouG,YangZ,XieS,etal.MixingMatrixEstimationFromSparseMixturesWithUnknownNumberOfSources[J].NeuralNetworks,IEEETransactionson,2011,22(2):211-21.

[9]XieS,YangL,YangJ-M,etal.Time-FrequencyApproachtoUnderdeterminedBlindSourceSeparation[J].IEEETransactionsonNeuralNetworksandLearningSystems,2012,23(2):306-16.

[10]毕晓君,宫汝江.基于混合聚类和网格密度的欠定盲矩阵估计[J].系统工程与电子技术,2012,34(3):614-8.

[11]付卫红,马丽芬,李爱丽.基于改进K-均值聚类的欠定混合矩阵盲估计[J].系统工程与电子技术,2014,(11):2143-8.

发明内容

本发明提供了一种稳健的欠定盲分离源数及混合矩阵估计方法及装置，本发明利用频谱校正法处理观测信号，降低了计算量；并且通过去除叠频率法扩大了应用范围，详见下文描述：

一种稳健的欠定盲分离源数及混合矩阵估计方法，所述方法包括以下步骤：

对观测频谱进行频谱校正；利用谱校正结果构造M×Q维的谐波参数矩阵，并对构造的矩阵进行归一化；

对归一化后的矩阵建立距离矩阵，判断距离矩阵的每一列中元素是否小于阈值ξ，元素都不小于阈值的列对应的频率为重叠频率，将该频率对应的列从归一化后的矩阵中剔除；

对处理后的矩阵计算势函数，并搜索势函数的峰值，确定源数目估计和混合矩阵估计。

所述方法还包括：对观测信号x_m(t),m＝1～M做加汉宁窗L点FFT变换，得到观测频谱X_m(k)。

所述对观测频谱进行频谱校正的步骤具体为：

逐个扫描观测频谱X_m(k)的Q_m簇谱线，逐个记录每簇谱线的最高谱线的位置；

求出谱线峰值与其左右侧相邻两谱线的最大值的比值；

求各簇谱线的频偏估计、频率校正值、幅值校正值和相位校正值。

所述利用谱校正结果构造M×Q维的谐波参数矩阵，并对构造的矩阵进行归一化的步骤具体为：

对各路观测幅值和相位的谱校正结果进行复数组合；

对谱校正得到的所有频率估计值的集合按从小到大排序、间距较近的多个频率取均值做合并，并聚类生成包含Q个频率值的向量F；

确定每个校正后的频率向量中的元素与向量F内某个子类的所属关系；

根据所属关系，把复幅度矢量中的每个复数据归入向量F中相应的子类中，构造谐波参数矩阵，并进行归一化处理。

所述对处理后的矩阵计算势函数，并搜索势函数的峰值，确定源数目估计和混合矩阵估计的步骤具体为：

取出处理后矩阵的实数部分，构造出M×Q'维实矩阵；

求出第q列的方向及模值；扫描角度，求角度上的势函数；

搜索势函数的峰值，得到峰值集合，峰值集合的元素个数即为源数目的估计。

一种稳健的欠定盲分离源数及混合矩阵估计装置，所述装置包括：

输入模块，用于对多路采集到的观测信号进行采样，以并行数字输入的形式进入处理模块；

处理模块，用于对观测频谱进行频谱校正；利用谱校正结果构造M×Q维的谐波参数矩阵，并对构造的矩阵进行归一化；对归一化后的矩阵建立距离矩阵，判断距离矩阵的每一列中元素是否小于阈值ξ，元素都不小于阈值的列对应的频率为重叠频率，将该频率对应的列从归一化后的矩阵中剔除；对处理后的矩阵计算势函数，并搜索势函数的峰值，确定源数目估计和混合矩阵估计；

输出模块，用于输出源数目估计和混合矩阵估计。

所述输入模块具体为：模数转换器。

所述处理模块具体为：DSP芯片。

本发明提出的稳健的欠定盲分离源数及混合矩阵估计，在实际应用中可以产生如下有益效果：

第一、本发明针对周期平稳信号的欠定盲分离系统，可高精度地估计其源信号数目和混合矩阵；

第二、本发明利用频谱校正法处理观测信号，使得势函数求取得以快速实现，大大降低了后续步骤的计算量，提高了时效性；

第三、本发明允许源信号之间存在轻度频率重叠，并针对源信号稀疏度的降低，提出了去除重叠频率，使本方法更加稳健，应用范围更广；

第四、本方法不仅可以处理理想的正弦组合信号，还可以处理近似周期信号的语音信号，具有很强的实用性；

第五、欠定盲分离具备足够的抗噪声性能，复杂度要求低，实效性要求高。

附图说明

图1为谐波源信号的时域波形及频谱图；

图2为谐波观测信号时域波形及频谱图；

图3为谐波信号势函数图；

图4为剔除重叠频率前的简化的势函数图；

图5为剔除重叠频率后的势函数图；

图6为一种稳健的欠定盲分离源数及混合矩阵估计方法的流程图；

图7为语音源信号时域波形及频谱图；

图8为语音观测信号时域波形及频谱图；

图9为语音观测信号势函数图；

图10为一种稳健的欠定盲分离源数及混合矩阵估计装置的硬件结构图；

图11为DSP内部流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步地详细描述。

实施例1

为了降低盲分离的计算复杂度，并且拓展盲信号分离的应用范围，本发明实施例提出一种稳健的欠定盲分离源数及混合矩阵估计方法，该方法针对周期或者近似周期的平稳信号，允许不同源信号在一定范围存在频率重叠，降低了对源信号稀疏度的要求，该方法包括以下步骤：

101：对观测频谱进行频谱校正；利用谱校正结果构造M×Q维的谐波参数矩阵，并对构造的矩阵进行归一化；

102：对归一化后的矩阵建立距离矩阵，判断距离矩阵的每一列中元素是否小于阈值ξ，元素都不小于阈值的列对应的频率为重叠频率，将该频率对应的列从归一化后的矩阵中剔除；

103：对处理后的矩阵计算势函数，并搜索势函数的峰值，确定源数目估计和混合矩阵估计。

综上所述，本发明实施例通过上述步骤实现了利用频谱校正法处理观测信号，降低了计算量；并且通过去除叠频率法扩大了应用范围，满足了实际应用中的多种需要。

实施例2

下面结合具体的计算公式、附图、例子对实施例1中的方案进行详细的描述，详见下文：

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

x(t)＝As(t)+n(t)(1)

其中，s(t)＝[s₁(t),s₂(t)...s_n(t)...s_N(t)]^T为源信号矢量，N为源信号数目；x(t)＝[x₁(t),x₂(t)...x_m(t)...x_M(t)]^T为观测信号矢量，M为观测信号数目；n(t)＝[n₁(t),n₂(t)...n_n(t)...n_N(t)]^T为噪声，A为M×N维未知混合矩阵。s_n(t)为第n个源信号(1≤n≤N)，n_n(t)为对应源信号的噪声，x_m(t)为第m个观测信号(1≤m≤M)。

考虑理想无噪情况，将(1)式表示为更易理解的向量形式：

其中，t代表时刻，a_n代表混合矩阵的n个列向量；s_n(t)为第n个源信号(1≤n≤N)。若M<N，则称上述混合系统为欠定盲分离系统。

本发明实施例针对源信号为周期平稳信号的欠定盲分离系统(M<N情况)，做稳健的源信号数目和混合矩阵估计。为方便理解，下面给出一个源信号为谐波组合的欠定盲分离系统的例子，详细说明欠定盲分离系统的源数和混合矩阵的估计过程。

例1：假定源数目N＝3，观测数目M＝2，采样速率f_s＝8000Hz，采样点数为L＝4096，令混合矩阵的方向为20°、80°、140°，则混合矩阵A为

这N＝3路源信号共包含Q＝10个频率，表1给出了各路源信号频率成分的幅值和相位参数表。各源信号的表达式为

其中，对于第q个频率成分(1≤q≤Q)f_q，a_n,q及分别为其在第n个源信号s_n(t)中对应的谐波成分的幅值及相位。

表1各路源信号频率成分的幅值和相位参数表

从表1可看出，整体上，s₁(t)～s₃(t)的频率成分分布是稀疏的，但却不是严格稀疏的。具体说来，s₁(t)与s₂(t)存在1个重合频率f_q＝250Hz,q＝4；s₂(t)和s₃(t)存在1个重合频率f_q＝300Hz,q＝5。根据表1可得如图1所示的源信号时域波形及其FFT频谱。

从图1可看出，源信号在时域内分布不是稀疏的，但经FFT转换到频域后，各源信号谱线分布呈现稀疏性。

假定观测信号采样率为f_s每路观测信号采集L个样点，则FFT频率分辨率△f＝f_s/L。则将线性混合欠定盲分离模型进行快速傅立叶变换后，(2)式变为

其中，X(k)为观测信号矢量x(t)的FFT，k为各谱线的位置；S_n(k)为第n(1≤n≤N)路源信号s_n(t)的FFT。由于总体上，源信号s₁(t)～s_N(t)的成分在FFT域中是稀疏分布的，也就是说，如图1所示，其FFT谱S₁(k)～S_N(k)的分布必然呈现少数几个明显大于零的稀疏分布特性，经过混合后，这种稀疏特征可以从各路观测结果后向逆推提取出来，即实现源数估计和混合矩阵估计。

仍以上述谐波信号为例，根据公式(1)对信号加均值为零，方差为0.01的高斯噪声，按照式(2)可得观测信号，其时域波形x₁(t)，x₂(t)及频谱图X₁(k)，X₂(k)如图2所示。从M＝2路观测波形x₁(t)～x₂(t)中估计源个数N和2×3维的混合矩阵A。

由于FFT存在固有的栅栏效应，信号进行FFT变换后，会产生频谱泄漏(即单个频率产生一簇谱线)，导致观测到的各路FFT峰值谱处对应的频率、幅值和相位值与真实值存在偏差，因而需要对各观测信号的FFT结果进行频谱校正。为提高校正精度，需要对观测信号加汉宁窗后再做FFT。

对观测信号进行频谱校正的方法为：在归一化后的加汉宁窗的FFT振幅谱线中选取相邻最大的两根进行比值(即将主谱线的幅值除以旁边幅值最大的一根谱线的幅值)，将该比值记为v，然后根据v求取比例偏差因子△k，如式(5)所示：

△k＝(2-v)/(v+1)(5)

根据△k进行频率、幅值和相位的校正，若当前谱线处于k＝m谱线处，则校正后的频率为：

其中，m为某谱峰的位置。

相位和幅值分别为

其中，|X(m)|为当前FFT谱峰的幅值，为k＝m谱线处的观测相位。

对于此例图2所示的2路观测信号，对其FFT频谱X₁(k),X₂(k)做上述频谱校正，所得结果见表2。信号进行频谱校正后，消除了频谱的泄漏，频率成分的位置从整数表示m△f，变为小数表示且相位非常接近表1的理想值。

另外，经过校正后，对于每路观测信号，仅仅需用表2的10列校正后的30个参数来描述观测信号。否则，若用图2所示的整幅FFT谱图来描述信号，则需L＝4096个复数据，故大大节省了内存空间以及后续步骤的计算量。观测信号的重叠频率是指被两个或者以上的源信号共有的频率。当信号中存在这种频率成分时，由于经过频谱校正后的信号数据量大大减少，利用各频率成分的估计的混合矩阵会产生较大误差。为剔除这些频率，需对重叠频率误差产生机理进行分析。

表2谐波观测信号频谱校正

其中，为第1路观测信号的第q个频率成分所对应的幅值，为其对应的相位。及同理。

对于源信号，假设其数目为N，按照谐波模型，假设所有信号的频率成分个数为Q，则s₁(t)～s_N(t)可表示为：

其中，f_n,q为第n路(1≤n≤N)源信号中第q(1≤q≤Q)个频率分量；c_n,q及为f_n,q对应的频率成分的幅值和相位。

由此可得出，源信号的N×Q谐波参数矩阵为：

其中，c_N,Q为第N路源信号的第Q个频率分量所对应的幅值。

源信号基波成分向量为：

其中，f_q为所有源信号第q个频率成分，Q为频率分量的个数。

因此，源信号矢量可表示为

s(t)＝[Cf(t)+conj(Cf(t))]/2(11)

其中，C为源信号谐波参数矩阵；conj(·)为求共轭运算(实部不变，虚部变为原来的相反数)。

根据式(1)，线性瞬时混合盲分离的表达式则变为：

x(t)＝Α[Cf(t)+conj(Cf(t))]/2(12)

同理，按照谐波模型，将观测信号表示为：

其中，为第M路观测，第q个频率成分的频率，q为1到Q间的任意数字；为对应频率成分的幅值；为对应频率成分的相位。可知

其中，为从观测估计的基波成分向量；β为观测信号的M×Q的谐波参数矩阵。

其中，为观测信号的谐波参数的第M行，Q列的元素。

显然，观测信号的基波成分向量为

考虑到共轭部分不含有新的信息，故只需考虑式(12)前半部分，不妨令变量 $\tilde{x}\left(t\right)=AC\hat{f}\left(t\right),$ 则变量符合 $x\left(t\right)=\tilde{x}\left(t\right)+conj\&lsqb;\tilde{x}\left(t\right)\&rsqb;.$ 因而问题转化为，如何根据观测的M×Q维的谐波参数矩阵β逆推出M×N维的混合矩阵A。很显然，这里源信号N×Q谐波参数矩阵C的稀疏特性，对该逆推过程有很大影响，因而需分情况讨论。

1)当系数矩阵C完全稀疏时(即C的任1列只有1个元素为非零值)：

该情况即对应着各路源信号的频率成分互不重叠的情况。此时，不妨假设C第n(n＝1,2,...,N)行第q(q＝1,2,...,Q)列的成分c_n,q非零，其他成分为零，则有

x′_q(t)＝[a₁...a_n...a_N][0...c_n,q...]^Tf_q(t)

(17)

＝c_n,qa_nf_q(t)

其中，f_q(t)为源信号基波成分向量f(t)第q行，x′_q(t)为观测信号的第q个基波成分所构成的信号成分。

因此联立(14)可得

由于从观测信号中估计的频率十分精确，可以认为则

β_q＝c_n,qa_n(19)

即观测矩阵的谐波参数矩阵的某列向量β_q的方向一定与混合矩阵的某列向量a_n的方向平行。由于盲分离允许存在幅值不确定性，这时仅由频谱校正得到的观测的M×Q维的谐波参数矩阵β即可完全确定混合矩阵A。

2)当系数矩阵C不满足完全稀疏条件时(即C的某些1列存在多个元素为非零值)：

这种情况即为本发明实施例针对的允许各路源信号的频率成分存在少部分重叠的情况。仍以第q(q＝1,2,...,Q)个分量为例，假设第n行及第l行均为非零值，即该成分在n及l个源中都存在，式(17)变为：

x'_q(t)＝[a₁...a_n...a_l...a_N][0...c_n,q...c_l,q...0]^Tf_q(t)

(20)

＝(c_n,qa_n+c_l,qa_l)f_q(t)

联立(14)可得

β_q＝c_n,qa_n+c_l,qa_l(21)

从该式可以看出，由于c_nq及c_lq的不确定性，不与任何混合矩阵列向量平行，故若把谐波参数矩阵β的这种列考虑进去就必然会导致估计混合矩阵A时产生误差。

因而，为解决因重叠频率导致的估计误差问题，就必须从M×Q维的谐波参数矩阵β中，把源信号中存在重叠的频率f_q剔除出来。从而根据删减部分列的β矩阵来估计源数和混合矩阵。

为寻找出存在重叠频率f_q的下标q，为减少源信号的幅值和相位对分类造成的影响，本发明实施例引入谐波参数矩阵β的归一化矩阵即对于每一列，做如下操作

其中，β_q为β的第q列；为的第q列。

进一步地，定义矩阵的各列之间的距离

其中，q₁,q₂为1～Q中任意两个不相等的数字；及为的不相同的两列。

显然，若距离足够小(假定小于某一阈值ξ)，则说明观测到的f_q1成分和f_q2成分相关性很大，来自同一个源的概率就很大。

进而可分两种情况(即某频率在各源信号中不重叠和重叠的两种情况)来讨论该频率的单位列向量与其他单位列向量的距离数值分布。

1)观测频率在各源信号中不重叠：

由式(19)、式(22)得知，该频率对应的单位列向量与混合矩阵的某一列a_n方向一致，因而必然是仅属于源信号s_n(t)的频率成分。再考虑任意一其他观测频率不外乎三种情况：

a)同一样，同属于源信号s_n(t)，且是一个不重叠的频率。

这时该频率对应的单位列向量也必然与源信号s_n(t)对应的混合矩阵的列向量a_n方向一致,所以一定是接近于0的值。

b)与同属于源信号s_n(t)，但却是一个与其他源重叠的频率，这时由式(21)可知，除了与a_n存在相关性外，还与至少另一个源对应的a_l存在相关性，这就使得的方向与混合矩阵的列向量a_n方向存在偏差，即方向与方向产生了偏差，所以 $d_{q_1,q_2}=\left|\right|{\widetilde{\mathrm{\&beta;}}}_{q_1}-{\widetilde{\mathrm{\&beta;}}}_{q_2}\left|\right|$ 是个较大的值。

c)f_q2同分别属于不同的信号源，这时方向与必然会产生很大的偏差,所以 $d_{q_1,q_2}=\left|\right|{\widetilde{\mathrm{\&beta;}}}_{q_1}-{\widetilde{\mathrm{\&beta;}}}_{q_2}\left|\right|$ 是一个很大的值。

2)观测频率在各源信号中是重叠的：

假定既属于源信号s_n(t)的成分，又属于源信号s_l(t)的成分。因而对应的单位列向量与混合矩阵对应的列a_n方向不一致。

从而再考虑任意一其他观测频率f_q2，则不管f_q2是不是属于源信号s_n(t)的成分,其距离都会是较大的值，即满足也就是说满足 $d_{q_1,q_2}=\left|\right|{\widetilde{\mathrm{\&beta;}}}_{q_1}-{\widetilde{\mathrm{\&beta;}}}_{q_2}\left|\right|<\mathrm{\&xi;}$ 的距离个数为0。

综上所述，判断某观测频率是不是重叠频率的依据，仅需判断的距离个数是不是为0就可以了，若为0，就应当把频率剔除出去。

对例1而言，根据表2的频谱校正的结果，按照式(15)可构造M×Q维的谐波参数矩阵β，进而根据式(22)可得归一化矩阵再根据式(23)可算出(Q-1)×(Q-1)＝81个距离值，如表3所示(假定阈值ξ＝0.1)。

可看出，表3的各列，除了f₄与f₅对应的两列小距离数目为0外，其余各列小距离数目都至少为1，故可判定检测到的两个频率成分是重叠的频率成分，这与表1列出的源成分数据是吻合的，故应把剔除出去。

另外，进一步还可从表3推出，应分别归属不同的信号源。

表3归一化矩阵的距离列表d_q1，q2(阈值ξ＝0.1)

这时谐波参数矩阵的维数由M×Q变为M×(Q-2),即为由(Q-2)个不重叠频率对应的谐波参数列向量组成的矩阵

如前所述，混合矩阵A每一列可由某个角度方向θ来描述，但仅从观测来看，具体落在那几个角度方向上是未知的。因而势函数法需对在[0,180°]范围内的所有角度方向θ进行扫描，构造每个角度方向θ上的势函数值，从而找出势函数的峰值个数，即可获得源数估计。找出势函数峰值所对应的角度方向θ，即可估计出混合矩阵A所有的列向量。

获得归一化谐波参数矩阵后，给定某个角度方向θ，势函数可由中所有列对角度方向θ的总贡献来决定。这里的贡献，包括角度贡献和幅值贡献两个方面。

对于角度贡献的定义，假定谐波参数矩阵中的某个列向量的角度为则可根据θ与的差值定义如下角度贡献：

式(24)中，θ与的角度差值越小时，即在扫描角度方向θ上的投影越大，则角度贡献必然越大(其贡献率可通过参数λ来调节)。反之，两者角度差异值越大，则角度贡献就越小，超过π/4时，则认为没有贡献。

对于幅值贡献的定义，列向量的2-范数即是1个合理的选择。故角度方向θ上的势函数可定义为

注意，在传统势函数求取中，由于没有做频谱校正提取谐波参数，式(25)中的统计平均只能在L个FFT谱线位置上进行(L×Q)。对例1求得的传统势函数如图3所示。

从图3可看出，由于存在频率重叠，传统势函数的峰值落在θ＝20°,θ＝46°,θ＝80°,θ＝108°,θ＝141°上，其峰值数目为5个(比实际N＝3多出2个)，与例1理想的源信号角度θ＝20°,θ＝80°,θ＝140°相比，势函数检测到的θ＝46°,θ＝108°为伪峰。而这个检测结果是L＝4096根FFT谱线做统计平均而得到的。

很显然，本发明实施例引入谱校正后，势函数求取耗费的L根FFT谱线统计平均可由少数几个谐波参数矩阵β的列向量平均所替代，故大大缩减了势函数的计算量。

对于例1，若把参数矩阵全部Q＝10个列向量都计入式(25)的统计平均计算，则可得到如图4所示的势函数。

从图4可看出，由于存在频率重叠，简化的势函数的峰值落在θ＝20°,θ＝45°,θ＝80°,θ＝107°4个位置(比实际N＝3多出1个)，与例1理想的源信号角度θ＝20°,θ＝80°,θ＝140°相比，势函数检测到的θ＝45°,θ＝107°为伪峰，剩下的θ＝20°,θ＝80°为正确检测到的峰值，漏检了峰值θ＝140°。也就是说，因为把两个重叠的频率计入式(25)的统计平均，对其他非重叠的频率产生了较大的干扰，造成了势函数Φ(θ)峰值的偏差和漏检。

因而，为提高势函数的估计精度，必须把重叠的频率剔除出去。从而，谐波参数矩阵由M×Q的β变为M×Q'(Q'<Q)的相应地，势函数的定义应改进为

其中，Q'为去掉重叠频率后谐波数目；为的第q列。

对于例1，Q'＝8,则根据式(26)计算而得到的势函数曲线如图5所示。

从图5可看出，剔除掉重叠频率后，简化的势函数的峰值落在期望的θ＝20°,θ＝80°,θ＝140°的3个位置，峰值个数正好为真实的N＝3个，既不增加，也不减小，而且峰值位置十分准确。体现了本发明实施例的稳健欠定盲分离源数及混合矩阵估计法的高精度、低计算量的优势。

由检测到的可获得估计矩阵为

与真实值完全一致。图6包含的5个具体步骤为：

201：对观测信号x_m(t),m＝1～M做加汉宁窗L点FFT变换，得到观测频谱X_m(k)；

202：对X_m(k)进行频谱校正，具体操作为：

1)逐个扫描X_m(k)的Q_m簇谱线(每簇谱线代表1个频率成分)，逐个记录每簇谱线的最高谱线的位置，将其记为 $k_{m,1},k_{m,2},...,k_{m,Q_m},m=1~M.$

2)求出上述谱线峰值X_m(k_m,q)，q＝1,2,...,Q_m与其左右侧相邻两谱线的最大值的比值v，即：

其中，X_m为第m路观测的FFT谱,k_m,q为X_m的第q个峰值谱线的位置。

3)求各簇谱线的频偏估计△k_m,q、频率校正值幅值校正值和相位校正值即

其中，为直接从第m路观测第q簇的峰值谱线处读取的相位值。

203：利用谱校正结果构造M×Q维谐波参数β矩阵，并做矩阵归一化；

1)对各路观测幅值和相位的谱校正结果进行复数组合：

其中，为第m路观测第q个基波频率成分的复幅度，包含了该成分的幅值及相位信息。从而得到M个复幅度矢量 ${\widehat{\mathrm{\&beta;}}}_m=\&lsqb;\begin{array}{cccc}{\mathrm{\&beta;}}_{m,1}& {\mathrm{\&beta;}}_{m,2}& ...& {\mathrm{\&beta;}}_{m,Q_m}\end{array}\&rsqb;,m=1,...,M.$

2)对谱校正得到所有频率估计值的集合按从小到大排序、间距较近的多个频率取均值做合并，而聚类生成包含Q个频率值的向量F：

F＝[F₁F₂,...,F_Q](33)

同时也确定每个校正后的频率向量中的元素与F＝[F₁F₂,...,F_Q]内某个子类的所属关系(Q_m<Q)。

3)根据上述所属关系，把复幅度矢量中的每个归入F＝[F₁F₂,...,F_Q]中相应子类中，进而设定一小阈值ε，根据下式

而构造谐波参数矩阵

4)将谐波参数矩阵β每列按照下式做归一化处理而构造矩阵

204：计算(Q-1)×(Q-1)个距离，并判别、剔除重叠频率。

按照下式计算(Q-1)×(Q-1)个距离，而实际构造一Q×Q维的距离矩阵D

设定一小值阈值ξ，统计距离矩阵D每一列中距离小于阈值ξ的数目，记录数目为0的列标号，把该列标号对应的频率认为是重叠频率，并且从矩阵中把该列标号对应的列剔除出去，从而M×Q维的归一化谐波参数矩阵缩减为M×Q'维的矩阵

205：计算势函数，并搜索峰值，确定源数目估计和混合矩阵估计。

1)取出的实数部分，构造出M×Q'维实矩阵

2)按照下式求出第q(q＝1,2,...,Q′)列的方向θ_q及模值l_q

θ_q＝tan^-1(B_2,q/B_1,q)q＝1,2,...,Q′(37)

3)扫描角度θ＝(1,2,...,180)×π/180，按照下式求角度θ上的势函数

其中，为对θ的角度贡献值，即

4)搜索势函数Φ(θ)的峰值，得到峰值集合其中集合的元素个数N即为源数目的估计。而混合矩阵估计则为：

综上所述，本发明实施例通过上述的计算公式、步骤实现了利用频谱校正法处理观测信号，降低了计算量；并且通过去除叠频率法扩大了应用范围，满足了实际应用中的多种需要。

实施例3

下面结合具体的例子、附图和计算公式对实施例1和2中的方案进可行性验证，详见下文描述：

本验证用以验证本发明遂犯对欠定混合盲源分离的源数目估计和混合矩阵的有效性，利用近似周期语音信号验证盲分离算法的可行性。

文献[4]中的六个长笛语音信号为六个相对稳定，近似周期，频率几乎不重叠的语音信号。利用六个长笛信号的线性叠加，构造复合信号，使信号间有部分频率重叠，s₁(t)～s₃(t)的构造方式如下

s₁(t)＝Sound₁(t)+Sound₄(t)

s₂(t)＝Sound₄(t)+Sound₅(t)+Sound₆(t)

s₃(t)＝Sound₂(t)+Sound₆(t)

其中，Sound₁(t)～Sound₆(t)为采样频率f_s＝8000Hz，采样点数为L＝32768的长笛语音信号。可以看出，s₁(t)与s₂(t)共有Sound₄(t)部分，s₂(t)与s₃(t)共有Sound₆(t)部分。s₁(t)～s₃(t)时域波形及其频谱图如图7。

令混合矩阵的方向为45°、105°、140°，得混合矩阵A为

根据公式(1)可得观测信号如图8。

下面根据观测信号对源数目及混合矩阵进行估计。

对观测信号进行频谱校正；再去除重合频率；将得到的单源频率作为势函数的输入，得到势函数图，如图9所示。

求出势函数的峰值为：45°105°140°，则混合矩阵估计为：

可见混合矩阵估计与真实混合矩阵估计A几乎完全一致。从而验证了本方法的正确。

实施例4

一种稳健的欠定盲分离源数及混合矩阵估计装置，该装置包括：

输入模块，用于对多路采集到的观测信号进行采样，以并行数字输入的形式进入处理模块；

处理模块，用于对观测频谱进行频谱校正；利用谱校正结果构造M×Q维的谐波参数矩阵，并对构造的矩阵进行归一化；对归一化后的矩阵建立距离矩阵，判断距离矩阵的每一列中元素是否小于阈值ξ，元素都不小于阈值的列对应的频率为重叠频率，将该频率对应的列从归一化后的矩阵中剔除；对处理后的矩阵计算势函数，并搜索势函数的峰值，确定源数目估计和混合矩阵估计；

输出模块，用于输出源数目估计和混合矩阵估计。

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

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

(1)调用核心算法，完成信号的幅值校正，得到正确的信号频率、幅值和相位特征；

(2)利用信号的幅值、相位特征去除重叠频率，绘制势函数图，并求其极值点，得到混合矩阵的估计，并将结果输出至驱动和显示模块；

DSP器件的内部程序流程如图11所示。本发明实施例将估计算法植入DSP器件内，基于此完成高精度、低复杂度、高效的源信号数目及混合矩阵的估计。

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

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

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

3)再次，对信号进行去直流处理，即令信号减去计算样点的平均值，消除待测信号中的直流成分的影响；

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

本发明实施例对各器件的型号除做特殊说明的以外，其他器件的型号不做限制，只要能完成上述功能的器件均可。

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种稳健的欠定盲分离源数及混合矩阵估计方法，其特征在于，所述方法包括以下步骤：

对观测频谱进行频谱校正；利用谱校正结果构造M×Q维的谐波参数矩阵，并对构造的矩阵进行归一化；

对归一化后的矩阵建立距离矩阵，判断距离矩阵的每一列中元素是否小于阈值ξ，元素都不小于阈值的列对应的频率为重叠频率，将该频率对应的列从归一化后的矩阵中剔除；

对处理后的矩阵计算势函数，并搜索势函数的峰值，确定源数目估计和混合矩阵估计。

2.根据权利要求1所述的一种稳健的欠定盲分离源数及混合矩阵估计方法，其特征在于，所述方法还包括：

对观测信号x_m(t),m＝1～M做加汉宁窗L点FFT变换，得到观测频谱X_m(k)。

3.根据权利要求1所述的一种稳健的欠定盲分离源数及混合矩阵估计方法，其特征在于，所述对观测频谱进行频谱校正的步骤具体为：

逐个扫描观测频谱X_m(k)的Q_m簇谱线，逐个记录每簇谱线的最高谱线的位置；

求出谱线峰值与其左右侧相邻两谱线的最大值的比值；

求各簇谱线的频偏估计、频率校正值、幅值校正值和相位校正值。

4.根据权利要求1所述的一种稳健的欠定盲分离源数及混合矩阵估计方法，其特征在于，所述利用谱校正结果构造M×Q维的谐波参数矩阵，并对构造的矩阵进行归一化的步骤具体为：

对各路观测幅值和相位的谱校正结果进行复数组合；

对谱校正得到的所有频率估计值的集合按从小到大排序、间距较近的多个频率取均值做合并，并聚类生成包含Q个频率值的向量F；

确定每个校正后的频率向量中的元素与向量F内某个子类的所属关系；

根据所属关系，把复幅度矢量中的每个复数据归入向量F中相应的子类中，构造谐波参数矩阵，并进行归一化处理。

5.根据权利要求1所述的一种稳健的欠定盲分离源数及混合矩阵估计方法，其特征在于，所述对处理后的矩阵计算势函数，并搜索势函数的峰值，确定源数目估计和混合矩阵估计的步骤具体为：

取出处理后矩阵的实数部分，构造出M×Q'维实矩阵；

求出第q列的方向及模值；扫描角度，求角度上的势函数；

搜索势函数的峰值，得到峰值集合，峰值集合的元素个数即为源数目的估计。

6.一种稳健的欠定盲分离源数及混合矩阵估计装置，其特征在于，所述装置包括：

输入模块，用于对多路采集到的观测信号进行采样，以并行数字输入的形式进入处理模块；

处理模块，用于对观测频谱进行频谱校正；利用谱校正结果构造M×Q维的谐波参数矩阵，并对构造的矩阵进行归一化；对归一化后的矩阵建立距离矩阵，判断距离矩阵的每一列中元素是否小于阈值ξ，元素都不小于阈值的列对应的频率为重叠频率，将该频率对应的列从归一化后的矩阵中剔除；对处理后的矩阵计算势函数，并搜索势函数的峰值，确定源数目估计和混合矩阵估计；

输出模块，用于输出源数目估计和混合矩阵估计。

7.根据权利要求6所述的一种稳健的欠定盲分离源数及混合矩阵估计装置，其特征在于，所述输入模块具体为：模数转换器。

8.根据权利要求6所述的一种稳健的欠定盲分离源数及混合矩阵估计装置，其特征在于，所述处理模块具体为：DSP芯片。