CN104217722A - 一种海豚哨声信号时频谱轮廓提取方法 - Google Patents

Abstract

本发明涉及一种海豚哨声信号时频谱轮廓提取方法。本发明包括：对真实的海豚哨声信号作短时傅里叶变换，获得海豚哨声时频信号；得到第m个数据块短时傅里叶变换结果的最大值e(m)、傅里叶变换后能量最大值e(m)对应的频率的采样点index(m)以及在index(m)/2处对应的能量值e(m)′，通过判决式：e(m)′＞βP₀确定第m个数据块对应的哨声信号基频频率；去除野点，得到平滑的基频轮廓；进行最小二乘多项式拟合，得到一个连续变化的光滑哨声信号时频谱轮廓线。由于海豚哨声信号时频谱轮廓曲线是进行海豚种群分类、识别的重要依据，本方法将对海豚的生物学研究提供帮助，仿照海豚哨声信号进行仿生通信，是实现水声隐蔽通信的有效方法。

Description

一种海豚哨声信号时频谱轮廓提取方法

技术领域

本发明涉及一种海豚哨声信号时频谱轮廓提取方法。

背景技术

海豚的听觉系统和感觉器官很发达，其声呐系统具有很强的回声定位和探测识别能力。为了交流，经过漫长的进化，海豚声呐具有通信的性能。海豚的通信本领十分发达，并且拥有着丰富的“语汇”。因此，对其声音数据的采集、监测、记录和分析具有重要的意义，可以更进一步了解它的声呐系统工作特性，为鲸豚动物仿生学提供新的灵感。

在海豚的社会交际中，主要产生两种声音：一种是哨声，海豚之间可以利用哨声保持联系，发出猎物坐标信息，以及捕猎的策略；另一种是突发脉冲，与哨声信号相比，突发脉冲更为复杂多变。突发脉冲的特点是重复次数很高，每组脉冲之间的间隔很短。研究发现突发脉冲和海豚的多种行为相联系，一般在警告、惊吓、恐吓、入侵或者两只海豚面对面张开嘴进行点头、摇头、弓起背等行为时会发出。

发明内容

本发明的目的在于提出一种海豚哨声信号时频谱轮廓提取方法。

本发明的目的是这样实现的：

(1)对真实的海豚哨声信号作短时傅里叶变换，获得海豚哨声时频信号X_m(k)；

(2)通过 $\left\{\begin{array}{c}[e\left(m\right),\mathrm{index}\left(m\right)]=\max \left|X_m\left(k\right)\right|\\ e{\left(m\right)}^{\′}=X_m\left(\frac{\mathrm{index}\left(m\right)}{2}\right)\end{array}\right.$ 得到第m个数据块短时傅里叶变换结果的最大值e(m)、傅里叶变换后能量最大值e(m)对应的频率的采样点index(m)以及在index(m)/2处对应的能量值e(m)′，通过判决式：e(m)′>βP₀确定第m个数据块对应的哨声信号基频频率，当判决成立时，提取的第m个数据块哨声信号基频频率采样点为index(m)/2，当判决不成立时，提取的第m个数据块哨声信号基频频率采样点为index(m)，其中，P₀表示噪声的能量，系数β的值大于等于2；

(3)对步骤(2)得到的短时傅里叶变化的基频频率值进行滤波，去除其中的野点，得到平滑的基频轮廓；

(4)对步骤(3)滤波后的基频曲线进行最小二乘多项式拟合，得到一个连续变化的光滑哨声信号时频谱轮廓线。

滤波方法为动窗加权平均法。

滤波方法为圆筛滤波法。

滤波方法为卡尔曼滤波法。

滤波方法为动窗加权平均法、圆筛滤波法、卡尔曼滤波法三种滤波法组合串行级联的方式。

本发明的有益效果在于：

由于海豚哨声信号时频谱轮廓曲线是进行海豚种群分类、识别的重要依据，因此，本方法将对海豚的生物学研究提供帮助。同时，利用本方法提取出的海豚哨声信号时频谱轮廓曲线，仿照海豚哨声信号进行仿生通信，是实现水声隐蔽通信的有效方法。

附图说明

图1为海豚哨声信号时频谱轮廓提取流程图；

图2为海豚哨声信号时频图；

图3为采用最大值提取法得到的基频频率值曲线图；

图4为判决式e(m)′>βP₀中β等于2时得到基频频率值曲线图；

图5为采用动窗加权平均滤波法得到的基频轮廓；

图6为采用圆筛滤波法得到的基频轮廓图；

图7为采用卡尔曼滤波法得到的基频轮廓图；

图8为采用圆筛滤波法和卡尔曼滤波法串行级联滤波得到的基频轮廓图；

图9为最小二乘多项式拟合处理后的结果和残差；

图10为基波时频谱轮廓与二次谐波时频谱轮廓。

具体实施方式

本发明涉及一种海豚哨声信号时频谱轮廓提取方法。

步骤1：对真实的海豚哨声信号作短时傅里叶变换，获得海豚哨声时频信号X_m(k)；

步骤2：通过表达式 $\left\{\begin{array}{c}[e\left(m\right),\mathrm{index}\left(m\right)]=\max \left|X_m\left(k\right)\right|\\ e{\left(m\right)}^{\′}=X_m\left(\frac{\mathrm{index}\left(m\right)}{2}\right)\end{array}\right.$ 得到第m个数据块短时傅里叶变换结果的最大值e(m)、傅里叶变换后能量最大值e(m)对应的频率的采样点值index(m)以及在index(m)/2处对应的能量值e(m)′，通过判决式：e(m)′>βP₀确定第m个数据块对应的哨声信号基频频率,当判决成立时，提取的第m个数据块对应的哨声信号基频频率采样点值为index(m)/2，当判决不成立时，提取的第m个数据块对应的哨声信号基频频率采样点值为index(m),其中，P₀表示噪声的能量，系数β的值大于等于2；

步骤3：对步骤2得到的短时傅里叶变化的基频频率值进行滤波，进一步去除其中的野点，得到平滑的基频轮廓；

步骤4：对滤波后的基频曲线进行最小二乘多项式拟合，得到一个连续变化的光滑哨声信号时频谱轮廓线。

由此可以得到一个连续变化的光滑哨声信号时频谱轮廓线。实现海豚哨声信号基波时频谱轮廓提取。

由得到的连续变化哨声信号基波时频谱轮廓曲线，可以依据谐波的倍数关系合成任意次谐波的时频轮廓曲线。

上述的海豚哨声信号时频谱轮廓提取方法，步骤3中使用的滤波方法可以是动窗加权平均法、圆筛滤波法、卡尔曼滤波法或动窗加权平均法、圆筛滤波法、卡尔曼滤波法三种滤波法组合串行级联的方式。使用上述滤波方法可以去除步骤2中得到的短时傅里叶变化的基频频率值中的野点，得到平滑的基频轮廓。

采用本发明所述海豚哨声信号时频谱轮廓提取方法，可以避免提取基频轮廓曲线时谐波造成的影响，从而得到正确的基频频点。同时，通过滤波、拟合可以去除野点，得到清晰、平滑、连续变化的哨声信号时频谱轮廓曲线。

由于采集到的真实海豚哨声信号通常信噪比较低，对短时傅里叶变换后的时频谱处理，提取信号时频谱轮廓线，通常会得到野点，所以结合海豚哨声信号具有谐波特点，应用判决式的判决结果提取谱轮廓线。提取出的轮廓线，若依然含有噪声，采用滤波的方法，实现平滑的效果，对处理后的结果采用最小二乘多项式拟合法，得到最终的哨声信号时频谱轮廓曲线。具体处理流程如图1所示。

步骤1：对真实的海豚哨声信号作短时傅里叶变换，获得海豚哨声时频信号，X_m[k]表示以窗长为L、第m个数据块短时傅里叶变换的结果，其中哨声信号时频图如图2所示；

步骤2：结合海豚哨声信号时频谱的特点，通过表达式

$\left\{\begin{array}{c}[e\left(m\right),\mathrm{index}\left(m\right)]=\max \left|X_m\left(k\right)\right|\\ e{\left(m\right)}^{\′}=X_m\left(\frac{\mathrm{index}\left(m\right)}{2}\right)\end{array}\right.---\left(1\right)$

得到第m个数据块短时傅里叶变换结果的最大值e(m)、傅里叶变换后能量最大值e(m)对应的频率的采样点值index(m)以及在index(m)/2处对应的能量值e(m)′。通过判决式：

e(m)′>βP₀           (2)

确定第m个数据块对应的哨声信号基频频率。当判决成立时，提取的第m个数据块对应的哨声信号基频频率采样点值为index(m)/2；当判决不成立时，提取的第m个数据块对应的哨声信号基频频率采样点值为index(m)。其中，P₀表示噪声的能量，系数β的值大于等于2。

上述的海豚哨声信号时频谱轮廓提取方法，步骤2的判决式中的参数β的值为无穷大。这时，判决e(m)′>βP₀恒不成立，认为能量值e(m)′是噪声。直接采用最大值提取法，利用公式：[e(m),index(m)]＝max|X_m(k)|，提取的第m个数据块对应的哨声信号基频频率采样点值为index(m)。

上述的海豚哨声信号时频谱轮廓提取方法，步骤2的判决式中的参数β等于2。当index(m)/2处对应的能量值e(m)′满足e(m)′>2P₀而大于噪声能量2倍时，可以认为能量值e(m)′不是噪声，而是实际的基频处能量。提取的第m个数据块对应的哨声信号基频频率采样点值为index(m)/2，避免在某些频率点，哨声信号二次谐波的能量大于基波的能量所导致的提取哨声信号基频频率采样点错误的情况。

在仿真实验中，对于同一海豚哨声信号时频谱轮廓进行提取，令步骤2的判决式中的参数β的值为无穷大，采用最大值提取法得到的短时傅里叶变化的基频频率值曲线提取结果图如图3所示；令步骤2的判决式中的参数β的值为2，得到的短时傅里叶变化的基频频率值曲线提取结果图如图4所示。

步骤3：对步骤2得到的短时傅里叶变化的基频频率值进行滤波，进一步去除其中的野点，得到平滑的基频轮廓。

滤波可以采用动窗加权平均法。

动窗加权平均法是一种信号平滑算法，要求随机误差服从高斯分布。设残差是输出信号s_m和信号中值的差，用r_m来表示。则

r_m＝s_m-median(s_m)          (3)

其中median(s_m)表示信号的中值，且有：

由于随机误差服从高斯分布，进而残差r_m也是服从高斯分布的，高斯分布的前提下，残差的标准差σ_m大约是mean(|r_m|)的1.3倍，并引入常数项10^-12防止标准差为0，则σ_m为：

σ_m＝1.3*mean(|r_m|)+10^-12          (5)

利用残差的标准差σ_m定义一个门限值th_m，且：

th_m＝4.6*σ_m          (6)

则对于动窗加权平均的权系数w_m，当|r_m(i)|≤th_m时有：

w_m＝(1-(r_m(i)/th_m)²)²          (7)

当|r_m(i)|>th_m时有：

w_m＝0          (8)

由式(7)可以看出，权系数w_m随着|r_m(i)|减小而趋向于1，反之趋向于0。由于是残差的标准偏差，当|r_m(i)|>th_m时，即r_m的绝对值大于或者等于4.6倍的σ_m时，说明r_m落在了-4.6σ_m到+4.6σ_m范围外。当不存在野值点时，数据落在-4.6σ_m到+4.6σ_m范围外的概率就很小了，只有0.000002。

根据动窗加权平均的权系数w_m，可以得到动窗加权平均的输出：

$z_m\left(j\right)=\frac{\underset{i=m-n+1}{\overset{m}{\mathrm{\Σ}}}{s}_m\left(i\right)\·w_m\left(i\right)}{\underset{i=m-n+1}{\overset{m}{\mathrm{\Σ}}}{w}_m\left(i\right)}---\left(9\right)$

根据上述方法可以去除野点，得到平滑的基频轮廓。

在仿真中采用动窗加权平均滤波法，对经步骤2得到的提取结果进行滤波，得到的滤波后的基频轮廓如图5所示。

滤波可以采用圆筛滤波法。

圆筛滤波法要求有效样本大于50％。若f(n)为需处理的样本，其中1<n<N，N为待处理样本总数。对于第k点样本数据f(n_k)，其中N₀≤k≤N，判断其是否为野点时，需用到f(n_i)得到判据，k-N₀<i<k-1。不妨令Q和Q₀为判断f(n_k)是否为野点的判据，则

$Q=(\underset{i=1}{\overset{N_0}{\mathrm{\&Sigma;}}}{(n_k-n_{k-i})}^2+{(f\left(n_k\right)-f\left(n_{k-i}\right))}^2),N_0\&le;k\&le;N---\left(10\right)$

$Q_0=\left(\underset{i=1}{\overset{N_0}{\mathrm{\&Sigma;}}}{(n_k-n_{k-i})}^2\right)+N_0*b^2,N_0\&le;k\&le;N---\left(11\right)$

上式中，b为筛选样本时的跳动容限。令

$\mathrm{\&mu;}=\frac{1}{Q},{\mathrm{\&mu;}}_0=\frac{1}{Q_0}---\left(12\right)$

μ₀为筛选样本的门限，若μ≥μ₀则该点为有效样本，否则判定为野点。通常在实际应用时，筛选样本的门限μ₀根据有效样本数据在总样本中的比重而加以适当权重，即μ₀*ε，0<ε≤1。有效样本数据越多，ε应趋向1，反之趋向0。

设l(n_k)为圆筛剔野点后的结果，对于第k点样本数据f(n_k)，经过圆筛剔野点后，有：

$\left\{\begin{array}{cc}l\left(n_k\right)=f\left(n_k\right),& \mathrm{\&mu;}\&GreaterEqual;{\mathrm{\&mu;}}_0*\mathrm{\&epsiv;}\\ l\left(n_k\right)=f\left(n_{k-1}\right)& \mathrm{\&mu;}\&le;{\mathrm{\&mu;}}_0*\mathrm{\&epsiv;}\end{array}\right.---\left(13\right)$

即f(n_k)的判据μ不小于设定的筛选样本门限时，圆筛剔野点后的结果为其本身，若f(n_k)的判据μ小于设定的筛选样本门限时，圆筛剔野点后的结果为该样本点的前一点f(n_k-1)。

根据上述方法可以去除野点，得到平滑的基频轮廓。

在仿真中采用圆筛滤波法，对经步骤2得到的提取结果进行滤波，得到的滤波后的基频轮廓如图6所示。

滤波可以采用卡尔曼滤波法。

引入一个离散时间过程，该过程可用线性随机差分方程描述，如下式所示：

X_k＝AX_k-1+BU_k-1+W_k-1          (14)

其中A为状态转移矩阵，B为系统输入量和系统状态之间的关系矩阵，U为系统的输入量，W_k-1表示过程激励噪声，其相应的量测方程如下式所示：

Z_k＝HX_k+V_k          (15)

其中H为观测量与系统状态之间的关系矩阵，V_k表示观测噪声。

卡尔曼滤波以反馈的方式估计系统的状态，它分为两个部分：时间更新过程和测量更新过程。

时间更新。时间更新过程又称为预测过程。该过程在上个历元状态估计值的基础上预测当前历元的状态值。离散卡尔曼滤波器的时间更新方程如式(16)和(17)所示。

${\hat{X}}_k^{-}=A{\hat{X}}_{k-1}+B{U}_{k-1}---\left(16\right)$

$P_k^{-}={\mathrm{AP}}_{k-1}{A}^T+Q---\left(17\right)$

其中^-代表先验，^∧代表估计，^T代表转置，Q为过程噪声的协方差。式(17)是用来衡量式(16)中状态估计值可靠性的误差协方差估计。时间更新将系统状态估计值和协方差估计从k-1时刻推算到k时刻，并将其作为k+1时刻为的先验估计。

测量更新。测量更新过程又称为校正过程。该过程利用实际测量值校正预测得到的状态先验估计值。离散卡尔曼滤波器的测量更新方程如式(18)、(19)和(20)所示。

${\hat{X}}_k={\hat{X}}_k^{-}+K_k[Z_k-H{\hat{X}}_k^{-}]---\left(18\right)$

$K_k=P_k^{-}{H}^T{({\mathrm{HP}}_k^{-}{H}^T+R)}^{-1}---\left(19\right)$

$P_k=(I-K_{k}H)P_k^{-}---\left(20\right)$

式(18)中的K_k为卡尔曼增益，被称为残余，它反映了预测值和实际值之间的不一致程度，该式表明了后验估计由先验估计值以及测量值结合构造而成。

式(19)是卡尔曼增益K的一种表示形式，其中R为量测方程中的噪声协方差。卡尔曼增益K作用是使达到最小，而将式(18)带入中，再将e_k带入P_k的定义式中，得到期望后，将P_k对K求导。令一阶导数为0，从而求得K。式(20)是用来衡量式(18)中后验估计值可靠性的后验误差协方差。

完成k时刻的时间更新和测量更新后，将k时刻得到的后验估计作为k+1时刻先验估计，整个过程再次重复。

在卡尔曼滤波器的实际实现过程中，测量噪声协方差R一般通过离线获取的系统观测值去进行计算得到。由于无法直接观测到过程信号X_k，所以过程噪声的协方差Q比较难确定，这种情况下，通常的做法是选择一个适当的Q值，并以此向X_k“注入”不确定性，建立一个简单的状态方程，当然该状态方程并不精准，但其结果可以接受，此时通常要求观测值可信。以上的两种情况下，不管选择系数的标准是否恰当合理，从统计学上看来都能通过调节滤波器系数达到需要的滤波效果。

根据上述方法可以去除野点，得到平滑的基频轮廓。

在仿真中采用卡尔曼滤波法，对经步骤2得到的提取结果进行滤波，得到的滤波后的基频轮廓如图7所示。

滤波可以采用动窗加权平均法、圆筛滤波法、卡尔曼滤波法三种滤波法组合串行级联的方式。

根据上述方法可以去除野点，得到平滑的基频轮廓。

在仿真中采用圆筛滤波法和卡尔曼滤波法串行级联进行滤波，对经步骤2得到的提取结果进行滤波，得到的滤波后的基频轮廓如图8所示。

步骤4：对滤波后的基频曲线进行最小二乘多项式拟合，得到一个连续变化的光滑哨声信号时频谱轮廓线。

由于哨声信号瞬时频率在短时间内被认为是连续变化的，所以对原始数据进行短时傅里叶变换时，由于噪声的影响和数据处理长度的限制，瞬时频率并不是很好的连续变换，所以对滤波后的曲线进行最小二乘多项式拟合，力求得到一个连续变化的哨声信号时频谱轮廓线。

设(x_i,y_i),i＝0,1,2,...,m为一组实验数据，引入函数p(x)，令：

$\left(x\right)=\underset{k=0}{\overset{n}{\mathrm{\&Sigma;}}}{a}_k{q}_k\left(x\right),n\&le;m---\left(21\right)$

最小二乘曲线拟合旨在确定p(x)使得的值最小，从几何意义上看，即确定与给定数据点的距离平方和为最小的曲线p(x)。p(x)的形式需要根据经验进行判断。本文中采用的是最小二乘多项式拟合方法，则：

$p_n\left(x\right)=\underset{k=0}{\overset{n}{\mathrm{\&Sigma;}}}{a}_k{x}^k\&Element;\mathrm{\&Phi;}---\left(22\right)$

其中，Φ为所有次数不大于n，其中n≤m的多项式构成的类。此时需要确定p(x)满足L最小：

$L=\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}[p_n\left(x_i\right)-y_i{]}^2=\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{(\underset{k=0}{\overset{n}{\mathrm{\&Sigma;}}}{a}_k{x}_i^k-y_i)}^2---\left(23\right)$

式(23)表明L为a₀,a₁,...a_n的函数，因此只要求得L的极值便可确定p(x)转化。对L求导，有：

$\frac{\&PartialD;L}{\&PartialD;a_j}=2\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}(\underset{k=0}{\overset{n}{\mathrm{\&Sigma;}}}{a}_k{x}_i^k-y_i)x_i^j=0,j=\mathrm{0,1},...,n---\left(24\right)$

即：

$\underset{k=0}{\overset{n}{\mathrm{\&Sigma;}}}\left(\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i^{k+j}\right)a_k=\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i^j{y}_i,j=\mathrm{0,1},...,n---\left(25\right)$

式(25)的矩阵形式为：

$\left[\begin{array}{cccc}m+1& \underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i& ...& \underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i^n\\ \underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i& \underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i^2& ...& \underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i^{n+1}\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ \underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i^n& \underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i^{n+1}& ...& \underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i^{2n}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ .\\ .\\ .\\ a_n\end{array}\right]=\left[\begin{array}{c}\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{y}_i\\ \underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i{y}_i\\ .\\ .\\ .\\ \underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{x}_i^n{y}_i\end{array}\right]---\left(26\right)$

式(26)的系数矩阵为对称正定矩阵，存在唯一解，将数据带入可解出a₀,a₁,...a_n，从而得到多项式p(x)。

由此可以得到一个连续变化的光滑哨声信号时频谱轮廓线。实现海豚哨声信号基波时频谱轮廓提取。

仿真中，图9给出了最小二乘多项式拟合处理后的结果和残差。

由得到的连续变化哨声信号基波时频谱轮廓曲线，可以依据谐波与基波的倍数关系合成任意次谐波的时频轮廓曲线。

仿真中，得到的基波时频谱轮廓与二次谐波时频谱轮廓如图10所示。

Claims (5)

1.一种海豚哨声信号时频谱轮廓提取方法，其特征在于：

(1)对真实的海豚哨声信号作短时傅里叶变换，获得海豚哨声时频信号X_m(k)；

(2)通过 $\left\{\begin{array}{c}[e\left(m\right),\mathrm{index}\left(m\right)]=\max \left|X_m\left(k\right)\right|\\ e{\left(m\right)}^{\&prime;}=X_m\left(\frac{\mathrm{index}\left(m\right)}{2}\right)\end{array}\right.$ 得到第m个数据块短时傅里叶变换结果的最大值e(m)、傅里叶变换后能量最大值e(m)对应的频率的采样点index(m)以及在index(m)/2处对应的能量值e(m)′，通过判决式：e(m)′>βP₀确定第m个数据块对应的哨声信号基频频率，当判决成立时，提取的第m个数据块哨声信号基频频率采样点为index(m)/2，当判决不成立时，提取的第m个数据块哨声信号基频频率采样点为index(m)，其中，P₀表示噪声的能量，系数β的值大于等于2；

(3)对步骤(2)得到的短时傅里叶变化的基频频率值进行滤波，去除其中的野点，得到平滑的基频轮廓；

(4)对步骤(3)滤波后的基频曲线进行最小二乘多项式拟合，得到一个连续变化的光滑哨声信号时频谱轮廓线。

2.根据权利要求1所述的一种海豚哨声信号时频谱轮廓提取方法，其特征在于：所述滤波方法为动窗加权平均法。

3.根据权利要求1所述的一种海豚哨声信号时频谱轮廓提取方法，其特征在于：所述滤波方法为圆筛滤波法。

4.根据权利要求1所述的一种海豚哨声信号时频谱轮廓提取方法，其特征在于：所述滤波方法为卡尔曼滤波法。

5.根据权利要求1所述的一种海豚哨声信号时频谱轮廓提取方法，其特征在于：所述滤波方法为动窗加权平均法、圆筛滤波法、卡尔曼滤波法三种滤波法组合串行级联的方式。