最全语音特征总结

1、能量
音频的能量通常指的是时域上每帧的能量，幅度的平方。在简单的语音活动检测（Voice Activity Detection，VAD）中，直接利用能量特征：能量大的音频片段是语音，能量小的音频片段是非语音（包括噪音、静音段等）。这种VAD的局限性比较大，正确率也不高，对噪音非常敏感。
2、短时能量
短时能量体现的是信号在不同时刻的强弱程度。
3、声强和声强级（声压和声压级）
单位时间内通过垂直于声波传播方向的单位面积的平均声能，称作声强，声强用P表示，单位为“瓦/平米”。实验研究表明，人对声音的强弱感觉并不是与声强成正比，而是与其对数成正比，所以一般声强用声强级来表示：
4、过零率
过零率体现的是信号过零点的次数，体现的是频率特性。
5、基频、基音周期、音高–有代码
基音周期反映了声门相邻两次开闭之间的时间间隔，基频（fundamental frequency， F0）则是基音周期的倒数，对应着声带振动的频率，代表声音的音高，声带振动越快，基频越高。它是语音激励源的一个重要特征，比如可以通过基频区分性别。一般来说，成年男性基频在 100-250Hz左右，成年女性基频在 150-350Hz左右，女声的音高一般比男声稍高。 人类可感知声音的频率大致在20-20000Hz之间，人类对于基频的感知遵循对数律，也就是说，人们会感觉100Hz到200Hz的差距，与200Hz到400Hz的差距相同。因此，音高常常用基频的对数来表示。
6、共振峰
声门处的准周期激励进入声道时会引起共振特性，产生一组共振频率，这一组共振频率称为共振峰频率或简称共振峰。共振峰包含在语音的频谱包络中，频谱极大值就是共振峰。频率最低的共振峰称为第一共振峰，对应的频率也称作基频，决定语音的F0，其它的共振峰统称为谐波，
7、语速（speaking rate）
特征表达了讲话速度的快慢，可以定义为单位时间内发音的词汇（或者音节）个数。语速受文化、环境、思维和表达能力多种因素的影响。和语速密切相关的因素还有停顿，是否考虑语段中的停顿对语速的计算数值有明显影响。
8、谱特征（spectral feature）
8.1、MFCC
MFCC，即梅尔倒谱系数（Mel-scaleFrequency Cepstral Coefficients）。是一种非线性映射，根据人耳对不同频率的声波有不同的听觉敏感度进行映射的。
原理：根据人耳听觉机理的研究发现，人耳对不同频率的声波有不同的听觉敏感度。从200HZ到5000HZ对语音的清晰度影响最大。两个响度不等的声音作用于人耳时，则响度较高的频率成分的存在会影响到对响度较低的频率成分的成分，使其变得不易察觉，这种现象称为掩蔽效应。由于频率较低的声音在内耳蜗基底膜上行波传递的距离大于频率较高的声音，故一般来说，低音容易掩蔽高音，而高音掩蔽低音较困难。在低频处的声音掩蔽临界带宽较高频要小。所以从低频到高频这一频带内按临界带宽的大小由密到疏安排一组带通滤波器，对输入信号进行滤波。
梅尔倒谱系数是在Mel标度频率域上提取出来的参数。
用途：可以看到，MFCC是根据掩蔽效应原理，在Mel标度上提取出来的参数，为了符合人耳听觉机理的。所以其常用于ASR上，语音识别中。当然，其它一些处理也会考虑到直接用线性谱运算量太大，然后采用MFCC。
具体特征提取过程：语音特征参数MFCC提取过程详解_James Zhang’s Blog-CSDN博客_mfcc特征提取后的结果
8.2、Bark谱
Bark谱与MFCC，Mel谱非常相似，都是将线性谱映射到非线性谱上的表征，而且都是低频带宽低，高频带宽高。但还是略有区别的：
上世纪，研究者发现人耳结构对24个频点产生共振，根据这一理论，Eberhard Zwicker在1961年针对人耳特殊结构提出：信号在频带上也呈现出24个临界频带，分别从1到24。这就是Bark域。
其实Mel谱和Bark谱两者的核心都是在掩蔽效应，人耳对不同频带听感不同。然后划分出的非线性表示。
用途：比较多的看到，Bark谱用于基频，降噪，编解码，特殊声音检测等领域。
8.3、CQT
CQT即恒Q变换，它是用一组恒Q滤波器对时域语音信号进行滤波，因为，滤波器是恒Q的，即中心频率与带宽比相同，则在低频时，带宽窄，高频时带宽高，从而得到非线性频域信号。
与MFCC，Bark谱非常相似，也是一种将线性谱转换到非线性谱的处理。
CQT更加符合乐理，在音乐中，所有的音都是由若干八度的12平均律共同组成的，12个半音等于一个八度，一个八度的跨度等于频率翻倍，所以一个半音等于 21/12
倍频。因此，音乐中的音调呈指数型跨度的，而CQT就很好的模拟了这种非线性度，以 log2
为底的非线性频谱。
用途：常用于音乐方向。但因深度学习的兴起，很多方向也会用这种特征。
9、特征集汇总
9.1 GeMAPS特征集
GeMAPS特征集总共62个特征，这62个都是HSF特征，是由18个LLD特征计算得到。下面先介绍18个LLD特征，然后介绍62个HSF特征。这里只简单介绍每个特征的概念，不涉及具体计算细节。

18个LLD特征包括6个频率相关特征，3个能量/振幅相关特征，9个谱特征。
6个频率相关特征包括：Pitch（log F0，在半音频率尺度上计算，从27.5Hz开始）；Jitter（单个连续基音周期内的偏差，偏差衡量的是观测变量与特定值的差，如果没有指明特定值通常使用的是变量的均值）；前三个共振峰的中心频率，第一个共振峰的带宽。3个能量/振幅的特征包括：Shimmer（相邻基音周期间振幅峰值之差），Loudness（从频谱中得到的声音强度的估计，可以根据能量来计算），HNR（Harmonics-to-noise）信噪比。9个谱特征包括，Alpha Ratio（50-1000Hz的能量和除以1-5kHz的能量和），Hammarberg Index（0-2kHz的最强能量峰除以2-5kHz的最强能量峰），Spectral Slope 0-500 Hz and 500-1500 Hz（对线性功率谱的两个区域0-500 Hz和500-1500 Hz做线性回归得到的两个斜率），Formant 1, 2, and 3 relative energy（前三个共振峰的中心频率除以基音的谱峰能量），Harmonic difference H1-H2（第一个基音谐波H1的能量除以第二个基音谐波的能量），Harmonic difference H1-A3（第一个基音谐波H1的能量除以第三个共振峰范围内的最高谐波能量）。
对18个LLD做统计，计算的时候是对3帧语音做symmetric moving average。首先计算算术平均和coefficient of variation（计算标准差然后用算术平均规范化），得到36个统计特征。然后对loudness和pitch运算8个函数，20百分位，50百分位，80百分位，20到80百分位之间的range，上升/下降语音信号的斜率的均值和标准差。这样就得到16个统计特征。上面的函数都是对voiced regions（非零的F0）做的。对Alpha Ratio，Hammarberg Index，Spectral Slope 0-500 Hz and 500-1500 Hz做算术平均得到4个统计特征。另外还有6个时间特征，每秒loudness峰的个数，连续voiced regions（F0>0）的平均长度和标准差，unvoiced regions（F0=0）的平均长度和标准差，每秒voiced regions的个数。36+16+4+6得到62个特征。
9.2 eGeMAPS特征集
（1）eGeMAPS是GeMAPS的扩展，在18个LLDs的基础上加了一些特征，包括5个谱特征：MFCC1-4和Spectral flux（两个相邻帧的频谱差异）和2个频率相关特征：第二个共振峰和第三个共振峰的带宽。
（2）对这扩展的7个LLDs做算术平均和coefficient of variation（计算标准差然后用算术平均规范化）可以得到14个统计特征。对于共振峰带宽只在voiced region做，对于5个谱特征在voiced region和unvoiced region一起做。
（3）另外，只在unvoiced region计算spectral flux的算术平均，然后只在voiced region计算5个谱特征的算术平均和coefficient of variation，得到11个统计特征。
（4）另外，还加多一个equivalent sound level 。
（5）所以总共得到14+11+1=26个扩展特征，加上原GeMAPS的62个特征，得到88个特征，这88个特征就是eGeMAPS的特征集。

9.3 ComParE特征集
（1）ComParE，Computational Paralinguistics ChallengE，是InterSpeech上的一个挑战赛，从13年至今（2018年），每年都举办，每年有不一样的挑战任务。
（2）从13年开始至今（2018年），ComParE的挑战都会要求使用一个设计好的特征集，这个特征集包含了6373个静态特征，是在LLD上计算各种函数得到的，称为ComParE特征集。
（3）可以通过openSmile开源包来获得，另外前面提到的eGeMAPS也可以用openSmile获得。
五：2009 InterSpeech挑战赛特征
（1）前面说的6373维特征集ComparE是13年至今InterSpeech挑战赛中用的。（2）有论文还用了09年InterSpeech上Emotion Challenge提到的特征，总共有384个特征，计算方法如下。
（3）首先计算16个LLD，过零率，能量平方根，F0，HNR（信噪比，有些论文也叫vp，voice probability 人声概率），MFCC1-12，然后计算这16个LLD的一阶差分，可以得到32个LLD。
（4）对这32个LLD应用12个统计函数，最后得到32x12 = 384个特征。
（5）同样可以通过openSmile来获得。

9.4 BoAW
BoAW，bag-of-audio-words，是特征的进一步组织表示，是根据一个codebook对LLDs做计算得到的。这个codebook可以是k-means的结果，也可以是对LLDs的随机采样。在论文会看到BoAW特征集的说法，指的是某个特征集的BoAW形式。比如根据上下文“使用特征集有ComparE和BoAW”，可以知道，这样的说法其实是指原来的特征集ComparE，和ComparE经过计算后得到的BoAW表示。可通过openXBOW开源包来获得BoAW表示。

9.5 YAAFE
使用YAAFE库提取到的特征，具体特征见YAAFE主页。

部分代码如下：

import argparse
import os, librosa,scipy,csv
import numpy as np


class audio:
    def __init__(self, input_file, sr=None, frame_len=512, n_fft=None, win_step=2 / 3, window="hamming"):
        """
        初始化
        :param input_file: 输入音频文件
        :param sr: 所输入音频文件的采样率，默认为None
        :param frame_len: 帧长，默认512个采样点(32ms,16kHz),与窗长相同
        :param n_fft: FFT窗口的长度，默认与窗长相同
        :param win_step: 窗移，默认移动2/3，512*2/3=341个采样点(21ms,16kHz)
        :param window: 窗类型，默认汉明窗
        """
        self.input_file = input_file
        self.frame_len = frame_len  # 帧长，单位采样点数
        self.wave_data, self.sr = librosa.load(self.input_file, sr=sr)
        self.window_len = frame_len  # 窗长512
        if n_fft is None:
            self.fft_num = self.window_len  # 设置NFFT点数与窗长相等
        else:
            self.fft_num = n_fft
        self.win_step = win_step
        self.hop_length = round(self.window_len * win_step)  # 重叠部分采样点数设置为窗长的1/3（1/3~1/2）,即帧移(窗移)2/3
        self.window = window


    def energy(self):
        """
        每帧内所有采样点的幅值平方和作为能量值
        :return: 每帧能量值，np.ndarray[shape=(1，n_frames), dtype=float64]
        """
        mag_spec = np.abs(librosa.stft(self.wave_data, n_fft=self.fft_num, hop_length=self.hop_length,
                                       win_length=self.frame_len, window=self.window))
        pow_spec = np.square(mag_spec)  # [frequency, time (n_frames)]
        energy = np.sum(pow_spec, axis=0)  # [n_frames]
        energy = np.where(energy == 0, np.finfo(np.float64).eps,
                          energy)  # 避免能量值为0，防止后续取log出错(eps是取非负的最小值), 即np.finfo(np.float64).eps = 2.220446049250313e-16
        return energy

    def short_time_energy(self):
        """
        计算语音短时能量：每一帧中所有语音信号的平方和
        :return: 语音短时能量列表(值范围0-每帧归一化后能量平方和，这里帧长512，则最大值为512)，
        np.ndarray[shape=(1，无加窗，帧移为0的n_frames), dtype=float64]
        """
        energy = []  # 语音短时能量列表
        energy_sum_per_frame = 0  # 每一帧短时能量累加和
        for i in range(len(self.wave_data)):  # 遍历每一个采样点数据
            energy_sum_per_frame += self.wave_data[i] ** 2  # 求语音信号能量的平方和
            if (i + 1) % self.frame_len == 0:  # 一帧所有采样点遍历结束
                energy.append(energy_sum_per_frame)  # 加入短时能量列表
                energy_sum_per_frame = 0  # 清空和
            elif i == len(self.wave_data) - 1:  # 不满一帧，最后一个采样点
                energy.append(energy_sum_per_frame)  # 将最后一帧短时能量加入列表
        energy = np.array(energy)
        energy = np.where(energy == 0, np.finfo(np.float64).eps, energy)  # 避免能量值为0，防止后续取log出错(eps是取非负的最小值)
        return energy

    def intensity(self):
        """
        计算声音强度，用声压级表示：每帧语音在空气中的声压级Sound Pressure Level(SPL)，单位dB
        公式：20*lg(P/Pref)，P为声压（Pa），Pref为参考压力(听力阈值压力)，一般为1.0*10-6 Pa
        这里P认定为声音的幅值：求得每帧所有幅值平方和均值，除以Pref平方，再取10倍lg
        :return: 每帧声压级，dB，np.ndarray[shape=(1，无加窗，帧移为0的n_frames), dtype=float64]
        """
        p0 = 1.0e-6  # 听觉阈限压力auditory threshold pressure: 2.0*10-5 Pa
        e = self.short_time_energy()
        spl = 10 * np.log10(1 / (np.power(p0, 2) * self.frame_len) * e)
        return spl

    def zero_crossing_rate(self):
        """
        计算语音短时过零率：单位时间(每帧)穿过横轴（过零）的次数
        :return: 每帧过零率次数列表，np.ndarray[shape=(1，无加窗，帧移为0的n_frames), dtype=uint32]
        """
        zcr = []  # 语音短时过零率列表
        counting_sum_per_frame = 0  # 每一帧过零次数累加和，即过零率
        for i in range(len(self.wave_data)):  # 遍历每一个采样点数据
            if i % self.frame_len == 0:  # 开头采样点无过零，因此每一帧的第一个采样点跳过
                continue
            if self.wave_data[i] * self.wave_data[i - 1] < 0:  # 相邻两个采样点乘积小于0，则说明穿过横轴
                counting_sum_per_frame += 1  # 过零次数加一
            if (i + 1) % self.frame_len == 0:  # 一帧所有采样点遍历结束
                zcr.append(counting_sum_per_frame)  # 加入短时过零率列表
                counting_sum_per_frame = 0  # 清空和
            elif i == len(self.wave_data) - 1:  # 不满一帧，最后一个采样点
                zcr.append(counting_sum_per_frame)  # 将最后一帧短时过零率加入列表
        return np.array(zcr, dtype=np.uint32)

    def pitch(self, ts_mag=0.25):
        """
        获取每帧音高，即基频，这里应该包括基频和各次谐波，最小的为基频（一次谐波），其他的依次为二次、三次...谐波
        各次谐波等于基频的对应倍数，因此基频也等于各次谐波除以对应的次数，精确些等于所有谐波之和除以谐波次数之和
        :param ts_mag: 幅值倍乘因子阈值，>0，大于np.average(np.nonzero(magnitudes)) * ts_mag则认为对应的音高有效,默认0.25
        :return: 每帧基频及其对应峰的幅值(>0)，
                 np.ndarray[shape=(1 + n_fft/2，n_frames), dtype=float32]，（257，全部采样点数/(512*2/3)+1）

        usage:
        pitches, mags = self.pitch()  # 获取每帧基频
        f0_likely = []  # 可能的基频F0
        for i in range(pitches.shape[1]):  # 按列遍历非0最小值，作为每帧可能的F0
            try:
                f0_likely.append(np.min(pitches[np.nonzero(pitches[:, i]), i]))
            except ValueError:
                f0_likely.append(np.nan)  # 当一列，即一帧全为0时，赋值最小值为nan
        f0_all = np.array(f0_likely)
        """
        mag_spec = np.abs(librosa.stft(self.wave_data, n_fft=self.fft_num, hop_length=self.hop_length,
                                       win_length=self.frame_len, window=self.window))
        pitches, magnitudes = librosa.piptrack(S=mag_spec, sr=self.sr, threshold=1.0, ref=np.mean,
                                               fmin=50, fmax=500)  # 人类正常说话基频最大可能范围50-500Hz
        ts = np.average(magnitudes[np.nonzero(magnitudes)]) * ts_mag
        pit_likely = pitches
        mag_likely = magnitudes
        pit_likely[magnitudes < ts] = 0
        mag_likely[magnitudes < ts] = 0
        return pit_likely.tolist(), mag_likely.tolist()

    def preemphasis(self, coef=0.97, zi=None, return_zf=False):
        """Pre-emphasize an audio signal with a first-order auto-regressive filter:

            y[n] -> y[n] - coef * y[n-1]


        Parameters
        ----------
        y : np.ndarray
            Audio signal

        coef : positive number
            Pre-emphasis coefficient.  Typical values of ``coef`` are between 0 and 1.

            At the limit ``coef=0``, the signal is unchanged.

            At ``coef=1``, the result is the first-order difference of the signal.

            The default (0.97) matches the pre-emphasis filter used in the HTK
            implementation of MFCCs [#]_.

            .. [#] http://htk.eng.cam.ac.uk/

        zi : number
            Initial filter state.  When making successive calls to non-overlapping
            frames, this can be set to the ``zf`` returned from the previous call.
            (See example below.)

            By default ``zi`` is initialized as ``2*y[0] - y[1]``.

        return_zf : boolean
            If ``True``, return the final filter state.
            If ``False``, only return the pre-emphasized signal.

        Returns
        -------
        y_out : np.ndarray
            pre-emphasized signal

        zf : number
            if ``return_zf=True``, the final filter state is also returned
        """
        y=self.wave_data
        b = np.asarray([1.0, -coef], dtype=y.dtype)
        a = np.asarray([1.0], dtype=y.dtype)

        if zi is None:
            # Initialize the filter to implement linear extrapolation
            zi = 2 * y[..., 0] - y[..., 1]

        zi = np.atleast_1d(zi)

        y_out, z_f = scipy.signal.lfilter(b, a, y, zi=np.asarray(zi, dtype=y.dtype))

        if return_zf:
            return y_out, z_f

        return y_out
def get_max_min_avg(data):
    '''

    :param data:
    :return: 最大值，最小值，平均值，方差，标准差，和，中值
    '''
    data=np.array(data)
    # print(data.max(),data.min(),data.mean(),data.var(),data.std(),data.sum(),np.median(data))
    return [data.max(),data.min(),data.mean(),data.var(),data.std(),data.sum(),np.median(data)]

def get_feature(filename):
    '''

    :param filename:
    :return: [文件路径，能量，短时能量，声音强度，过零率， 基频， 各次协波]
    '''
    audio_process=audio(filename)
    pit,mag=audio_process.pitch()
    data = [filename]
    data.extend(get_max_min_avg(audio_process.energy()))
    data.extend(get_max_min_avg(audio_process.short_time_energy()))
    data.extend(get_max_min_avg(audio_process.intensity()))
    data.extend(get_max_min_avg(audio_process.zero_crossing_rate()))
    data.extend(get_max_min_avg(pit))
    data.extend(get_max_min_avg(mag))
    return data

def get_features_by_folder(args):
    csv_result=[]
    for filename in os.listdir(args.audio_path):
        if filename[-3:]=='wav':
            features=get_feature(os.path.join(args.audio_path,filename))
            csv_result.append(features)

    with open(args.csv_path,'w',encoding='utf-8') as fp:
        csv_writer=csv.writer(fp,delimiter='|')
        csv_writer.writerows(csv_result)

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('--csv_path', type=str, default='/home/wangyuke_i/tmp/feature.csv')
    parser.add_argument('--audio_path', type=str, default='/home/florence/data/amr_file_0502_trimed/')
    args = parser.parse_args()

    get_features_by_folder(args)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225