CN103714810A - 基于Gammatone滤波器组的车型特征提取方法 - Google Patents

Abstract

本发明公开了基于Gammatone滤波器组的车型特征提取方法，属于模式识别领域，涉及车辆辐射声信号的特征提取方法，具体来讲是一种通过计算车辆声信号在Gammatone滤波器组下的倒谱系数模拟人耳听觉特性的特征提取方法。该方法利用Gammatone滤波器组可模拟人耳非线性频率分解的特性，将车辆声信号滤波划分为不同子带信号并求取倒谱系数。本发明基于听觉特性中频率分解原理，对车辆声信号提取Gammatone倒谱系数，得到原始信号频带-能量特征，其中涉及计算均为常用信号处理技术，原理简单，步骤明了，便于编程实现，适用性广。

Description

基于Gammatone滤波器组的车型特征提取方法

技术领域

本发明属于模式识别领域，涉及车辆辐射声信号的特征提取方法，具体来讲是一种计算车辆声信号在Gammatone滤波器组下的倒谱系数的特征提取方法，所得特征可模拟人耳的听觉特性。

背景技术

车型识别的过程主要包括特征提取与分类器的训练识别两个部分。特征提取是模式识别的关键技术之一，利用车辆运动产生声信号对其进行分类识别是车型识别的一种重要方式。车型识别被广泛应用于交通管理、区域防护与要地警戒等领域。

目前，车辆声信号的特征提取主要基于信号的频率特征，AR参数模型是广泛采用的计算信号功率谱的频谱分析方法，有严格的理论支持与成熟的技术实现，但对于不同环境条件下的应用，提高特征分类精度的难度较大；此外还有基于小波及小波包的能量系数特征提取方法，进一步提高了特征的分类精度，但其原理深奥、实现复杂；近几年还出现了采用高阶谱分析、Mel倒谱系数、经验模式分解等新原理的基于频率的特征提取方法，可以得到相对精确的声信号特征，但整体而言，计算量大，实现复杂，且有的提取方法对背景噪声依赖较强，应用具有局限性。

发明内容

本发明目的是提供一种较为简单易行的基于Gammatone滤波器组的车型特征提取方法，该方法基于人耳非线性频率分解的听觉特性，利用听觉外周模型广泛采用的Gammatone滤波器组，过滤、划分车辆声信号并计算滤波得到子带信号的倒谱系数。

本发明的技术解决方案是：

一种基于Gammatone滤波器组的车型特征提取方法，其特殊之处在于：包括以下步骤：

1】采集原始车辆声信号s(n)，n表示采样数据序列中数据点的编号，采样率f_s满足奈奎斯特采样定理，即f_s≥2f_max，f_max为信号的最高频率；对采样率为f_s的原始车辆声信号s(n)进行预滤波、归一化、加窗分帧的预处理，得到时域的短时信号x(n)；相应计算如下所示：

信号预滤波：y(n)=s(n)+0.9375·s(n-1)  (1)

信号归一化： $\stackrel{\‾}{y}\left(n\right)=\frac{y\left(n\right)}{\left|y{\left(n\right)}_{\max }\right|}---\left(2\right)$

加窗分帧： $x\left(n\right)=\stackrel{\‾}{y}\left(n\right)\·w\left(n\right)---\left(3\right)$

所述加窗分帧采用重叠分帧，w(n)为加窗函数，采用hamming窗，其函数表达式如下：

$w\left(n\right)=h\left(n\right)=\left\{\begin{array}{cc}0.54-0.46\cos [2\mathrm{\πn}/(N-1)],& 0\≤n\≤N-1\\ 0,& \mathrm{else}\end{array}\right.---\left(4\right)$

其中，N为窗函数数据序列长度，即序列包含数据点数；n为此数据序列中任一数据点编号；

2】确定滤波器组应用的频率范围[f_L,f_H](f_L为滤波范围的频率下界，f_H为滤波范围的频率上界)、滤波器个数M与滤波器的阶数r，实现M个一组的Gammatone滤波器组：

$g_m\left(n\right)=\frac{1}{2\mathrm{\πj}}{\∫G}_m\left(z\right)z^{n-1}\mathrm{dz}---\left(5\right)$

其中，G_m(z)为滤波器在离散系统中的对应函数，z为Z变换中的复变量，n为离散系统中数据序列的任一数据点编号(n为整数)，j为虚数单位，m表示滤波器组中任一滤波器的编号，

3】短时信号x(n)通过Gammatone滤波器组，划分为M个子带信号(M为滤波器个数)x_m(n)：

x_m(n)=x(n)*g_m(n)  1≤m≤M  (6)

4】对子带信号x_m(n)作N点的FFT计算，N为子带信号数据序列长度，得到子带信号的功率谱X_m(k)，进一步对其取模平方，得到子带信号的能量谱E_m(k)；再除以子带信号帧长度N得到信号的功率谱G_m(k)；对各子带信号平均功率谱G_m(k)加和并取对数，得到子带信号的对数功率值，相应计算如下所示：

子带信号能量谱：E_m(k)=|X_m(k)|²  (7)

子带信号功率谱： $G_m\left(k\right)=\frac{E_m\left(k\right)}{N}---\left(8\right)$

子带信号对数功率值： $e_m=\log \left(\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{G}_m\left(k\right)\right),1\≤m\≤M---\left(9\right)$

其中，k代表子带信号数据序列中任一点编号，m代表子带信号编号，与各滤波器编号对应，

5】得到车型识别特征系数，实现车型识别：

5.1】对各子带信号的对数功率值组成的数据序列e(m)按照定义进行离散余弦变换，得到子带信号的原始p阶倒谱系数C(n)，转换公式为：

$C\left(n\right)=\sqrt{\frac{2}{M}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}e\left(m\right)\cos \left[\frac{\mathrm{\πn}}{M}(m-0.5)\right],n=\mathrm{1,2},...,p---\left(10\right)$

式中：

e(m)={e₁,…,e_m,…,e_M}；

n代表得到倒谱系数数据序列点的编号；

p为倒谱系数的阶数；

5.2】对得到的原始p阶倒谱系数C(n)，根据半正弦窗函数表达式(11)进行升半正弦倒谱提升，得到提升后的车型识别特征系数，如式（12）所示：

w(i)=1+6×sin(πi/N),1≤i≤N  (11)

C_G(n)=C(n)×w(i)  (12)

其中，N对应倒谱系数的阶数p，N=p，i代表1到N的正整数，式(12)中’×’表示C(n)与w(i)两数据序列作点乘运算，即对应位置数据点相乘。

上述步骤2实现M个一组的Gammatone滤波器组，具体如下：

2.1】根据式(13)、式(14)计算得到各滤波器的中心频率及带宽：

中心频率f_m为：

$f_m=(f_H+228.7)\exp (-\frac{v_i}{9.26})-228.7---\left(13\right)$

其中，f_H为滤波器的截止频率上界，v_i是滤波器重叠因子，用于指定相邻滤波器之间的重叠百分比，

再由中心频率f_m计算带宽，表达式为

$b_m=\mathrm{ERB}\left(f_m\right)=24.7\×(4.37\frac{f_m}{1000}+1)---\left(14\right)$

2.2】对各参数对应的Gammatone滤波器，其冲击响应的典型模式为：

其中，A为滤波器增益，r为滤波器阶数，b_m为滤波器带宽，f_m为滤波器的中心频率，

为相位，U(t)为阶跃函数，t为时域变量符号，M为滤波器个数，m表示滤波器组中任一滤波器的编号，简化模型中，取

对上式中g_m(t)按照拉普拉斯变换定义，计算得到滤波器在复频域的对应函数G_m(s)：

$G_m\left(s\right)={\∫}_{-\∞}^{\∞}{g}_m\left(t\right)e^{-\mathrm{st}}\mathrm{dt}$

$=\frac{A}{2}{\∫}_0^{\∞}{t}^{r-1}{e}^{(-2\mathrm{\πERB}\left(f_m\right)t)}(e^{j2\mathrm{\π}{f}_{m}t}+e^{-j2\mathrm{\π}{f}_{m}t})e^{-\mathrm{st}}\mathrm{dt}$

$=\frac{A}{2}[\frac{(r-1)!}{{(s+b-\mathrm{j\ω})}^n}+\frac{(r-1)!}{{(s+b+\mathrm{j\ω})}^n}]---\left(16\right)$

其中，A为滤波器增益，r为滤波器阶数，f_m是中心频率，

是相位，U(t)为阶跃函数，s为复频率，j为虚数单位，b=2πERB(f_m),ω=2πf_m,m表示滤波器组中任一滤波器的编号；

2.3】根据拉普拉斯变换中s平面与Z变换中z平面的映射关系，将G_m(s)转换为离散系统中Z变换的G_m(z)，得到滤波器在离散系统中z变换域的表示形式，再由逆z变换得定义，计算得到离散系统中一组M个的Gammatone滤波器的单位冲击响应：

$g_m\left(n\right)=\frac{1}{2\mathrm{\πj}}\∫G_m\left(z\right)z^{n-1}\mathrm{dz}---\left(17\right)$

其中，G_m(z)为滤波器在离散系统中的对应函数，z为Z变换中的复变量，n为离散系统中数据序列的任一数据点编号，n为整数，j为虚数单位，m表示滤波器组中任一滤波器的编号。

本发明所述的车辆声信号特征方法，基于听觉特征原理，通过计算滤波后子带信号倒谱系数，得到原始信号频带-能量特征，其中所涉及计算均为常用的信号处理技术，原理简单，步骤明了，便于编程实现，适用性广。根据该方法得到的特征具有低维数、高可分性的优点

附图说明

图1是基于Gammatone滤波器组的车辆声信号特征提取方法计算流程图。

图2是一组17个的Gammatone滤波器幅频响应曲线图。

图3是两个不同类型的车辆声信号Gammatone特征向量示意图。

具体实施方式

本发明应用听觉特性中频率分解原理，对车辆声信号进行Gammatone倒谱系数提取的流程为：

参阅图1，基于Gammatone滤波器组提取车辆声信号特征的方法为：

1）对采样率为f_s（满足奈奎斯特采样定理，即f_s≥2f_max，f_max为信号的最高频率）的原始车辆声信号s(n)进行预滤波、归一化、加窗分帧的预处理，得到时域的短时信号帧x(n)。相应计算如下所示：

信号低频加强：y(n)=s(n)+0.9375·s(n-1)  (1)

信号归一化： $\stackrel{\‾}{y}\left(n\right)=\frac{y\left(n\right)}{\left|y{\left(n\right)}_{\max }\right|}---\left(2\right)$

加窗分帧： $x\left(n\right)=\stackrel{\‾}{y}\left(n\right)\·w\left(n\right)---\left(3\right)$

车辆声信号是非平稳随机信号，在基于信号短时平稳的基础上，再应用平稳信号处理方法分析。车辆信号在0.2s～1s时间段内可近似平稳，通过对信号加窗，将原始信号分割为帧片段。此处，采用“重叠分帧”的方法，即前一帧和后一帧包含重叠数据，更好的保证帧间的连续性。

上式中，w(n)为加窗函数，常采用hamming窗，它的函数表达式如下：

$w\left(n\right)=h\left(n\right)=\left\{\begin{array}{cc}0.54-0.46\cos [2\mathrm{\πn}/(N-1)],& 0\≤n\≤N-1\\ 0,& \mathrm{else}\end{array}\right.---\left(4\right)$

2）根据滤波器覆盖频率范围[f_L,f_H]、滤波器个数M及滤波器阶数r(r=4)，实现M个一组的Gammatone滤波器组（在信号采集过程中，信号采样率是确定的，以上其余参数由信号处理者根据实际需求设定），具体步骤如下。

首先，根据式(5)、式(6)计算得到各滤波器的f_m及带宽b_m。

$f_m=(f_H+228.7)\exp (-\frac{v_i}{9.26})-228.7---\left(5\right)$

$b_m=\mathrm{ERB}\left(f_m\right)=24.7\×(4.37\frac{f_m}{1000}+1)---\left(6\right)$

再根据Gammatone滤波器的典型冲击响应函数（式(7)）做拉普拉斯变换，得到滤波器复频域的表达式G_m(s)。

$G_m\left(s\right)={\∫}_{-\∞}^{\∞}{g}_m\left(t\right)e^{-\mathrm{st}}\mathrm{dt}$

$=\frac{A}{2}{\∫}_0^{\∞}{t}^{r-1}{e}^{(-2\mathrm{\πERB}\left(f_m\right)t)}(e^{j2\mathrm{\π}{f}_{m}t}+e^{-j2\mathrm{\π}{f}_{m}t})e^{-\mathrm{st}}\mathrm{dt}$

$=\frac{A}{2}[\frac{(r-1)!}{{(s+b-\mathrm{j\ω})}^n}+\frac{(r-1)!}{{(s+b+\mathrm{j\ω})}^n}]---\left(8\right)$

其中，A为滤波器增益，r为滤波器阶数，f_m是中心频率，

是相位，U(t)为阶跃函数，s为复频率，j为虚数单位，b=2πERB(f_m),ω=2πf_m,m表示滤波器组中任一滤波器的编号。

最后，将G_m(s)转换为Z变换的G_m(z)形式，再按照逆Z变换定义计算得到离散系统下一组M个的Gammatone滤波器的单位冲击响应：

$g_m\left(n\right)=\frac{1}{2\mathrm{\πj}}\∫G_m\left(z\right)z^{n-1}\mathrm{dz}---\left(9\right)$

图2所示是10～2500Hz范围内的17通道4阶Gammatone滤波器组的对数幅频响应曲线。在线性频率上，随着滤波器编号的增大，相邻滤波器中心频率的带宽逐渐增大。

3）短时信号帧x(n)通过Gammatone滤波器组，划分为M个子带信号(M为滤波器个数)x_m(n)：

x_m(n)=x(n)*g_m(n)  1≤m≤M  (10)

4）求子带功率谱。对子带信号x_m(n)进行快速傅里叶变换(Fast FourierTransform，FFT)，得到子带信号的功率谱X_m(k)，进一步取模平方，除以信号长度N得到信号的功率谱G_m(k)。取对数，得到子带信号的对数功率谱。相应计算如下所示：

子带信号能量谱：E_m(k)=|X_m(k)|²  (11)

子带信号平均功率谱： $G_m\left(k\right)=\frac{E_m\left(k\right)}{N}---\left(12\right)$

子带信号对数功率值： $e_m=\log \left(\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{G}_m\left(k\right)\right),1\≤m\≤M---\left(13\right)$

5）求子带倒谱系数。将各子带信号的对数能量谱进行离散余弦变换(Discrete Cosine Transform，DCT)。转换公式为：

$C\left(n\right)=\sqrt{\frac{2}{M}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}e\left(m\right)\cos \left[\frac{\mathrm{\πn}}{M}(m-0.5)\right],n=\mathrm{1,2},...,p---\left(14\right)$

式中，e(m)为各子带信号对数功率值组成序列，即e(m)={e₁,…,e_m,…,e_M}；p为倒谱特征系数的阶数。C(n)为得到的原始p阶倒谱系数。再对DCT得到的特征进行升半正弦倒谱提升，式（15）为半正弦窗函数，得到提升后的特征如式（16）所示：

w(i)=1+6×sin(πi/N),1≤i≤N  (15)

C_G(n)=C(n)×w(i)  (16)

选取两类车辆目标水陆两用作战车（Assault AmphibianVehicle，AAV）、牵引式货车（Dragon Wagon，DW）运动辐射声信号作为样本库，信号采样率f_s=4960Hz。基于Gammatone滤波器组提取两类车辆声信号特征的方法为：参阅图1：

1）声信号经预滤波、分帧、加窗等预处理后，得到时域的短时信号帧x(n)。相应计算如下所示：

信号预滤波：y(n)=s(t/f_s)+0.9375·s(t/f_s-1)  (17)

信号归一化： $\stackrel{\‾}{y}\left(n\right)=\frac{y\left(n\right)}{\left|\stackrel{\‾}{y}{\left(n\right)}_{\max }\right|}---\left(18\right)$

加窗分帧： $x\left(n\right)=\stackrel{\‾}{y}\left(n\right)\·w\left(n\right)---\left(3\right)$

本例中，Matlab计算如下：

预滤波：y(n)=filter([10.9375],1,s(n));

信号归一化：y(n)=y(n)/max(abs(y(n)));

分帧：xx=enframe(y(n),N1,N2);

其中s(n)为采集待处理的原始车辆声信号，y(n)为信号滤波处理中间结果，N1为信号帧长度，取α=0.2065，N1=FrameSize=α·f_s=1024点；N2为帧重叠点数，取β=0.5，N2=β·N1=512点。

2）根据滤波器组应用的频率范围[f_L,f_H](f_L为滤波范围的频率下界，f_H为滤波范围的频率上界)、滤波器个数M与滤波器的阶数r(r=4)，计算得到滤波器的中心频率f_m及等效带宽b_m，计算M个一组的Gammatone滤波器组滤波参数。本例中，取M=17，[f_L,f_H]=[10,2500]Hz，其他输入采用函数默认值，得到一组17个的Gammatone滤波器组的滤波系数、中心频率、等效带宽及群延时，调用Matlab中voicebox工具包中gammabank函数，计算如下：

[b,a,fx,bx,gd]=gammabank(17,4960,′′,[102500])；

gammabank函数原型为：

[b,a,fx,bx,gd]=gammabank(n,fs,w,fc,bw,ph,k)其中，输入参数n为滤波器个数，fs为信号采样率（单位：Hz），w为可选频率转化方式，fc为滤波器中心频率，bw为滤波器带宽，ph为滤波器冲击函数响应的相位，k为滤波器阶数。此处，需确定滤波器个数n，采样频率fs及中心频率覆盖范围[f_L,f_H]，其他取函数默认参数值。函数输出包括系统模型中滤波器系数a,b,滤波器中心频率点fx，带宽值bx与各滤波器群延时gd。

gammabank函数采用传递函数模型描述Gammatone滤波器组，传递函数H(z)的一般表达式如下所示：

$H\left(z\right)=\frac{Y\left(z\right)}{X\left(z\right)}=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}+\·\·\·+b_{N-1}{z}^{-(N-1)}}{1+a_1{z}^{-1}+a_2{z}^{-2}+\·\·\·+a_{M-1}{z}^{-(M-1)}}$

式中，X(z)表示输入x(n)的Z变换，Y(z)表示输出y(n)的Z变换，b_i(i=1,…,m)和a_i(i=1,…,n)是滤波器的分子分母系数。由Gammatone滤波器的时域冲击响应长度可知，Gammatone滤波器是无限长冲击响应（IIR）滤波器，其阶数为H(z)的极点个数。gammabank计算滤波器组系数的主要步骤如下：

(1)将线性Hz频率表示的频带范围转化为ERB频率刻度表示范围：

g=abs(frq);

erb=11.17268*sign(frq).*log(1+46.06538*g./(g+14678.49));

其中frq为待转换线性Hz频率，分别对滤波器组覆盖最小、最大频率变换得到ERB刻度下的频率分布范围[fx1,fx2];

(2)在ERB刻度下划分排列各滤波器中心频率：

fx=linspace(fx1,fx2,n);

上式表示取[fx1,fx2]范围的n等分点，其中n为滤波器个数，得到fx为n个滤波器的中心频率。

(3)计算各滤波器带宽：

bnd=6.23e-6*g.^2+93.39e-3*g+28.52;

bx=1.019bnd;

其中g为各中心频率在ERB刻度下表示频率，bnd为ERB刻度下各中心频率对应带宽，bx为转化到线性Hz频率下对应带宽。

(4)计算滤波器系数a，b：

ww=exp((1i*fx-bx)*2*pi/fs);%计算滤波器起始频率对应复变量

a=round([1cumprod((-k:-1)./(1:k))]);%构造k(k=4)阶的二项式系数

b=conv(a,(0:k-1).^(k-1));%通过卷积计算得到k阶二项式对应初始分子系数

b=exp(1i*ph)*b(1:k);%计算复数空间下滤波器相位对应各阶系数值

wwp=repmat(ww,1,k+1).^repmat(0:k,nf,1);%计算各通道中

滤波器起始频率对应各阶系数

denc=repmat(a,nf,1).*wwp;%计算各通道中滤波器系数a对应分式方程系数

numc=b.*wwp(:,1:k);%计算各通道中滤波器系数b对应分式方程系数

ww=exp(2i*fx*pi/fs);%计算各通道滤波器中心频率对应复变量

u=polyval(b(i,:),ww(i));%构造系数b下复数空间二项式

v=polyval(a(i,:),ww(i));%构造系数a下复数空间二项式

b_i=real(conv(numc_i,conj(denc_i)))*abs(v-u)%计算各通道滤波器分子系数bi

a_i=real(conv(denc_i,conj(denc_i)))%计算各通道滤波器分母系数ai

式中，i为虚数单位，ph为滤波器相位取为0，k为滤波器阶数取为4，ww表示传递函数中复变量z，fx为滤波器中心频率，fs为信号采样率。

3）短时信号经Gammatone滤波器组，划分为M个子信号(M为滤波器个数)。本例中，取两类车辆信号各一帧x₁(n)，x₂(n)，在Matlab平台下，调用filter函数计算信号帧经过Gammatone滤波器组输出。

y=filterbank(b,a,xx(i,:),gd);

4）求子带功率谱。本例中，对两类车辆信号帧经Gammatone滤波后的子带信号分别进行快速傅里叶变换，取模平方得到子信号的能量谱G₁(k)，G₂(k)。将子信号能量谱在时间上平均后取对数，得到子带信号的对数功率谱。Matlab中计算如下：

r=abs(fft(y,1024));

e=log(sum(r.^2/1024));

5）求子带倒谱系数。将各子带信号的对数能量谱进行离散余弦变换(Discrete Cosine Transform，DCT)。转换公式为：

$C\left(n\right)=\sqrt{\frac{2}{M}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}e\left(m\right)\cos \left[\frac{\mathrm{\πn}}{M}(m-0.5)\right],n=\mathrm{1,2},...,p---\left(20\right)$

式中，e(m)为各子带信号对数功率值组成序列，即e(m)={e₁,…,e_m,…,e_M}，p为倒谱特征系数的阶数。C(n)为得到的原始p阶倒谱系数。再对DCT得到的特征进行升半正弦倒谱提升，式（21）为半正弦窗函数，得到提升后的特征如式（22）所示：

w(i)=1+6×sin(πi/N),1≤i≤N  (21)

C_G(n)=C(n)×w(i)  (22)

本例中，p=N=17，由于第一维数据主要代表信号能量值将其去除，余下的16维组成两类车辆信号的特征向量C_G1(n)，C_G2(n)。Matlab中计算如下：

c=dct(log(sum(r.^2、1024)));

c1=c.*w;

图3中为AAV、DW两类车辆信号各一帧提取得到特征，其中，横轴为特征维数，纵轴为特征幅度。可见，两类特征各维取值间距大，明显可分，进一近实验结果证明，相比传统Mel倒谱系数特征提取方法，在特征维数相同时，本发明方法可提高识别正确率5%～10%。

两类车辆信号Gammatone特征提取完整Matlab程序：

本发明原理：

本发明提取车辆声信号特征的方法，模拟人耳对声信号的非线性频率分解特性，提取的特征参数C_G(n)各维系数代表了车辆声信号按其组成频率的能量分布情况，有效地表征了车型特征，实现基于声信号的车型识别。

本发明基于人耳非线性频率分解的听觉特性，利用听觉外周模型广泛采用的Gammatone滤波器组，过滤、划分车辆声信号并计算滤波得到子带信号的倒谱系数，提供一种较为简单易行的车型特征提取方法，根据该方法得到的特征具有低维数、高可分性的优点。

Gammatone滤波器是一个标准的耳蜗听觉滤波器，该滤波器组冲击响应的典型模式为：

式中，A为滤波器增益，r为滤波器阶数，f_m是中心频率，是相位，b_m为等效带宽，M为滤波器个数，U(t)为阶跃函数。简化模型中，取A=1，r=4，

b_m=ERB(f_m)。

ERB(f_m)为等效矩形带宽(Equivalent Rectangular Bandwidth,ERB)，它决定了脉冲响应的衰减速度，与滤波器带宽有关，而每个滤波器带宽与人耳听觉临界频带(Critical Band,CB)有关，听觉心理学中，ERB(f_m)可由式(24)计算得到：

$\mathrm{ERB}\left(f_m\right)=24.7\×(4.37\frac{f_m}{1000}+1)---\left(24\right)$

其中，中心频率f_m为：

$f_m=(f_H+228.7)\exp (-\frac{v_i}{9.26})-228.7---\left(25\right)$

式中，f_H为滤波器的截止频率，v_i是滤波器重叠因子，用来指定相邻滤波器之间重叠的百分比。每个滤波器的中心频率确定后，相应的带宽可由式(24)计算得出。

Gammatone滤波器组是应用最多的耳蜗基底膜模型，在语音信号的处理中，主要用于语音信号增强及特征提取。由各种机动车辆目标产生的声信号是一种多声源噪声，在应用Gammatone滤波器组提取特征时，主要存在以下难点：(1)车辆运动声源包括气动噪声和机械噪声，信号成分较语音更复杂；(2)车辆声信号对环境变化十分敏感，路况、车速、自然环境都会引起信号变化；(3)需根据目标车辆声信号的优势频带范围，动态调整Gammatone滤波器覆盖频带。目前，对Gammatone滤波器处理其他目标声信号的研究相对较少，曾有文献提出利用听觉模型（Gammatone滤波器及高/低通滤波器组成混合滤波器组）输出谱特征用于声目标识别，并对水下目标进行试验测试，得到优于传统方法的识别结果。

Claims (2)

1.一种基于Gammatone滤波器组的车型特征提取方法，其特征在于：

包括以下步骤：

1】采集原始车辆声信号s(n)，n表示采样数据序列中数据点的编号，采样率f_s满足奈奎斯特采样定理，即f_s≥2f_max，f_max为信号的最高频率；对采样率为f_s的原始车辆声信号s(n)进行预滤波、归一化、加窗分帧的预处理，得到时域的短时信号x(n)；相应计算如下所示：

信号预滤波：y(n)=s(n)+0.9375·s(n-1)  (1)

信号归一化： $\stackrel{\‾}{y}\left(n\right)=\frac{y\left(n\right)}{\left|y{\left(n\right)}_{\max }\right|}---\left(2\right)$

加窗分帧： $x\left(n\right)=\stackrel{\‾}{y}\left(n\right)\·w\left(n\right)---\left(3\right)$

所述加窗分帧采用重叠分帧，w(n)为加窗函数，采用hamming窗，其函数表达式如下：

$w\left(n\right)=h\left(n\right)=\left\{\begin{array}{cc}0.54-0.46\cos [2\mathrm{\πn}/(N-1)],& 0\≤n\≤N-1\\ 0,& \mathrm{else}\end{array}\right.---\left(4\right)$

其中，N为窗函数数据序列长度，即序列包含数据点数；n为此数据序列中任一数据点编号；

2】确定滤波器组应用的频率范围[f_L,f_H](f_L为滤波范围的频率下界，f_H为滤波范围的频率上界)、滤波器个数M与滤波器的阶数r，实现M个一组的Gammatone滤波器组：

$g_m\left(n\right)=\frac{1}{2\mathrm{\πj}}{\∫G}_m\left(z\right)z^{n-1}\mathrm{dz}---\left(5\right)$

其中，G_m(z)为滤波器在离散系统中的对应函数，z为Z变换中的复变量，n为离散系统中数据序列的任一数据点编号(n为整数)，j为虚数单位，m表示滤波器组中任一滤波器的编号，

3】短时信号x(n)通过Gammatone滤波器组，划分为M个子带信号(M为滤波器个数)x_m(n)：

x_m(n)=x(n)*g_m(n)1≤m≤M  (6)

4】对子带信号x_m(n)作N点的FFT计算，N为子带信号数据序列长度，得到子带信号的功率谱X_m(k)，进一步对其取模平方，得到子带信号的能量谱E_m(k)；再除以子带信号帧长度N得到信号的功率谱G_m(k)；对各子带信号平均功率谱G_m(k)加和并取对数，得到子带信号的对数功率值，相应计算如下所示：

子带信号能量谱：E_m(k)=|X_m(k)|²  (7)

子带信号功率谱： $G_m\left(k\right)=\frac{E_m\left(k\right)}{N}---\left(8\right)$

子带信号对数功率值： $e_m=\log \left(\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{G}_m\left(k\right)\right),1\≤m\≤M---\left(9\right)$

其中，k代表子带信号数据序列中任一点编号，m代表子带信号编号，与各滤波器编号对应，

5】得到车型识别特征系数，实现车型识别：

5.1】对各子带信号的对数功率值组成的数据序列e(m)按照定义进行离散余弦变换，得到子带信号的原始p阶倒谱系数C(n)，转换公式为：

$C\left(n\right)=\sqrt{\frac{2}{M}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}e\left(m\right)\cos \left[\frac{\mathrm{\πn}}{M}(m-0.5)\right],n=\mathrm{1,2},...,p---\left(10\right)$

式中：

e(m)={e₁,…,e_m,…,e_M}；

n代表得到倒谱系数数据序列点的编号；

p为倒谱系数的阶数；

5.2】对得到的原始p阶倒谱系数C(n)，根据半正弦窗函数表达式(11)进行升半正弦倒谱提升，得到提升后的车型识别特征系数，如式（12）所示：

w(i)=1+6×sin(πi/N),1≤i≤N  (11)

C_G(n)=C(n)×w(i)  (12)

其中，N对应倒谱系数的阶数p，N=p，i代表1到N的正整数，式(12)中’×’表示C(n)与w(i)两数据序列作点乘运算，即对应位置数据点相乘。

2.根据权利要求1所述基于Gammatone滤波器组的车型特征提取方法，其特征在于：所述步骤2实现M个一组的Gammatone滤波器组，具体如下：

2.1】根据式(13)、式(14)计算得到各滤波器的中心频率及带宽：

中心频率f_m为：

$f_m=(f_H+228.7)\exp (-\frac{v_i}{9.26})-228.7---\left(13\right)$

其中，f_H为滤波器的截止频率上界，v_i是滤波器重叠因子，用于指定相邻滤波器之间的重叠百分比，

再由中心频率f_m计算带宽，表达式为

$b_m=\mathrm{ERB}\left(f_m\right)=24.7\×(4.37\frac{f_m}{1000}+1)---\left(14\right)$

2.2】对各参数对应的Gammatone滤波器，其冲击响应的典型模式为：

其中，A为滤波器增益，r为滤波器阶数，b_m为滤波器带宽，f_m为滤波器的中心频率，

为相位，U(t)为阶跃函数，t为时域变量符号，M为滤波器个数，m表示滤波器组中任一滤波器的编号，简化模型中，取

对上式中g_m(t)按照拉普拉斯变换定义，计算得到滤波器在复频域的对

应函数G_m(s)：

$G_m\left(s\right)={\∫}_{-\∞}^{\∞}{g}_m\left(t\right)e^{-\mathrm{st}}\mathrm{dt}$

$=\frac{A}{2}{\∫}_0^{\∞}{t}^{r-1}{e}^{(-2\mathrm{\πERB}\left(f_m\right)t)}(e^{j2\mathrm{\π}{f}_{m}t}+e^{-j2\mathrm{\π}{f}_{m}t})e^{-\mathrm{st}}\mathrm{dt}$

$=\frac{A}{2}[\frac{(r-1)!}{{(s+b-\mathrm{j\ω})}^n}+\frac{(r-1)!}{{(s+b+\mathrm{j\ω})}^n}]---\left(16\right)$

其中，A为滤波器增益，r为滤波器阶数，f_m是中心频率，

是相位，U(t)为阶跃函数，s为复频率，j为虚数单位，b=2πERB(f_m),ω=2πf_m,m表示滤波器组中任一滤波器的编号；

2.3】根据拉普拉斯变换中s平面与Z变换中z平面的映射关系，将G_m(s)转换为离散系统中Z变换的G_m(z)，得到滤波器在离散系统中z变换域的表示形式，再由逆z变换得定义，计算得到离散系统中一组M个的Gammatone滤波器的单位冲击响应：

$g_m\left(n\right)=\frac{1}{2\mathrm{\πj}}\∫G_m\left(z\right)z^{n-1}\mathrm{dz}---\left(17\right)$

其中，G_m(z)为滤波器在离散系统中的对应函数，z为Z变换中的复变量，n为离散系统中数据序列的任一数据点编号，n为整数，j为虚数单位，m表示滤波器组中任一滤波器的编号。

Images