CN113409819B - 一种基于听觉谱特征提取的直升机声信号识别方法 - Google Patents

一种基于听觉谱特征提取的直升机声信号识别方法 Download PDF

Info

Publication number
CN113409819B
CN113409819B CN202110951550.1A CN202110951550A CN113409819B CN 113409819 B CN113409819 B CN 113409819B CN 202110951550 A CN202110951550 A CN 202110951550A CN 113409819 B CN113409819 B CN 113409819B
Authority
CN
China
Prior art keywords
frequency
auditory
filter
helicopter
feature extraction
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.)
Active
Application number
CN202110951550.1A
Other languages
English (en)
Other versions
CN113409819A (zh
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.)
Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center
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 Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202110951550.1A priority Critical patent/CN113409819B/zh
Publication of CN113409819A publication Critical patent/CN113409819A/zh
Application granted granted Critical
Publication of CN113409819B publication Critical patent/CN113409819B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/48—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use
    • G10L25/51—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use for comparison or discrimination

Landscapes

  • Engineering & Computer Science (AREA)
  • Computational Linguistics (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)
  • Auxiliary Devices For Music (AREA)
  • Compression, Expansion, Code Conversion, And Decoders (AREA)

Abstract

本发明涉及直升机声信号识别，具体公开了一种基于听觉谱特征提取的直升机声信号识别方法，包括如下步骤：步骤1：分段加窗；步骤2：FFT分析；步骤3：尺度变换；步骤4：听觉滤波；步骤5：对数压缩；步骤6：求取均值；步骤7：分类识别。本发明的有益效果为：将非线性频率尺度变换和听觉滤波器引入到FFT分析与对数压缩之间，借助听觉计算模型的非线性频率选择能力以及更强的中低频分辨率和分析处理能力，使不易察觉的直升机声信号个性特征在若干分析频带内显露出来，提升直升机辨识的有效性和鲁棒性。

Description

一种基于听觉谱特征提取的直升机声信号识别方法
技术领域
本发明涉及直升机声信号识别，具体涉及一种基于听觉谱特征提取的直升机声信号识别方法。
背景技术
直升机因其独特的垂直起降、高机动性和低空突防能力，在灾难救治、局部打击和斩首行动中得到了越来多的应用。随着“黑飞”问题的逐渐突出和低空防御的日渐重要，如何对直升机进行探测、识别、定位和跟踪也引起了越来越多的关注和研究。直升机声信号，特别是其主桨和尾桨周期性扰动空气产生的强中低频噪声，是识别直升机的重要特征。
直升机声信号识别能在恶劣气象条件和高度遮挡情况下（如云雾、山岳和丛林等）有效弥补雷达、红外、光学等传统探测识别手段的不足，其关键环节是有效地提取出隐含在声信号中的能反应直升机目标类型的个性特征。典型的直升机声信号特征提取方法大致分为时域方法、频域方法、时频域方法和倒谱域方法等。时域特征提取直接对传声器采集得到的原始时域声信号进行统计分析，提取过零率、峰值位置、波形结构等多维特征，速度快，实时性好，但是低信噪比和复杂环境下，特征提取困难且识别性能迅速降低。其它的特征提取方法本质上均以频谱特性分析为基础，谱特征提取能力的提升对这些方法的性能改进具有重要意义。
发明内容
本发明的目的在于，针对上述问题，提出一种基于听觉谱特征提取的直升机声信号识别方法。
一种基于听觉谱特征提取的直升机声信号识别方法，包括如下步骤：
步骤1：分段加窗：将传声器采集的直升机原始声信号划分为长度为L的若干信号段，得到分段声信号
Figure 236683DEST_PATH_IMAGE001
，
Figure 259872DEST_PATH_IMAGE002
；将所述分段声信号
Figure DEST_PATH_IMAGE003
乘以窗函数
Figure 718535DEST_PATH_IMAGE004
，得到加窗声信号
Figure DEST_PATH_IMAGE005
；
步骤2：FFT分析：对步骤1中得到的加窗信号
Figure 404862DEST_PATH_IMAGE006
进行FFT分析，得到加窗声信号频谱，进而取加窗声信号频谱绝对值得到幅值谱；
步骤3：尺度变换：在频率分析范围
Figure DEST_PATH_IMAGE007
内，计算听觉滤波器组在非线性频率尺度变换下的中心频率f i ，
Figure 722711DEST_PATH_IMAGE008
；
Figure DEST_PATH_IMAGE009
为频率分析范围的下界频率，大于FFT分析的最低频率；
Figure 300323DEST_PATH_IMAGE010
为频率分析范围的上界频率，小于FFT分析的最高频率及奈奎斯特频率
Figure DEST_PATH_IMAGE011
的最小值，
Figure 868620DEST_PATH_IMAGE012
为声信号采样频率；所述听觉滤波器组包括Mel滤波器组或Gammatone滤波器组；所述非线性尺度变换包括Mel尺度变换或ERB尺度变换或Bark尺度变换；
步骤4：听觉滤波：先根据步骤3求解得到的中心频率
Figure DEST_PATH_IMAGE013
确定听觉滤波器的表达式，再用听觉滤波器组对步骤2中得到的幅值谱进行带通滤波处理，输出滤波处理后的幅值谱；
步骤5：对数压缩：对步骤4中的听觉滤波器组输出的幅值谱进行对数压缩；
步骤6：求取均值：求取步骤5中得到的对数压缩结果的均值得到最终的听觉谱特征；
步骤7：分类识别：按照上述步骤1-6分别求取训练集和测试集的听觉谱特征，先将训练集的听觉谱特征送入分类器进行训练，再将测试集的听觉谱特征送入分类器进行识别，确定直升机类型。
优选的，所述窗函数
Figure 975116DEST_PATH_IMAGE014
包括Hamming窗；窗函数
Figure DEST_PATH_IMAGE015
为：
Figure 780261DEST_PATH_IMAGE016
，
n为与分段声信号
Figure DEST_PATH_IMAGE017
对应的离散时间点，
Figure 646717DEST_PATH_IMAGE018
；L为分段声信号的长度。
优选的，所述Mel尺度变换与中心频率
Figure DEST_PATH_IMAGE019
的对应关系为：
Figure 548814DEST_PATH_IMAGE020
；
其中，
Figure DEST_PATH_IMAGE021
为频率分析范围的下界频率，
Figure 29474DEST_PATH_IMAGE022
为频率分析范围的上界频率，i表示滤波器的序号，
Figure DEST_PATH_IMAGE023
，N为滤波器个数；
Figure 899078DEST_PATH_IMAGE024
表示中心频率
Figure DEST_PATH_IMAGE025
对应的Mel尺度；
Figure 818493DEST_PATH_IMAGE026
为分析范围下界频率对应的Mel尺度，
Figure DEST_PATH_IMAGE027
为分析范围上界频率对应的Mel尺度；
所述ERB尺度变换与中心频率
Figure 325829DEST_PATH_IMAGE028
的对应关系为：
Figure DEST_PATH_IMAGE029
；
其中，
Figure 977390DEST_PATH_IMAGE030
为频率分析范围的下界频率，
Figure DEST_PATH_IMAGE031
为频率分析范围的上界频率，i表示滤波器的序号，
Figure 819444DEST_PATH_IMAGE032
，N为滤波器个数；
Figure DEST_PATH_IMAGE033
表示中心频率
Figure 791817DEST_PATH_IMAGE034
对应的ERB尺度；
Figure DEST_PATH_IMAGE035
为分析范围下界频率对应的ERB尺度，
Figure 402927DEST_PATH_IMAGE036
为分析范围上界频率对应的ERB尺度；
所述 Bark尺度变换与中心频率
Figure DEST_PATH_IMAGE037
的对应关系为：
Figure 225389DEST_PATH_IMAGE038
；
其中，
Figure DEST_PATH_IMAGE039
为频率分析范围的下界频率，
Figure 305472DEST_PATH_IMAGE040
为频率分析范围的上界频率，i表示滤波器的序号，
Figure DEST_PATH_IMAGE041
，N为滤波器个数；
Figure 566689DEST_PATH_IMAGE042
表示中心频率
Figure DEST_PATH_IMAGE043
对应的Bark尺度；
Figure 547152DEST_PATH_IMAGE044
为分析范围下界频率对应的Bark尺度，
Figure DEST_PATH_IMAGE045
为分析范围上界频率对应的Bark尺度。
由上述公式可知，进行逆变换即可反算出相应尺度下的中心频率
Figure 540516DEST_PATH_IMAGE046
。
优选的，所述Mel滤波器的传递函数为：
Figure DEST_PATH_IMAGE047
；
其中，i表示滤波器的序号，i=1，…，N，N为滤波器个数；
Figure 622742DEST_PATH_IMAGE048
表示Mel滤波器的传递函数；
Figure DEST_PATH_IMAGE049
为某一分析频率，
Figure 172803DEST_PATH_IMAGE050
，
Figure DEST_PATH_IMAGE051
。
优选的，所述Gammatone滤波器的时域表达式为：
Figure 758505DEST_PATH_IMAGE052
；
其中，
Figure DEST_PATH_IMAGE053
为Gammatone滤波器的时域表达式，A为滤波器增益；m为滤波器阶数；
Figure 922770DEST_PATH_IMAGE054
为滤波器带宽；
Figure DEST_PATH_IMAGE055
为中心频率，t为时间。
优选的，所述分类器包括基于欧式距离的最近邻分类器。
优选的，所述非线性频率尺度变换与听觉滤波器组组合得到听觉谱特征提取算法，包括M-M算法、E-M算法、B-M算法、M-G算、E-G算法和B-G算法；所述M-M算法是指Mel尺度变换和Mel滤波器组成的听觉谱特征提取算法；所述E-M算法是指ERB尺度变换和Mel滤波器组成的听觉谱特征提取算法；所述B-M算法是指Bark尺度变换和Mel滤波器组成的听觉谱特征提取算法；所述M-G算法是指Mel尺度变换和Gammatone滤波器组成的听觉谱特征提取算法；所述E-G算法是指ERB尺度变换和Gammatone滤波器组成的听觉谱特征提取算法；所述B-G算法是指Bark尺度变换和Gammatone滤波器组成的听觉谱特征提取算法。
优选的，所述分段声信号的长度L需使得FFT分析的频率分辨率
Figure 735700DEST_PATH_IMAGE056
与直升机的基频为同一数量级；
Figure DEST_PATH_IMAGE057
为声信号采样频率，L为分段声信号长度。
优选的，所述听觉滤波器组中的滤波器个数Ｎ需使得尺度变换的最小频率分辨率与直升机的基频为同一数量级。
与现有技术相比，采用上述技术方案的有益效果为：将非线性频率尺度变换和听觉滤波器引入到FFT分析与对数压缩之间，借助听觉计算模型的非线性频率选择能力以及更强的中低频分辨率和分析处理能力，使不易察觉的声信号个性特征在若干分析频带内显露出来，提升直升机辨识的有效性和鲁棒性。
附图说明
图1是基于传统谱特征提取的直升机声信号识别流程。
图2是本发明的基于听觉谱特征提取的直升机声信号识别流程。
图3为非线性频率尺度变换中各滤波器对应的中心频率图。
图4为不同信噪比下的直升机声信号频谱图。
图5为基于本发明得出的不同分段长度下的识别率云图。
图6为基于本发明得出的不同下界频率
Figure 73140DEST_PATH_IMAGE058
时的识别率云图。
图7为基于本发明得出的不同上界频率
Figure 451032DEST_PATH_IMAGE059
时的识别率云图。
图8为基于本发明得出的不同滤波器个数时的识别率云图，非线性频率尺度变换采用Mel尺度，滤波器组采用Mel滤波器组。
图9为基于本发明得出的不同滤波器个数时的识别率云图，非线性频率尺度变换采用Mel尺度，滤波器组采用Gammatone滤波器组。
图10为基于本发明得出的不同滤波器个数时的识别率云图，非线性频率尺度变换采用ERB尺度，滤波器组采用Mel滤波器组。
图11为基于本发明得出的不同滤波器个数时的识别率云图，非线性频率尺度变换采用ERB尺度，滤波器组采用Gammatone滤波器组。
图12为基于本发明得出的不同滤波器个数时的识别率云图，非线性频率尺度变换采用Bark尺度，滤波器组采用Mel滤波器组。
图13为基于本发明得出的不同滤波器个数时的识别率云图，非线性频率尺度变换采用Bark尺度，滤波器组采用Gammatone滤波器组。
图14为尺度变换的最小频率分辨率及识别准确率随滤波器个数的变化图。
图15为尺度变换的频率分辨率随滤波器序号的变化图。
图16 为Mel滤波器和Gammatone滤波器的系数曲线随频率的变化图。
具体实施方式
下面结合附图对本发明作进一步描述。
如图1所示，基于传统谱特征提取的直升机声信号识别流程。对分段加窗的时域声信号数据进行快速傅里叶变换（FFT），并通过对数压缩得到声压级频谱，再求取各分段数据的频谱均值得到待识别的谱特征，最后送入分类器识别目标的类型，这种方法得到的频谱具有单一的分辨率
Figure 864827DEST_PATH_IMAGE060
，
Figure 328169DEST_PATH_IMAGE061
为声信号采样频率，L为分段声信号长度。
如图2所示，一种基于听觉谱特征提取的直升机声信号识别方法，包括如下步骤：
步骤1：分段加窗：将传声器采集的直升机原始声信号划分为长度为L的若干信号段，得到分段声信号
Figure 734880DEST_PATH_IMAGE062
，
Figure 967278DEST_PATH_IMAGE063
；将所述分段声信号
Figure 50509DEST_PATH_IMAGE062
乘以窗函数
Figure 1148DEST_PATH_IMAGE064
以减小分段信号首尾不连续引起的“频谱泄漏”，得到加窗声信号
Figure 883653DEST_PATH_IMAGE065
；为了避免相邻数据段之间信号的变化过大，通常设置50%的重叠区域；
步骤2：FFT分析：对步骤1中得到的加窗声信号
Figure 564033DEST_PATH_IMAGE006
进行FFT分析，得到加窗声信号频谱，进而取加窗声信号频谱绝对值得到幅值谱；
步骤3：尺度变换：在频率分析范围
Figure 54052DEST_PATH_IMAGE007
内，计算听觉滤波器组在非线性频率尺度变换下的中心频率
Figure 491986DEST_PATH_IMAGE066
，
Figure 178182DEST_PATH_IMAGE008
；
Figure 447490DEST_PATH_IMAGE009
为频率分析范围的下界频率，大于FFT分析的最低频率；
Figure 560939DEST_PATH_IMAGE010
为频率分析范围的上界频率，小于FFT分析的最高频率及奈奎斯特频率
Figure 486170DEST_PATH_IMAGE067
的最小值，
Figure 287642DEST_PATH_IMAGE012
为声信号采样频率；所述听觉滤波器组包括Mel滤波器组或Gammatone 滤波器组；所述非线性尺度变换包括Mel尺度变换或ERB尺度变换或Bark尺度变换；
步骤4：听觉滤波：先根据步骤3求解得到的中心频率f i 确定听觉滤波器的表达式，再用听觉滤波器组对步骤2中得到的幅值谱进行带通滤波处理，输出滤波处理后的幅值谱；
步骤5：对数压缩：对步骤4中的听觉滤波器组输出的幅值谱进行对数压缩，得到听觉谱的声压级表示；
步骤6：求取均值：求取步骤5中得到的对数压缩结果的均值得到最终的听觉谱特征；
步骤7：分类识别：按照上述步骤1-6分别求取训练集和测试集的听觉谱特征，先将训练集的听觉谱特征送入分类器进行训练，再将测试集的听觉谱特征送入分类器进行识别，确定直升机类型。
所述窗函数
Figure 614718DEST_PATH_IMAGE014
包括Hamming窗；所述窗函数
Figure 430227DEST_PATH_IMAGE014
为：
Figure 842754DEST_PATH_IMAGE068
，
n为与分段声信号
Figure 870753DEST_PATH_IMAGE017
对应的离散时间点，
Figure 865385DEST_PATH_IMAGE018
；L为分段声信号的长度。
所述Mel尺度变换与中心频率
Figure 320637DEST_PATH_IMAGE069
的对应关系为：
Figure 954880DEST_PATH_IMAGE070
；
其中，
Figure 583308DEST_PATH_IMAGE021
为频率分析范围的下界频率，
Figure 884976DEST_PATH_IMAGE071
为频率分析范围的上界频率，i表示滤波器的序号，
Figure 822714DEST_PATH_IMAGE023
，N为滤波器个数；
Figure 944254DEST_PATH_IMAGE072
表示中心频率
Figure 314055DEST_PATH_IMAGE073
对应的Mel尺度；
Figure 798126DEST_PATH_IMAGE074
为分析范围下界频率对应的Mel尺度，
Figure 595181DEST_PATH_IMAGE027
为分析范围上界频率对应的Mel尺度；
所述ERB尺度变换与中心频率
Figure 204017DEST_PATH_IMAGE075
的对应关系为：
Figure 659400DEST_PATH_IMAGE076
；
其中，
Figure 670081DEST_PATH_IMAGE077
为频率分析范围的下界频率，
Figure 700354DEST_PATH_IMAGE078
为频率分析范围的上界频率，i表示滤波器的序号，
Figure 796486DEST_PATH_IMAGE032
，N为滤波器个数；
Figure 508090DEST_PATH_IMAGE079
表示中心频率
Figure 956301DEST_PATH_IMAGE080
对应的ERB尺度；
Figure 95158DEST_PATH_IMAGE081
为分析范围下界频率对应的ERB尺度，
Figure 678587DEST_PATH_IMAGE082
为分析范围上界频率对应的ERB尺度；
所述 Bark尺度变换与中心频率
Figure 990619DEST_PATH_IMAGE083
的对应关系为：
Figure 710313DEST_PATH_IMAGE084
；
其中，
Figure 20072DEST_PATH_IMAGE085
为频率分析范围的下界频率，
Figure 903846DEST_PATH_IMAGE086
为频率分析范围的上界频率，i表示滤波器的序号，
Figure 957252DEST_PATH_IMAGE041
，N为滤波器个数；
Figure 859349DEST_PATH_IMAGE042
表示中心频率
Figure 340009DEST_PATH_IMAGE087
对应的Bark尺度；
Figure 898029DEST_PATH_IMAGE044
为分析范围下界频率对应的Bark尺度，
Figure 801132DEST_PATH_IMAGE088
为分析范围上界频率对应的Bark尺度。
由上述公式可知，利用非线性频率尺度变换与实际频率
Figure 495419DEST_PATH_IMAGE089
的关系式进行逆变换即可反算出相应尺度下的中心频率
Figure 412559DEST_PATH_IMAGE090
；非线性频率尺度变换与实际频率
Figure 520192DEST_PATH_IMAGE091
的关系式如下：
Mel尺度变换与实际频率
Figure 915402DEST_PATH_IMAGE092
的关系式
Figure 729774DEST_PATH_IMAGE093
为：
Figure 99707DEST_PATH_IMAGE094
；
ERB尺度变换与实际频率
Figure 366740DEST_PATH_IMAGE095
的关系式
Figure 627957DEST_PATH_IMAGE096
为：
Figure 296836DEST_PATH_IMAGE097
；
Bark尺度变换与实际频率
Figure 290199DEST_PATH_IMAGE098
的关系式
Figure 621692DEST_PATH_IMAGE099
为：
Figure 358704DEST_PATH_IMAGE100
。
以频率分析范围[0, FS/2]（FS=44100 Hz为直升机声信号的采样频率）、100个滤波器为例（每隔5个显示1个中心频率）求出了三种尺度变换下，每个滤波器对应的中心频率，如图3所示；由图3可知，尺度变换后的感知频率均与普通频率呈非线性关系，三种尺度变换具有不同的频率特性：在相同的滤波器编号下，Mel尺度对应的中心频率最高但变化更平缓，而 Bark尺度对应的中心频率较低但对应高频段的滤波器更陡峭。
需要说明的是，为了模拟人耳动态、非线性的冲激响应和幅频特性，采用Mel滤波器组或Gammatone滤波器组，是常用的带通听觉滤波器组。
所述Mel滤波器的传递函数为：
Figure 882090DEST_PATH_IMAGE101
；
其中，i表示滤波器的序号，i=1，…，N，N为滤波器个数；
Figure 108672DEST_PATH_IMAGE102
表示Mel滤波器的传递函数；
Figure 350297DEST_PATH_IMAGE103
为某一分析频率，
Figure 625421DEST_PATH_IMAGE104
，
Figure 81941DEST_PATH_IMAGE105
。
所述Gammatone滤波器的时域表达式为：
Figure 417107DEST_PATH_IMAGE106
；
其中，
Figure 208346DEST_PATH_IMAGE107
为Gammatone滤波器的时域表达式，A为滤波器增益；m为滤波器阶数；
Figure 287160DEST_PATH_IMAGE108
为滤波器带宽，
Figure 785138DEST_PATH_IMAGE109
为中心频率，t为时间；需要说明的是，当m=4时更能模拟人耳的听觉特性。
所述分类器采用基于欧式距离的最近邻分类器，但是不限于最近邻分类器。
所述非线性频率尺度变换与听觉滤波器组组合得到听觉谱特征提取算法，包括M-M算法、E-M算法、B-M算法、M-G算法、E-G算法和B-G算法。
需要说明的是，M-M算法是指Mel尺度变换和Mel滤波器组成的听觉谱特征提取算法；E-M算法是指ERB尺度变换和Mel滤波器组成的听觉谱特征提取算法；B-M算法是指Bark尺度变换和Mel滤波器组成的听觉谱特征提取算法；M-G算法是指Mel尺度变换和Gammatone滤波器组成的听觉谱特征提取算法；E-G算法是指ERB尺度变换和Gammatone滤波器组成的听觉谱特征提取算法；B-G算法是指Bark尺度变换和Gammatone滤波器组成的听觉谱特征提取算法。
所述分段声信号的长度L需使得FFT分析的频率分辨率
Figure 868369DEST_PATH_IMAGE060
与直升机的基频为同一数量级；
Figure 819007DEST_PATH_IMAGE110
为声信号采样频率，L为分段声信号长度。
所述听觉滤波器组中的滤波器个数N需使得尺度变换的最小频率分辨率与直升机的基频为同一数量级。
需要说明的是，尺度变换的最小频率分辨率与直升机的基频同一数量级还不够，需要与基频基本相当或还需要略小一些。如基频为20Hz，则最小频率分辨率需在10-25Hz左右。
实施例1
采用上述基于听觉谱特征提取的直升机声信号识别方法，考察在噪声环境下的有效性和鲁棒性。
首先，利用传声器采集得到10种不同型号直升机在野外飞行时主桨、尾桨和发动机等辐射的声信号数据，采样率和采样时间为HS=44100Hz和t=18s；再将每种型号直升机的声信号分为互不重叠的1s时间段的数据，得到共180段数据样本。
其次，从总样本中随机挑选50%数据作为训练集，剩下的50%数据作为测试集，即训练集和测试集的样本数均为90。
最后，在总样本中添加不同强度的高斯白噪声干扰，得到如图4所示的某段直升机声信号在不同信噪比下的功率谱密度图。
由图4可知，直升机声信号的能量集中在 1000 Hz 以下的中低频段，特别是旋翼噪声基频（约 20 Hz）及其谐波对应频率处具有尖锐的、高幅值的声能量。随着含噪数据的信噪比（Signal-to-Noise Ratio, SNR）从SNR=40 dB降低到SNR=-40dB，声信号频谱的幅值呈整体下降的趋势。当信噪比达到SNR=-20 dB和SNR=-40 dB时，高频段的声能量逐渐超过中低频段，从频谱中不再能明显感知到旋翼噪声基频及其谐波频率，直升机的个性特征逐渐淹没于噪声之中。
实施例2
采用上述基于听觉谱特征提取的直升机声信号识别方法，考察分段长度对直升机声信号识别有效性和鲁棒性的影响程度，得出如图5所示不同分段长度下的识别率云图。
如图5所示，给出了不同分段长度L下的识别率云图。从图5可知，随着信噪比的降低，识别率也逐渐降低；但在所有分段长度L下，SNR≥0dB时，识别准确率均在60%以上。另一方面，随着分段长度的增加，识别率呈增大的趋势；但并不是L越大越好，而是有个合适的中间值——高信噪比下，L=1024～4096时，识别准确率接近100%。其原因在于：分段长度较小时，分段数据中包含的有用信息较小且 FFT的频率分辨率较低，（见表 1，FS/ L值越大，频率分辨率越低），而典型直升机的基频通常只有几十Hz，这导致尺度变换前的频谱特性并不能有效地区分不同的直升机；尽管非线性尺度变换和听觉滤波能改善分辨率，但 FFT基础较差时亦不能明显提升识别率。另一方面，分段长度较大时，分段数据更容易由近似平稳变为非平稳，进而影响特征提取和识别性能。当L的长度使得频率分辨率FS/L与直升机基频同一量级时——即几十Hz，识别率较高。
表1，不同数据分段长度下的频率分辨率
L 频率分辨率
64 689.06
128 344.53
256 172.27
512 86.13
1024 43.07
2048 21.53
4096 10.77
8192 5.38
16384 2.69
需要说明的是，实施例2中非线性频率尺度变换采用Mel尺度，滤波器组采用Mel滤波器组。
实施例3
采用上述基于听觉谱特征提取的直升机声信号识别方法，考察频率分析范围对直升机声信号识别有效性和鲁棒性的影响程度，得出如图6-7所示不同频率分析范围下的识别率云图。
如图6-7所示，给出了不同频率分析范围的下界
Figure 701513DEST_PATH_IMAGE058
和上界
Figure 381893DEST_PATH_IMAGE111
设置下的识别率云图。从图中可以看出，频率下界
Figure 324441DEST_PATH_IMAGE112
对识别率的影响大于频率上界
Figure 309846DEST_PATH_IMAGE113
，尽管其最大参数变化范围2560 Hz远小于后者的最大参数变化范围12800 Hz。这印证了直升机声信号识别的主要特征是能量大、衰减慢、传播远的强中低频信号。从图 6中还可以看出，随着频率下界
Figure 730463DEST_PATH_IMAGE112
的增大，识别率呈降低的趋势，特别是
Figure 203032DEST_PATH_IMAGE114
大于1000 Hz后，识别率的下降非常明显。另一方面，随着频率上界
Figure 378799DEST_PATH_IMAGE115
的减小，识别率略有提升。其原因在于：当频率分析的下界
Figure 304030DEST_PATH_IMAGE112
增大时，中低频段的信号特征将被尺度变换和听觉滤波丢弃，数据中所蕴含的直升机个性特征逐渐减小，因此识别率下降；而频率分析的上界
Figure 528337DEST_PATH_IMAGE116
减小时，被丢弃的是容易受干扰影响的高频信号，而保留了中低频的有效信号，因此识别率并不会受太大的影响，反而略有提高。
需要说明的是，实施例3中非线性频率尺度变换采用Mel尺度，滤波器组采用Mel滤波器组。
实施例4
采用上述基于听觉谱特征提取的直升机声信号识别方法，考察滤波器个数对直升机声信号识别有效性和鲁棒性的影响程度，得出如图8-13所示，不同滤波器个数时，应用不同尺度变换和滤波器组的识别率云图。
由图8-13可知，识别准确率均随着滤波器个数的增加呈先迅速提高再逐渐平稳的趋势。
如图14所示，为尺度变换的最小频率分辨率及识别准确率随滤波器个数的变化情况，尺度变换的最小频率分辨率与识别率随滤波器个数的变化趋势刚好相反：随着滤波器个数从1增加到200，尺度变换的频率分辨率从 3000 Hz迅速减小到 100 Hz以内，并在很大的滤波器个数范围内稳定在几十Hz量级。这与听觉谱特征提取中引入非线性尺度变换的初衷——加强直升机声信号的中低频分辨率——是相符的：滤波器个数较少时，尺度变换自身的频率分辨率也较差，其与FFT频谱耦合时并不能提升频谱分析的分辨率；同时，滤波器个数少时，听觉滤波后得到的用于分类识别的听觉谱特征也少，因此识别准确率明显较低。另一方面，滤波器个数足够多时，尺度变换提升中低频分辨率及听觉滤波加强中低频分析的能力逐渐显现，因此识别率增加。
由图8-13还可以看出，整体而言，Bark尺度的识别率略优于ERB尺度，两者均优于Mel尺度；同时，Mel滤波器组的识别率优于Gammatone滤波器组。
如图15所示，为尺度变换的频率分辨率随滤波器序号的变化情况，滤波器个数N=100，由图可知在相同频率分析范围和滤波器个数时，Bark尺度以牺牲高频段的分辨率（滤波器序号大时分辨率值大）为代价，使得中低频分辨率低于ERB尺度和Mel尺度。Mel尺度的分辨率曲线整体比较平缓但多数情况下分辨率低于其它两个尺度变换。
如图16所示，为Mel滤波器和Gammatone滤波器的系数曲线随频率的变化情况，滤波器个数N=100，由图可知，三角形的Mel带通滤波器只在其相邻的两个滤波器中心频率之间具有大于0的系数，而类似gamma函数的Gammatone带通滤波器在其中心频率处具有最大的系数、中心频率两侧的系数曲线较陡但拖尾较长。更重要的是，在滤波器序号相同的情况下，Mel滤波器系数曲线更窄，从而使其具有更尖锐的频率选择能力。
需要说明的是，本发明提出了一种基于听觉谱特征提取的直升机声信号识别方法，将听觉模型引入谱特征提取中以改进其性能；基于三种非线性频率尺度变换（Mel尺度、ERB尺度和 Bark尺度）和两种听觉滤波器组（Mel滤波器和Gammatone滤波器），给出了六种具体的听觉谱特征提取算法：M-M算法、E-M算法、B-M算法、M-G算法、E-G算法和B-G算法。数值仿真实验验证了所提出算法的有效性和鲁棒性，并为算法的参数设置和进一步优化提供了指导：
（1）恰当选取数据分段长度 L和滤波器组个数N以使得频率分辨率与直升机的基频（约几十Hz）同一量级时，听觉谱特征提取的识别率和噪声鲁棒性较高；
（2）频率分析范围
Figure 184576DEST_PATH_IMAGE117
的下界
Figure 468927DEST_PATH_IMAGE118
对识别率和噪声鲁棒性的影响大于上界
Figure 881454DEST_PATH_IMAGE119
，一般可设置
Figure 706190DEST_PATH_IMAGE120
、
Figure 887773DEST_PATH_IMAGE121
（信号的Nyquist频率）以减少参数；
（3）设计更接近于真实人耳听觉感知的计算模型，增强尺度变换的中低频分辨率以及听觉滤波器组的频率选择尖锐性，有助于提升识别性能。
本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程步骤或任何新的组合。如果本领域技术人员，在不脱离本发明的精神所做的非实质性改变或改进，都应该属于本发明权利要求保护的范围。

Claims (7)

1.一种基于听觉谱特征提取的直升机声信号识别方法，其特征在于，包括如下步骤：
步骤1：分段加窗：将传声器采集的直升机原始声信号划分为长度为L的若干信号段，得到分段声信号
Figure 854766DEST_PATH_IMAGE001
，
Figure 968215DEST_PATH_IMAGE002
；将所述分段声信号
Figure 893446DEST_PATH_IMAGE003
乘以窗函数
Figure 445650DEST_PATH_IMAGE004
，得到加窗声信号
Figure 507147DEST_PATH_IMAGE005
；
步骤2：FFT分析：对步骤1中得到的加窗声信号
Figure 791498DEST_PATH_IMAGE006
进行FFT分析，得到加窗声信号频谱，进而取加窗声信号频谱绝对值得到幅值谱；
步骤3：尺度变换：在频率分析范围
Figure 204025DEST_PATH_IMAGE007
内，计算听觉滤波器组在非线性频率尺度变换下的中心频率
Figure 232023DEST_PATH_IMAGE008
，
Figure 413606DEST_PATH_IMAGE009
；
Figure 150749DEST_PATH_IMAGE010
为频率分析范围的下界频率，大于FFT分析的最低频率；
Figure 50572DEST_PATH_IMAGE011
为频率分析范围的上界频率，小于FFT分析的最高频率及奈奎斯特频率
Figure 616683DEST_PATH_IMAGE012
的最小值，
Figure 918351DEST_PATH_IMAGE013
为声信号采样频率；所述听觉滤波器组包括Mel滤波器组或Gammatone滤波器组；所述非线性频率尺度变换包括Mel尺度变换或ERB尺度变换或Bark尺度变换；
步骤4：听觉滤波：先根据步骤3求解得到的中心频率
Figure 278925DEST_PATH_IMAGE014
确定听觉滤波器的表达式，再用听觉滤波器组对步骤2中得到的幅值谱进行带通滤波处理，输出滤波处理后的幅值谱；
步骤5：对数压缩：对步骤4中的听觉滤波器组输出的幅值谱进行对数压缩；
步骤6：求取均值：求取步骤5中得到的对数压缩结果的均值得到最终的听觉谱特征；
步骤7：分类识别：按照上述步骤1-6分别求取训练集和测试集的听觉谱特征；先将训练集的听觉谱特征送入分类器进行训练，再将测试集的听觉谱特征送入分类器进行识别，确定直升机类型；
所述分段声信号的长度L需使得FFT分析的频率分辨率
Figure 728361DEST_PATH_IMAGE015
与直升机的基频为同一数量级；
Figure 98162DEST_PATH_IMAGE016
为声信号采样频率，L为分段声信号长度；所述听觉滤波器组中的滤波器个数N需使得尺度变换的最小频率分辨率与直升机的基频为同一数量级。
2.如权利要求1所述的一种基于听觉谱特征提取的直升机声信号识别方法，其特征在于，所述窗函数
Figure 988758DEST_PATH_IMAGE017
包括Hamming窗；窗函数
Figure 785813DEST_PATH_IMAGE018
为：
Figure 394649DEST_PATH_IMAGE019
，
n为与分段声信号
Figure 643840DEST_PATH_IMAGE020
对应的离散时间点，
Figure 654521DEST_PATH_IMAGE021
；L为分段声信号的长度。
3.如权利要求1所述的一种基于听觉谱特征提取的直升机声信号识别方法，其特征在于，所述Mel尺度变换与中心频率
Figure 622477DEST_PATH_IMAGE022
的对应关系为：
Figure 718609DEST_PATH_IMAGE023
；
其中，
Figure 430213DEST_PATH_IMAGE024
为频率分析范围的下界频率，
Figure 92139DEST_PATH_IMAGE025
为频率分析范围的上界频率，i表示滤波器的序号，
Figure 230996DEST_PATH_IMAGE026
，N为滤波器个数；
Figure 814424DEST_PATH_IMAGE027
表示中心频率
Figure 64140DEST_PATH_IMAGE028
对应的Mel尺度；
Figure 49413DEST_PATH_IMAGE029
为分析范围下界频率对应的Mel尺度，
Figure 172221DEST_PATH_IMAGE030
为分析范围上界频率对应的Mel尺度；
所述ERB尺度变换与中心频率
Figure 242945DEST_PATH_IMAGE031
的对应关系为：
Figure 296352DEST_PATH_IMAGE032
；
其中，
Figure 136132DEST_PATH_IMAGE033
为频率分析范围的下界频率，
Figure 679109DEST_PATH_IMAGE034
为频率分析范围的上界频率，i表示滤波器的序号，
Figure 971550DEST_PATH_IMAGE035
，N为滤波器个数；
Figure 828648DEST_PATH_IMAGE036
表示中心频率
Figure 522934DEST_PATH_IMAGE037
对应的ERB尺度；
Figure 174495DEST_PATH_IMAGE038
为分析范围下界频率对应的ERB尺度，
Figure 32861DEST_PATH_IMAGE039
为分析范围上界频率对应的ERB尺度；
所述 Bark尺度变换与中心频率
Figure 428070DEST_PATH_IMAGE040
的对应关系为：
Figure 976863DEST_PATH_IMAGE041
；
其中，
Figure 64905DEST_PATH_IMAGE042
为频率分析范围的下界频率，
Figure 331938DEST_PATH_IMAGE043
为频率分析范围的上界频率，i表示滤波器的序号，
Figure 858735DEST_PATH_IMAGE044
，N为滤波器个数；
Figure 262034DEST_PATH_IMAGE045
表示中心频率
Figure 520977DEST_PATH_IMAGE046
对应的Bark尺度；
Figure 275307DEST_PATH_IMAGE047
为分析范围下界频率对应的Bark尺度，
Figure 12318DEST_PATH_IMAGE048
为分析范围上界频率对应的Bark尺度。
4.如权利要求1所述的一种基于听觉谱特征提取的直升机声信号识别方法，其特征在于，所述Mel滤波器的传递函数为：
Figure 535704DEST_PATH_IMAGE049
；
其中，i表示滤波器的序号，i=1，…，N，N为滤波器个数；
Figure 513018DEST_PATH_IMAGE050
表示Mel滤波器的传递函数；
Figure 489064DEST_PATH_IMAGE051
为某一分析频率，
Figure 29767DEST_PATH_IMAGE052
，
Figure 673238DEST_PATH_IMAGE053
。
5.如权利要求1所述的一种基于听觉谱特征提取的直升机声信号识别方法，其特征在于，所述Gammatone滤波器的时域表达式为：
Figure 8404DEST_PATH_IMAGE054
；
其中，
Figure 799643DEST_PATH_IMAGE055
为Gammatone滤波器的时域表达式，A为滤波器增益；m为滤波器阶数；
Figure 878457DEST_PATH_IMAGE056
为滤波器带宽；
Figure 376435DEST_PATH_IMAGE057
为中心频率，t为时间。
6.如权利要求1所述的一种基于听觉谱特征提取的直升机声信号识别方法，其特征在于，所述分类器包括基于欧式距离的最近邻分类器。
7.如权利要求1所述的一种基于听觉谱特征提取的直升机声信号识别方法，其特征在于，所述非线性频率尺度变换与听觉滤波器组组合得到听觉谱特征提取算法，包括M-M算法、E-M算法、B-M算法、M-G算法、E-G算法和B-G算法；所述M-M算法是指Mel尺度变换和Mel滤波器组成的听觉谱特征提取算法；所述E-M算法是指ERB尺度变换和Mel滤波器组成的听觉谱特征提取算法；所述B-M算法是指Bark尺度变换和Mel滤波器组成的听觉谱特征提取算法；所述M-G算法是指Mel尺度变换和Gammatone滤波器组成的听觉谱特征提取算法；所述E-G算法是指ERB尺度变换和Gammatone滤波器组成的听觉谱特征提取算法；所述B-G算法是指Bark尺度变换和Gammatone滤波器组成的听觉谱特征提取算法。
CN202110951550.1A 2021-08-19 2021-08-19 一种基于听觉谱特征提取的直升机声信号识别方法 Active CN113409819B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110951550.1A CN113409819B (zh) 2021-08-19 2021-08-19 一种基于听觉谱特征提取的直升机声信号识别方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110951550.1A CN113409819B (zh) 2021-08-19 2021-08-19 一种基于听觉谱特征提取的直升机声信号识别方法

Publications (2)

Publication Number Publication Date
CN113409819A CN113409819A (zh) 2021-09-17
CN113409819B true CN113409819B (zh) 2022-01-25

Family

ID=77688545

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110951550.1A Active CN113409819B (zh) 2021-08-19 2021-08-19 一种基于听觉谱特征提取的直升机声信号识别方法

Country Status (1)

Country Link
CN (1) CN113409819B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101829689A (zh) * 2010-03-31 2010-09-15 北京科技大学 一种基于声信号的热轧带钢甩尾故障识别方法
CN104157290A (zh) * 2014-08-19 2014-11-19 大连理工大学 一种基于深度学习的说话人识别方法
CN104778948A (zh) * 2015-04-29 2015-07-15 太原理工大学 一种基于弯折倒谱特征的抗噪语音识别方法
CN110189757A (zh) * 2019-06-27 2019-08-30 电子科技大学 一种大熊猫个体识别方法、设备及计算机可读存储介质

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103714810B (zh) * 2013-12-09 2016-05-25 西北核技术研究所 基于Gammatone滤波器组的车型特征提取方法
EP3790006A4 (en) * 2018-06-29 2021-06-09 Huawei Technologies Co., Ltd. VOICE COMMAND PROCESS, PORTABLE DEVICE AND TERMINAL
CA3081962A1 (en) * 2019-06-07 2020-12-07 PredictMedix Inc. Systems and methods for detecting impairment of an individual
CN112086105B (zh) * 2020-08-31 2022-08-19 中国船舶重工集团公司七五0试验场 一种基于Gammatone分频带连续谱特征的目标识别方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101829689A (zh) * 2010-03-31 2010-09-15 北京科技大学 一种基于声信号的热轧带钢甩尾故障识别方法
CN104157290A (zh) * 2014-08-19 2014-11-19 大连理工大学 一种基于深度学习的说话人识别方法
CN104778948A (zh) * 2015-04-29 2015-07-15 太原理工大学 一种基于弯折倒谱特征的抗噪语音识别方法
CN110189757A (zh) * 2019-06-27 2019-08-30 电子科技大学 一种大熊猫个体识别方法、设备及计算机可读存储介质

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"基于Gammatone倒谱系数的直升机声信号识别";王勇;《湖南大学学报（自然科学版）》;20210625;第48卷(第6期);第77-79页 *

Also Published As

Publication number Publication date
CN113409819A (zh) 2021-09-17

Similar Documents

Publication Publication Date Title
CN111261189B (zh) 一种车辆声音信号特征提取方法
CN106653032A (zh) 低信噪比环境下基于多频带能量分布的动物声音检测方法
CN113707176B (zh) 一种基于声信号及深度学习技术的变压器故障检测方法
CN110570880B (zh) 一种鼾声信号识别方法
CN110970042B (zh) 一种电子听诊器的肺部啰音人工智能实时分类方法、系统、装置及可读存储介质
CN106772331A (zh) 目标识别方法和目标识别装置
CN113763986B (zh) 一种基于声音分类模型的空调内机异常声音检测方法
CN114863937A (zh) 基于深度迁移学习与XGBoost的混合鸟鸣识别方法
CN113409819B (zh) 一种基于听觉谱特征提取的直升机声信号识别方法
CN109300486B (zh) 基于PICGTFs和SSMC增强的腭裂语音咽擦音自动识别方法
CN113345443A (zh) 基于梅尔频率倒谱系数的海洋哺乳动物发声检测识别方法
CN112735468A (zh) 一种基于mfcc的汽车座椅电机异常噪声检测方法
CN112652312A (zh) 声纹相似度智能识别系统、方法及存储介质
CN114743562B (zh) 一种飞机声纹识别方法、系统、电子设备及存储介质
CN114626420A (zh) 一种基于CWT-Xception-RF的φ-OTDR振动事件分类方法
CN116597853A (zh) 一种音频消噪方法
CN112201226B (zh) 一种发声方式判别方法及系统
CN114093385A (zh) 一种无人机检测方法及装置
CN115331678A (zh) 利用Mel频率倒谱系数的广义回归神经网络声信号识别方法
CN112233693B (zh) 一种音质评估方法、装置和设备
CN113053351A (zh) 一种基于听觉感知的飞机舱内噪声合成方法
CN112086105A (zh) 一种基于Gammatone分频带连续谱特征的目标识别方法
Sánchez-García et al. A novel image-processing based method for the automatic detection, extraction and characterization of marine mammal tonal calls
CN112908344A (zh) 一种鸟鸣声智能识别方法、装置、设备和介质
CN111920390A (zh) 一种基于嵌入式终端的鼾声检测方法

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
GR01 Patent grant
GR01 Patent grant