Python实现语音信号的短时平均过零率

冷不防 2022-10-05 05:58 236阅读 0赞

### 短时平均过零率 ###

平均过零率：  
短时平均过零率表示的是一帧语音中语音信号波形穿过横轴也就是零电平的次数，也可以理解为样本数值改变符号的次数。  
若语音信号分帧以后为  y i ( n ) \{y\}\_\{i\} (n) yi(n)，帧长为L，则短时平均过零率可以表示为：  
 Z ( i ) = 1 2 ∑ n = 0 L − 1 ∣ s g n \[ y i ( n ) \] − s g n \[ y i ( n − 1 ) \] ∣ 1 < = i < = f n Z(i) = \\frac\{1\}\{2\} \\sum\_\{n=0\}^\{L-1\} |sgn\[\{y\}\_\{i\}(n)\]-sgn\[\{y\}\_\{i\}(n-1)\]| \\quad 1<=i<=fn Z(i)=21n=0∑L−1∣sgn\[yi(n)\]−sgn\[yi(n−1)\]∣1<=i<=fn  
其中，sgn\[\]是符号函数：  
 s g n \[ x \] = \{ 1 x>=0 − 1 x<0 sgn\[x\]= \\begin\{cases\} 1& \\text\{x>=0\}\\\\ -1& \\text\{x<0\} \\end\{cases\} sgn\[x\]=\{ 1−1x>=0x<0  
因为是过零，当离散信号相邻的取样值改变符号时，他们的乘积就为负数：  
 y i ( n ) y i ( n + 1 ) < 0 \{y\}\_\{i\}(n) \{y\}\_\{i\}(n+1) < 0 yi(n)yi(n\+1)<0

程序代码：

''' 短时平均过零率 '''
    
    import numpy as np
    import wave
    import scipy
    import scipy.signal as signal
    import matplotlib.pyplot as plt  # 导入绘图工作的函数集合
    # 设置字体，使可以输入中文
    plt.rcParams['font.sans-serif']=['STSong']
    
    # 读取音频文件
    f = wave.open('D:\任务\Bluesky3.wav','rb')
    '''
    nchannels = f.getnchannels()  # 返回音频通道的数量（1对于单声道，2对于立体声）---1
    sampwidth = f.getsampwidth()  # 以字节为单位返回样本宽度 ---2
    framerate = f.getframerate()  # 返回采样频率 ---8000
    nframes = f.getnframes()  # 返回音频帧数 ---19000
    print(nchannels)
    print(sampwidth)
    print(framerate)
    print(nframes)
    '''
    # 上面的四个参数也可以通过下面getparams()函数一次性读取出来
    params = f.getparams()  # 一次性返回所有的音频参数，声道数、量化位数、采样频率、采样点数
    print(params)
    nchannels, sampwidth, framerate, nframes = params[:4]
    
    str_data = f.readframes(nframes)  # 指定需要读取的长度(以取样点为单位)，返回字符串类型数据
    # 将读取的字符串数据转换为一维int16类型的数组
    wave_data = np.frombuffer(str_data, dtype=np.int16)
    wave_data = wave_data*1.0/(max(abs(wave_data)))  # wave幅值归一化
    f.close()  # 关闭文件流wave
    # 设置帧长和帧移
    lframe = 200
    mframe = 80
    
    # 定义分帧函数
    def frame(wave_data, lframe, mframe):
        signal_length = len(wave_data)  # 获取语音信号的长度
        if signal_length<=lframe:  # 判断信号的长度是否小于帧长
            fns = 1  # 若小于则帧数为1
        else:
            # 计算总共有多少帧
            fns = (1.0*signal_length-lframe+mframe)/mframe
            fns = int(np.ceil(fns))  # 将帧数向上取整，如果是浮点型则加一
        all_length = int((fns - 1) * mframe + lframe)  # 所有帧加起来总的铺平后的长度
        # 不够的长度使用0填补,填补数为fns帧总长减去信号总长
        fillzeros = np.zeros((all_length - signal_length,))
        # 填补后的信号：原信号加上填补的0构成所有的帧数总信号
        fillsignal = np.concatenate((wave_data, fillzeros))  # concatenate连接两个维度相同的矩阵
        # 对所有帧的时间点进行抽取，得到fn1*lframe长度的矩阵d
        d = np.tile(np.arange(0, lframe), (fns, 1)) + np.tile(np.arange(0, fns * mframe, mframe), (lframe, 1)).T
        d = np.array(d, dtype=np.int32)  # 将d转换为矩阵形式（数据类型为int类型）
        signal = fillsignal[d]  # 得到帧信号
        return (signal, fns)
    
    # 函数调用，把采样数组、帧长、帧移等参数传递进函数frame，并返回帧信号endframe、帧数fns
    endframe, fns = frame(wave_data, lframe, mframe)
    
    # 加窗
    HanningWindow = scipy.signal.windows.hann(lframe)  # 调用汉明窗,把参数帧长传进去
    win = np.tile(HanningWindow, (fns, 1))  # 将窗函数HanningWindow平铺成fns行1列的数组
    # 利用mat把数组解释为矩阵，然后再进行汉明窗与信号数组的相乘
    everyframeaddwindow = np.multiply(np.mat(endframe), np.mat(win))
    # 将进行加窗后的每一帧变成数组形式
    signal1 = np.array(everyframeaddwindow)
    print('查看加窗后数组的形状：\n', np.shape(signal1))  # 显示加窗以后的矩阵形状：帧数*帧长(236, 200)
    
    signal2 = signal1.T  # 转置使函数signal2矩阵是（200, 236）
    print('查看转置后数组的形状：\n', np.shape(signal2))
    framenums = signal2.shape[1]  # 返回矩阵列数，获取帧数
    print(framenums)  # (200)
    
    # 寻找过零率
    # 初始化，补零,返回一个帧数fns行1列的零数组
    zeropoint = np.zeros((framenums, 1))
    for i in range (framenums):  # for循环从0到帧数framenums
        # 消除直流分量，提高短时过零率的正确估算
        signal2[:, i] = signal2[:, i] - np.mean(signal2[:, i])
        # 利用公式寻找过零点
        zeropoint[i] = sum(signal2[0:-1, i] * signal2[1:, i] < 0)
    print(zeropoint)
    
    # 语音采样点对应时间
    t1 = np.arange(0, len(wave_data)) * (1.0 / framerate)
    # 过零点对应时间刻度
    t2 = np.arange(0, len(zeropoint)) * (len(wave_data)/len(zeropoint) / framerate)
    print('过零点：\n', t2)
    # 画出语音波形
    plt.plot(t1, wave_data)
    plt.xlabel('times/s')  # x轴时间
    plt.ylabel('Amplitude-幅值')  # y轴振幅
    plt.title('bluesky3.wav—语音波形', fontsize=12, color='black')  # 标题名称、字体大小、颜色
    plt.show()
    # 画出语音平均过零率图形
    plt.plot(t2, zeropoint/lframe)
    plt.xlabel('times/s')  # x轴时间
    plt.ylabel('过零率')  # y轴平均过零率
    plt.title('bluesky3.wav—短时平均过零率', fontsize=12, color='black')  # 标题名称、字体大小、颜色
    plt.show()

输出语音信号的波形：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80NTI3MjkwOA_size_16_color_FFFFFF_t_70]  
输出短时过零率：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80NTI3MjkwOA_size_16_color_FFFFFF_t_70 1]

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80NTI3MjkwOA_size_16_color_FFFFFF_t_70]: /images/20221005/9e0bdb090b6743f0914d2d83b9a2ceae.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80NTI3MjkwOA_size_16_color_FFFFFF_t_70 1]: /images/20221005/1f9b649e437c4a868153e8722d595ee4.png