python做语音信号处理

时间：12-08来源：作者：点击数：41

音频信号的读写、播放及录音

python已经支持WAV格式的书写，而实时的声音输入输出需要安装pyAudio(http://people.csail.mit.edu/hubert/pyaudio)。最后我们还将使用pyMedia(http://pymedia.org)进行Mp3的解码和播放。

读取音频文件

librosa库 (推荐)

这是我最常用也是最喜欢的语音库，librosa是python第三方库，我们在使用前需要在cmd终端运行： pip install librosa 关于librosa的介绍请看文章librosa语音信号处理。

import librosa

y, sr = librosa.load(path, sr=fs)

该函数是会改变声音的采样频率的。如果 sr 缺省，librosa.load()会默认以22050的采样率读取音频文件，高于该采样率的音频文件会被下采样，低于该采样率的文件会被上采样。因此，如果希望以原始采样率读取音频文件，sr 应当设为 None。具体做法为 y, sr = librosa(filename, sr=None)。

音频数据 y 是直接经过归一化的数组

soundfile库(推荐)

soundfile库也是我常用的读取语音的库，有时候他的读取速度会比librosa更快，他只能读取原始音频，并不会做重采样。

import soundfile as sf 

wav, wav_sr = sf.read(wav_path, always_2d=True, dtype='float32')

wave库

wave库是python的标准库，对于python来说相对底层，wave不支持压缩/解压，但支持单声道/立体声语音的读取。

wave_read = wave.open(file,mode="rb")

参数：

f：语音文件名或文件路径
mode：读或写
- "rb"：只读模式
- "wb"：只写模式

返回：读取的文件流

该open()函数可用于with声明中。当with块完成时，wave_read.close()或wave_write.close()方法被调用

文件路径：

例如voice.wav文件在路径C:\Users\Never\Desktop\code for the speech的文件夹里

则file有以下三种填写格式：

r"C:\Users\Never\Desktop\code for the speech\voice.wav"

"C:/Users/Never/Desktop/code for the speech/voice.wav"

"C:\\Users\\Never\\Desktop\\code for the speech\\voice.wav"

三者等价，右划线\为转意字符，如果要表达\则需要\\，引号前面加r表示原始字符串。

wave_read.getparams()：一次性返回所有的音频参数，返回的是一个元组(声道数，量化位数(byte单位)，采样频率，采样点数，压缩类型，压缩类型的描述)。(nchannels, sampwidth, framerate, nframes, comptype, compname)wave模块只支持非压缩的数据，因此可以忽略最后两个信息。

str_data = wave_read.readframes(nframes)：读取的长度(以取样点为单位)，返回的是字符串类型的数据

wave_data = np.fromstring(str_data, dtype=np.float16)：将上面字符串类型数据转换为一维float16类型的数组。

现在的wave_data是一个一维的short类型的数组，但是因为我们的声音文件是双声道的，因此它由左右两个声道的取样交替构成：LR

wave_data.shape = (-1, 2)# -1的意思就是没有指定,根据另一个维度的数量进行分割，得到n行2列的数组。

wave_read.close（）关闭文件流wave
wave_read.getnchannels（）返回音频通道的数量（1对于单声道，2对于立体声）。
wave_read.getsampwidth（）以字节为单位返回样本宽度
wave_read.getframerate（）返回采样频率。
wave_read.getnframes（）　返回音频帧数。
wave_read.rewind（）将文件指针倒回到音频流的开头。
wave_read.tell（）返回当前文件指针位置。

# -*- coding: utf-8 -*-
# 读Wave文件并且绘制波形
import wave
import matplotlib.pyplot as plt
import numpy as np

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示符号

# 打开WAV音频
f = wave.open(r"C:\Windows\media\Windows Background.wav", "rb")

# 读取格式信息
# (声道数、量化位数、采样频率、采样点数、压缩类型、压缩类型的描述)
# (nchannels, sampwidth, framerate, nframes, comptype, compname)
params = f.getparams()
nchannels, sampwidth, framerate, nframes = params[:4]
# nchannels通道数 = 2
# sampwidth量化位数 = 2
# framerate采样频率 = 22050
# nframes采样点数 = 53395

# 读取nframes个数据，返回字符串格式
str_data = f.readframes(nframes)

f.close()

# 将字符串转换为数组，得到一维的short类型的数组
wave_data = np.fromstring(str_data, dtype=np.short)

# 赋值的归一化
wave_data = wave_data * 1.0 / (max(abs(wave_data)))

# 整合左声道和右声道的数据
wave_data = np.reshape(wave_data, [nframes, nchannels])
# wave_data.shape = (-1, 2)   # -1的意思就是没有指定,根据另一个维度的数量进行分割

# 最后通过采样点数和取样频率计算出每个取样的时间
time = np.arange(0, nframes) * (1.0 / framerate)

plt.figure()
# 左声道波形
plt.subplot(2, 1, 1)
plt.plot(time, wave_data[:, 0])
plt.xlabel("时间/s",fontsize=14)
plt.ylabel("幅度",fontsize=14)
plt.title("左声道",fontsize=14)
plt.grid()  # 标尺

plt.subplot(2, 1, 2)
# 右声道波形
plt.plot(time, wave_data[:, 1], c="g")
plt.xlabel("时间/s",fontsize=14)
plt.ylabel("幅度",fontsize=14)
plt.title("右声道",fontsize=14)

plt.tight_layout()  # 紧密布局
plt.show()

读取通道数为2的音频信号

scipy库

from scipy.io import wavfile

sampling_freq, audio = wavfile.read("***.wav")

audio 是直接经过归一化的数组

写音频文件

soundfile库 (推荐)

在0.8.0以后的版本，librosa都会将这个函数删除，推荐用下面的函数：

import soundfile as sf

sf.write(file, data, samplerate)

参数：

file：保存输出wav文件的路径
data：音频数据
samplerate：采样率

wave库

在写入第一帧数据时，先通过调用setnframes()设置好帧数，setnchannels()设置好声道数，setsampwidth()设置量化位数，setframerate()设置好采样频率，然后writeframes(wave.tostring())用于写入帧数据。

wave_write = wave.open(file,mode="wb")

wave_write是写文件流，

wave_write.setnchannels（n）设置通道数。
wave_write.setsampwidth（n）将样本宽度设置为n个字节，量化位数
wave_write.setframerate（n）将采样频率设置为n。
wave_write.setnframes（n）将帧数设置为n
wave_write.setparams（tuple）以元组形式设置所有参数(nchannels, sampwidth, framerate, nframes,comptype, compname)
wave_write.writeframes（data）写入data个长度的音频，以采样点为单位
wave_write.tell（）返回文件中的当前位置

# Author:凌逆战
# -*- coding:utf-8 -*-
import wave
import numpy as np
import scipy.signal as signal

framerate = 44100   # 采样频率
time = 10           # 持续时间

t = np.arange(0, time, 1.0/framerate)

# 调用scipy.signal库中的chrip函数，
# 产生长度为10秒、取样频率为44.1kHz、100Hz到1kHz的频率扫描波
wave_data = signal.chirp(t, 100, time, 1000, method='linear') * 10000

# 由于chrip函数返回的数组为float64型，
# 需要调用数组的astype方法将其转换为short型。
wave_data = wave_data.astype(np.short)

# 打开WAV音频用来写操作
f = wave.open(r"sweep.wav", "wb")

f.setnchannels(1)           # 配置声道数
f.setsampwidth(2)           # 配置量化位数
f.setframerate(framerate)   # 配置取样频率
comptype = "NONE"
compname = "not compressed"

# 也可以用setparams一次性配置所有参数
# outwave.setparams((1, 2, framerate, nframes,comptype, compname))

# 将wav_data转换为二进制数据写入文件
f.writeframes(wave_data.tostring())
f.close()

写wav文件

# Author:凌逆战
# -*- coding:utf-8 -*-
import wave
import numpy as np
import struct

f = wave.open(r"C:\Windows\media\Windows Background.wav", "rb")
params = f.getparams()
nchannels, sampwidth, framerate, nframes = params[:4]
strData = f.readframes(nframes)
waveData = np.fromstring(strData,dtype=np.int16)
f.close()
waveData = waveData*1.0/(max(abs(waveData)))

# wav文件写入
# 待写入wav的数据，这里仍然取waveData数据
outData = waveData
outwave = wave.open("write.wav", 'wb')
nchannels = 1   # 通道数设置为1
sampwidth = 2   # 量化位数设置为2
framerate = 8000    # 采样频率8000
nframes = len(outData)    # 采样点数

comptype = "NONE"
compname = "not compressed"
outwave.setparams((nchannels, sampwidth, framerate, nframes,
    comptype, compname))

for i in outData:
        outwave.writeframes(struct.pack('h', int(i * 64000 / 2)))

        # struct.pack(FMT, V1)将V1的值转换为FMT格式字符串
outwave.close()

写WAV文件方法2

scipy库

from scipy.io.wavfile import write

write(output_filename, freq, audio)

import numpy as np
import matplotlib.pyplot as plt
from scipy.io.wavfile import write

# 定义存储音频的输出文件
output_file = 'output_generated.wav'

# 指定音频生成的参数
duration = 3            # 单位秒
sampling_freq = 44100   # 单位Hz
tone_freq = 587         # 音调的频率
min_val = -2 * np.pi
max_val = 2 * np.pi

# 生成音频信号
t = np.linspace(min_val, max_val, duration * sampling_freq)
audio = np.sin(2 * np.pi * tone_freq * t)

# 添加噪声(duration * sampling_freq个(0,1]之间的随机值)
noise = 0.4 * np.random.rand(duration * sampling_freq)
audio += noise

scaling_factor = pow(2,15) - 1  # 转换为16位整型数
audio_normalized = audio / np.max(np.abs(audio))    # 归一化
audio_scaled = np.int16(audio_normalized * scaling_factor)  # 这句话什么意思

write(output_file, sampling_freq, audio_scaled) # 写入输出文件

audio = audio[:300] # 取前300个音频信号

x_values = np.arange(0, len(audio), 1) / float(sampling_freq)
x_values *= 1000    # 将时间轴单位转换为秒

plt.plot(x_values, audio, color='blue')
plt.xlabel('Time (ms)')
plt.ylabel('Amplitude')
plt.title('Audio signal')
plt.show()

写WAV文件

合成有音调的音乐

import json
import numpy as np
from scipy.io.wavfile import write
import matplotlib.pyplot as plt

# 定义合成音调
def Synthetic_tone(freq, duration, amp=1.0, sampling_freq=44100):
    # 建立时间轴
    t = np.linspace(0, duration, duration * sampling_freq)
    # 构建音频信号
    audio = amp * np.sin(2 * np.pi * freq * t)
    return audio.astype(np.int16)


# json文件中包含一些音阶以及他们的频率
tone_map_file = 'tone_freq_map.json'

# 读取频率映射文件
with open(tone_map_file, 'r') as f:
    tone_freq_map = json.loads(f.read())
    print(tone_freq_map)
# {'A': 440, 'Asharp': 466, 'B': 494, 'C': 523, 'Csharp': 554, 'D': 587, 'Dsharp': 622, 'E': 659, 'F': 698, 'Fsharp': 740, 'G': 784, 'Gsharp': 831}

# 设置生成G调的输入参数
input_tone = 'G'
duration = 2             # seconds
amplitude = 10000        # 振幅
sampling_freq = 44100    # Hz
# 生成音阶
synthesized_tone = Synthetic_tone(tone_freq_map[input_tone], duration, amplitude, sampling_freq)

# 写入输出文件
write('output_tone.wav', sampling_freq, synthesized_tone)

# 音阶及其连续时间
tone_seq = [('D', 0.3), ('G', 0.6), ('C', 0.5), ('A', 0.3), ('Asharp', 0.7)]

# 构建基于和弦序列的音频信号
output = np.array([])
for item in tone_seq:
    input_tone = item[0]
    duration = item[1]
    synthesized_tone = Synthetic_tone(tone_freq_map[input_tone], duration, amplitude, sampling_freq)
    output = np.append(output, synthesized_tone, axis=0)

# 写入输出文件
write('output_tone_seq.wav', sampling_freq, output)

合成音调

{
    "A": 440,
    "Asharp": 466,
    "B": 494,
    "C": 523,
    "Csharp": 554,
    "D": 587,
    "Dsharp": 622,
    "E": 659,
    "F": 698,
    "Fsharp": 740,
    "G": 784,
    "Gsharp": 831
}

tone_freq_map

音频播放

wav音频播放用到的是pyaudio库

p = pyaudio.PyAudio()
stream = p.open(format = p.get_format_from_width(sampwidth) , channels ,rate ,output = True)
stream.write(data)# 播放data数据

以下列出pyaudio对象的open()方法的主要参数：

rate：取样频率
channels：声道数
format：取样值的量化格式 (paFloat32, paInt32, paInt24, paInt16, paInt8 ...)。在上面的例子中，使用get_format_from_width方法将wf.sampwidth()的返回值2转换为paInt16
input：输入流标志，如果为True的话则开启输入流
output：输出流标志，如果为True的话则开启输出流
input_device_index：输入流所使用的设备的编号，如果不指定的话，则使用系统的缺省设备
output_device_index：输出流所使用的设备的编号，如果不指定的话，则使用系统的缺省设备
frames_per_buffer：底层的缓存的块的大小，底层的缓存由N个同样大小的块组成
start：指定是否立即开启输入输出流，缺省值为True

# -*- coding: utf-8 -*-
import pyaudio
import wave

chunk = 1024

wf = wave.open(r"c:\WINDOWS\Media\Windows Background.wav", 'rb')

p = pyaudio.PyAudio()

# 打开声音输出流
stream = p.open(format = p.get_format_from_width(wf.getsampwidth()),
                channels = wf.getnchannels(),
                rate = wf.getframerate(),
                output = True)

# 写声音输出流到声卡进行播放
while True:
    data = wf.readframes(chunk)
    if data == "":
        break
    stream.write(data)

stream.stop_stream()
stream.close()
p.terminate()   # 关闭PyAudio

播放wav音频

录音

以SAMPLING_RATE为采样频率，每次读入一块有NUM_SAMPLES个采样的数据块，当读入的采样数据中有COUNT_NUM个值大于LEVEL的取样的时候，将数据保存进WAV文件，一旦开始保存数据，所保存的数据长度最短为SAVE_LENGTH个块。WAV文件以保存时的时刻作为文件名。

从声卡读入的数据和从WAV文件读入的类似，都是二进制数据，由于我们用paInt16格式(16bit的short类型)保存采样值，因此将它自己转换为dtype为np.short的数组。

'''
以SAMPLING_RATE为采样频率，
每次读入一块有NUM_SAMPLES个采样点的数据块，
当读入的采样数据中有COUNT_NUM个值大于LEVEL的取样的时候，
将采样数据保存进WAV文件，
一旦开始保存数据，所保存的数据长度最短为SAVE_LENGTH个数据块。

从声卡读入的数据和从WAV文件读入的类似，都是二进制数据，
由于我们用paInt16格式(16bit的short类型)保存采样值，
因此将它自己转换为dtype为np.short的数组。
'''


from pyaudio import PyAudio, paInt16
import numpy as np
import wave

# 将data中的数据保存到名为filename的WAV文件中
def save_wave_file(filename, data):
    wf = wave.open(filename, 'wb')
    wf.setnchannels(1)          # 单通道
    wf.setsampwidth(2)          # 量化位数
    wf.setframerate(SAMPLING_RATE)  # 设置采样频率
    wf.writeframes(b"".join(data))  # 写入语音帧
    wf.close()


NUM_SAMPLES = 2000      # pyAudio内部缓存块的大小
SAMPLING_RATE = 8000    # 取样频率
LEVEL = 1500           # 声音保存的阈值，小于这个阈值不录
COUNT_NUM = 20 # 缓存快类如果有20个大于阈值的取样则记录声音
SAVE_LENGTH = 8 # 声音记录的最小长度：SAVE_LENGTH * NUM_SAMPLES 个取样

# 开启声音输入
pa = PyAudio()
stream = pa.open(format=paInt16, channels=1, rate=SAMPLING_RATE, input=True,
                frames_per_buffer=NUM_SAMPLES)

save_count = 0  # 用来计数
save_buffer = []    #

while True:
    # 读入NUM_SAMPLES个取样
    string_audio_data = stream.read(NUM_SAMPLES)
    # 将读入的数据转换为数组
    audio_data = np.fromstring(string_audio_data, dtype=np.short)
    # 计算大于LEVEL的取样的个数
    large_sample_count = np.sum( audio_data > LEVEL )
    print(np.max(audio_data))
    # 如果个数大于COUNT_NUM，则至少保存SAVE_LENGTH个块
    if large_sample_count > COUNT_NUM:
        save_count = SAVE_LENGTH
    else:
        save_count -= 1

    if save_count < 0:
        save_count = 0

    if save_count > 0:
        # 将要保存的数据存放到save_buffer中
        save_buffer.append( string_audio_data )
    else:
        # 将save_buffer中的数据写入WAV文件，WAV文件的文件名是保存的时刻
        if len(save_buffer) > 0:
            filename = "recorde" + ".wav"
            save_wave_file(filename, save_buffer)
            print(filename, "saved")
            break

语音信号处理

语音信号是一个非平稳的时变信号，但语音信号是由声门的激励脉冲通过声道形成的，而声道(人的口腔、鼻腔)的肌肉运动是缓慢的，所以“短时间内可以认为语音信号是平稳时不变的，一般 10~30ms。由此构成了语音信号的“短时分析技术”。在短时分析中，将语音信号分为一段一段的语音帧，每一帧一般取10~30ms，我们的研究就建立在每一帧的语音特征分析上。

提取的不同的语音特征参数对应着不同的语音信号分析方法：时域分析、频域分析、倒谱域分析...

分帧

加窗的代价是一帧信号两端的部分被削弱了，没有像中央的部分那样得到重视。弥补的办法就是帧重叠。相邻两帧的起始位置的时间差叫做帧移(或者理解为后一帧第前一帧的偏移量)，常见的取法是取为帧长的一半。帧长(wlen) = 重叠(overlap)+帧移(inc)。fn表示一段语音信号的分帧数。

$frame_num=\frac{N-overlap}{inc}=\frac{N-wlen+inc}{inc}$

# librosa有专门的分帧函数
librosa.util.frame(x,frame_length=512,hop_length=256)

加窗

对连续的语音分帧做STFT处理，等价于截取一段时间信号，对其进行周期性延拓，从而变成无限长序列，并对该无限长序列做FFT变换，这一截断并不符合傅里叶变换的定义。因此，会导致频谱泄漏和混叠

频谱泄漏：如果不加窗，默认就是矩形窗，时域的乘积就是频域的卷积，使得频谱以实际频率值为中心，以窗函数频谱波形的形状向两侧扩散，指某一频点能量扩散到相邻频点的现象，会导致幅度较小的频点淹没在幅度较大的频点泄漏分量中
频谱混叠：会在分段拼接处引入虚假的峰值，进而不能获得准确的频谱情况

加窗的目的是：让一帧信号的幅度在两端渐变到 0，渐变对傅里叶变换有好处，可以让频谱上的各个峰更细，不容易糊在一起，从而减轻频谱泄漏和混叠的影响

加窗的代价是：一帧信号两端的部分被削弱了，没有像中央的部分那样得到重视。弥补的办法就是相互重叠。相邻两帧的起始位置的时间差叫做帧移，常见的取法是取为帧长的一半。

对于语音，窗函数常选汉宁窗(Hanning)、汉明窗(Hamming)、凯撒窗(Kaiser)及其改进窗，他们的时域波形和幅频响应如下所示：

1、汉宁窗(Hann)

$w(n) = 0.5 - 0.5 \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1$

2、汉明窗(Hamming)

$w(n) = 0.54 - 0.46 \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1$

# -*- coding:utf-8 -*-
# Author:凌逆战 | Never
# Date: 2023/1/1
"""
绘制 窗函数和对应的频率响应
"""
import numpy as np
from numpy.fft import rfft
import matplotlib.pyplot as plt

window_len = 60


# frequency response
def frequency_response(window, window_len=window_len, NFFT=2048):
    A = rfft(window, NFFT) / (window_len / 2)  # (513,)
    mag = np.abs(A)
    freq = np.linspace(0, 0.5, len(A))
    # 忽略警告
    with np.errstate(divide='ignore', invalid='ignore'):
        response = 20 * np.log10(mag)
    response = np.clip(response, -150, 150)
    return freq, response


def Rectangle_windows(win_length):
    # 矩形窗
    return np.ones((win_length))


def Voibis_windows(win_length):
    """ Voibis_windows窗函数，RNNoise使用的是它，它满足Princen-Bradley准则。
    :param x:
    :param win_length: 窗长
    :return:
    """
    x = np.arange(0, win_length)
    return np.sin((np.pi / 2) * np.sin((np.pi * x) / win_length) ** 2)


def sqrt_hanning_windows(win_length, mode="periodic"):
    # symmetric: 对称窗，主要用于滤波器的设计
    # periodic: 周期窗，常用于频谱分析
    if mode == "symmetric":
        haning_window = np.hanning(win_length)
        sqrt_haning_window = np.sqrt(haning_window)
    elif mode == "periodic":
        haning_window = np.hanning(win_length+1)
        sqrt_haning_window = np.sqrt(haning_window)
        sqrt_haning_window = sqrt_haning_window[0:-1].astype('float32')
    return sqrt_haning_window


Rectangle_windows = Rectangle_windows(window_len)
hanning_window = np.hanning(M=window_len)
print(np.argmax(hanning_window))
sqrt_hanning_windows = sqrt_hanning_windows(window_len)
hamming_window = np.hamming(M=window_len)
Voibis_windows = Voibis_windows(window_len)
blackman_window = np.blackman(M=window_len)
bartlett_window = np.bartlett(M=window_len)
kaiser_window = np.kaiser(M=window_len, beta=14)

plt.figure()
plt.plot(Rectangle_windows, label="Rectangle")
plt.plot(hanning_window, label="hanning")
plt.plot(sqrt_hanning_windows, label="sqrt_hanning")
plt.plot(hamming_window, label="hamming")
plt.plot(Voibis_windows, label="Voibis")
plt.plot(blackman_window, label="blackman")
plt.plot(bartlett_window, label="bartlett")
plt.plot(kaiser_window, label="kaiser")

plt.legend()
plt.tight_layout()
plt.show()

freq, Rectangle_FreqResp = frequency_response(Rectangle_windows, window_len)
freq, hanning_FreqResp = frequency_response(hanning_window, window_len)
freq, sqrt_hanning_FreqResp = frequency_response(sqrt_hanning_windows, window_len)
freq, hamming_FreqResp = frequency_response(hamming_window, window_len)
freq, Voibis_FreqResp = frequency_response(Voibis_windows, window_len)
freq, blackman_FreqResp = frequency_response(blackman_window, window_len)
freq, bartlett_FreqResp = frequency_response(bartlett_window, window_len)
freq, kaiser_FreqRespw = frequency_response(kaiser_window, window_len)

plt.figure()
plt.title("Frequency response")
plt.plot(freq, Rectangle_FreqResp, label="Rectangle")
plt.plot(freq, hanning_FreqResp, label="hanning")
plt.plot(freq, sqrt_hanning_FreqResp, label="sqrt_hanning")
plt.plot(freq, hamming_FreqResp, label="hamming")
plt.plot(freq, Voibis_FreqResp, label="Voibis")
plt.plot(freq, blackman_FreqResp, label="blackman")
plt.plot(freq, bartlett_FreqResp, label="bartlett")
plt.plot(freq, kaiser_FreqRespw, label="kaiser")
plt.ylabel("Magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.legend()
plt.tight_layout()
plt.show()

绘制窗函数和对应的频率响应

想要更进一步了解窗函数可以移步文章：语音信号处理中的“窗函数”

overlap and add

将分帧好的语音拼接回完整的语音，当前帧的前半部分+下一帧的后半部分=1

def overlap_add_2(win_array):
    wav_sys = np.zeros(window_len + frame_len * (frame_num - 1))
    for frame_index in range(frame_num):
        ytmp = win_array[:, frame_index]
        wav_sys[frame_index * frame_len: (frame_index * frame_len + window_len)] += ytmp
    wav_sys = wav_sys[frame_len: -frame_len]
    print("wav_sys.shape", wav_sys.shape)
    return wav_sys

# -*- coding:utf-8 -*-
# Author:凌逆战 | Never
# Date: 2023/1/2
"""

"""
import librosa
import numpy as np
import matplotlib.pyplot as plt
import soundfile
from librosa.filters import get_window

sr = 16000
frame_len = 256
window_len = 512
NFFT = 512
fft_window = get_window("hann", Nx=window_len, fftbins=True)  # 用于频率分析
wav = librosa.load("./p225_001.wav", sr=sr)[0]
wav = wav[:len(wav) - len(wav) % frame_len]
print("wav.shape", wav.shape)

# 如果不补零的话，前半帧和后半帧 会因为加窗而无法恢复
wav_pad = np.pad(wav, (frame_len, frame_len), mode="constant")  # center=True
print("wav_pad.shape", wav_pad.shape)
# librosa有专门的分帧函数
frame_array = librosa.util.frame(wav_pad, frame_length=window_len, hop_length=frame_len)  # (帧长,帧数) (512, 129)
frame_num = frame_array.shape[1]
# 加窗、FFT
win_array = np.zeros_like(frame_array)
for frame_index in range(frame_num):
    win_array[:, frame_index] = frame_array[:, frame_index] * fft_window  # (512, 129)


# ifft、加窗、overlap_add
def overlap_add_1(win_array):
    sys_frame = []
    previous_frame = np.zeros((frame_len))
    for frame_index in range(frame_num):
        current_frame = win_array[:, frame_index]  # 当前窗 (512,)
        sys_frame.append(previous_frame + current_frame[:frame_len])
        previous_frame = current_frame[frame_len:]
        if frame_index == frame_num - 1:
            sys_frame.append(current_frame[frame_len:])

    wav_sys = np.concatenate(sys_frame, axis=0)
    wav_sys = wav_sys[frame_len: -frame_len]
    print("wav_sys.shape", wav_sys.shape)
    return wav_sys


def overlap_add_2(win_array):
    wav_sys = np.zeros(window_len + frame_len * (frame_num - 1))
    for frame_index in range(frame_num):
        ytmp = win_array[:, frame_index]
        wav_sys[frame_index * frame_len: (frame_index * frame_len + window_len)] += ytmp
    wav_sys = wav_sys[frame_len: -frame_len]
    print("wav_sys.shape", wav_sys.shape)
    return wav_sys


# https://github.com/miralv/Deep-Learning-for-Speech-Enhancement/blob/b2f3d4e33fdc8a1d75b774f009aadf95616efc99/recoverSignal.py
def overlap_add_3(win_array):
    wav_sys = np.zeros(window_len + frame_len * (frame_num - 1))
    wav_sys[0:frame_len] = win_array[0:frame_len, 0]
    start_point = frame_len
    for i in range(0, (frame_num - 1)):
        # Add the elements corresponding to the current half window

        wav_sys[start_point:start_point + frame_len] = np.add(win_array[frame_len:, i], win_array[0:frame_len, i + 1])
        start_point += frame_len

    # Add the last half window manually
    wav_sys[start_point:] = win_array[frame_len:, frame_num - 1]
    return wav_sys


wav_sys = overlap_add_1(win_array)
soundfile.write("./overlap_add.wav", data=wav_sys, samplerate=sr)

plt.subplot(2, 1, 1)
plt.plot(wav)
plt.subplot(2, 1, 2)
plt.plot(wav_sys)
plt.show()

语音信号的短时时域处理

短时能量和短时平均幅度

短时能量和短时平均幅度的主要用途：

区分浊音和清音段，因为浊音的短时能量 $E(i)$ 比清音大很多；
区分声母和韵母的分界和无话段和有话段的分界

短时平均过零率

对于连续语音信号，过零率意味着时域波形通过时间轴，对于离散信号，如果相邻的取样值改变符号，则称为过零。

作用：

发浊音时由于声门波引起谱的高频跌落，所以语音信号能量约集中在3kHz以下
发清音时多数能量集中在较高的频率上，

因为高频意味着高的短时平均过零率，低频意味着低的短时平均过零率，所以浊音时具有较低的过零率，而清音时具有较高的过零率。

利用短时平均过零率可以从背景噪声中找出语音信号，
可以用于判断寂静无话段与有话段的起点和终止位置。
在背景噪声较小的时候，用平均能量识别较为有效，在背景噪声较大的时候，用短时平均过零率识别较为有效。

短时自相关函数

短时自相关函数主要应用于端点检测和基音的提取，在韵母基因频率整数倍处将出现峰值特性，通常根据除R(0)外的第一峰值来估计基音，而在声母的短时自相关函数中看不到明显的峰值。

短时平均幅度差函数

用于检测基音周期，而且在计算上比短时自相关函数更加简单。

语音信号的短时频域处理

在语音信号处理中，在语音信号处理中，信号在频域或其他变换域上的分析处理占重要的位置，在频域上研究语音可以使信号在时域上无法表现出来的某些特征变得十分明显，一个音频信号的本质是由其频率内容决定的，

将时域信号转换为频域信号一般对语音进行短时傅里叶变换。

fft_audio = np.fft.fft(audio)

import numpy as np
from scipy.io import wavfile
import matplotlib.pyplot as plt

sampling_freq, audio = wavfile.read(r"C:\Windows\media\Windows Background.wav")   # 读取文件

audio = audio / np.max(audio)   # 归一化，标准化

# 应用傅里叶变换
fft_signal = np.fft.fft(audio)
print(fft_signal)
# [-0.04022912+0.j         -0.04068997-0.00052721j -0.03933007-0.00448355j
#  ... -0.03947908+0.00298096j -0.03933007+0.00448355j -0.04068997+0.00052721j]

fft_signal = abs(fft_signal)
print(fft_signal)
# [0.04022912 0.04069339 0.0395848  ... 0.08001755 0.09203427 0.12889393]

# 建立时间轴
Freq = np.arange(0, len(fft_signal))

# 绘制语音信号的
plt.figure()
plt.plot(Freq, fft_signal, color='blue')
plt.xlabel('Freq (in kHz)')
plt.ylabel('Amplitude')
plt.show()

绘制语音信号的频谱图

提取频域特征

将信号转换为频域之后，还需要将其转换为有用的形式，梅尔频率倒谱系数(MFCC)，MFCC首先计算信号的功率谱，然后用滤波器组和离散余弦变换的组合来提取特征。

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from python_speech_features import mfcc, logfbank

# 读取输入音频文件
sampling_freq, audio = wavfile.read("input_freq.wav")

# 提取MFCC和滤波器组特征
mfcc_features = mfcc(audio, sampling_freq)
filterbank_features = logfbank(audio, sampling_freq)

print('\nMFCC:\n窗口数 =', mfcc_features.shape[0])
print('每个特征的长度 =', mfcc_features.shape[1])
print('\nFilter bank:\n窗口数 =', filterbank_features.shape[0])
print('每个特征的长度 =', filterbank_features.shape[1])

# 画出特征图，将MFCC可视化。转置矩阵，使得时域是水平的
mfcc_features = mfcc_features.T
plt.matshow(mfcc_features)
plt.title('MFCC')
# 将滤波器组特征可视化。转置矩阵，使得时域是水平的
filterbank_features = filterbank_features.T
plt.matshow(filterbank_features)
plt.title('Filter bank')

plt.show()

提取MFCC特征

语谱图

绝大部分信号都可以分解为若干不同频率的正弦波。

这些正弦波中，频率最低的称为信号的基波，其余称为信号的谐波。

基波只有一个，可以称为一次谐波，谐波可以有很多个，每次谐波的频率是基波频率的整数倍。谐波的大小可能互不相同。

以谐波的频率为横坐标，幅值（大小）为纵坐标，绘制的系列条形图，称为频谱。频谱能够准确反映信号的内部构造。

语谱图综合了时域和频域的特点，明显的显示出来了语音频率随时间的变化情况，语谱图的横轴为时间，纵轴为频率任意给定频率成分在给定时刻的强弱用颜色深浅表示。颜色深表示频谱值大，颜色浅表示频谱值小，语谱图上不同的黑白程度形成不同的纹路，称为声纹，不用讲话者的声纹是不一样的，可以用做声纹识别。

其实得到了分帧信号，频域变换取幅值，就可以得到语谱图，如果仅仅是观察，matplotlib.pyplot有specgram指令：

import wave
import matplotlib.pyplot as plt
import numpy as np

f = wave.open(r"C:\Windows\media\Windows Background.wav", "rb")
params = f.getparams()
nchannels, sampwidth, framerate, nframes = params[:4]
strData = f.readframes(nframes)#读取音频，字符串格式
waveData = np.fromstring(strData,dtype=np.int16)#将字符串转化为int
waveData = waveData*1.0/(max(abs(waveData)))#wave幅值归一化
waveData = np.reshape(waveData,[nframes,nchannels]).T
f.close()

plt.specgram(waveData[0],Fs = framerate, scale_by_freq = True, sides = 'default')
plt.ylabel('Frequency(Hz)')
plt.xlabel('Time(s)')
plt.show()

语谱图

[Y,FS]=audioread('p225_355_wb.wav');

% specgram(Y,2048,44100,2048,1536);
    %Y1为波形数据
    %FFT帧长2048点(在44100Hz频率时约为46ms)
    %采样频率44.1KHz
    %加窗长度，一般与帧长相等
    %帧重叠长度，此处取为帧长的3/4
specgram(Y,2048,FS,2048,1536);
xlabel('时间(s)')
ylabel('频率(Hz)')
title('语谱图')

MATLAB语谱图

语音识别

import os
import numpy as np
import scipy.io.wavfile as wf
import python_speech_features as sf
import hmmlearn.hmm as hl


# 1. 读取training文件夹中的训练音频样本，每个音频对应一个mfcc矩阵，每个mfcc都有一个类别(apple...)
def search_file(directory):
    """
    :param directory: 训练音频的路径
    :return: 字典{'apple':[url, url, url ... ], 'banana':[...]}
    """
    # 使传过来的directory匹配当前操作系统
    directory = os.path.normpath(directory)
    objects = {}
    # curdir：当前目录
    # subdirs: 当前目录下的所有子目录
    # files: 当前目录下的所有文件名
    for curdir, subdirs, files in os.walk(directory):
        for file in files:
            if file.endswith('.wav'):
                label = curdir.split(os.path.sep)[-1]  # os.path.sep为路径分隔符
                if label not in objects:
                    objects[label] = []
                # 把路径添加到label对应的列表中
                path = os.path.join(curdir, file)
                objects[label].append(path)
    return objects


# 读取训练集数据
train_samples = search_file('../machine_learning_date/speeches/training')

"""
2. 把所有类别为apple的mfcc合并在一起，形成训练集。
    训练集:
    train_x：[mfcc1,mfcc2,mfcc3,...],[mfcc1,mfcc2,mfcc3,...]...
    train_y：[apple],[banana]...
由上述训练集样本可以训练一个用于匹配apple的HMM。"""

train_x, train_y = [], []
# 遍历字典
for label, filenames in train_samples.items():
    # [('apple', ['url1,,url2...'])
    # [("banana"),("url1,url2,url3...")]...
    mfccs = np.array([])
    for filename in filenames:
        sample_rate, sigs = wf.read(filename)
        mfcc = sf.mfcc(sigs, sample_rate)
        if len(mfccs) == 0:
            mfccs = mfcc
        else:
            mfccs = np.append(mfccs, mfcc, axis=0)
    train_x.append(mfccs)
    train_y.append(label)

# 3.训练模型，有7个句子，创建了7个模型
models = {}
for mfccs, label in zip(train_x, train_y):
    model = hl.GaussianHMM(n_components=4, covariance_type='diag', n_iter=1000)
    models[label] = model.fit(mfccs)  # # {'apple':object, 'banana':object ...}

"""
4. 读取testing文件夹中的测试样本，
    测试集数据：
        test_x  [mfcc1, mfcc2, mfcc3...]
        test_y  [apple, banana, lime]
"""
test_samples = search_file('../machine_learning_date/speeches/testing')

test_x, test_y = [], []
for label, filenames in test_samples.items():
    mfccs = np.array([])
    for filename in filenames:
        sample_rate, sigs = wf.read(filename)
        mfcc = sf.mfcc(sigs, sample_rate)
        if len(mfccs) == 0:
            mfccs = mfcc
        else:
            mfccs = np.append(mfccs, mfcc, axis=0)
    test_x.append(mfccs)
    test_y.append(label)


# 5.测试模型
#    1. 分别使用7个HMM模型，对测试样本计算score得分。
#    2. 取7个模型中得分最高的模型所属类别作为预测类别。
pred_test_y = []
for mfccs in test_x:
    # 判断mfccs与哪一个HMM模型更加匹配
    best_score, best_label = None, None
    # 遍历7个模型
    for label, model in models.items():
        score = model.score(mfccs)
        if (best_score is None) or (best_score < score):
            best_score = score
            best_label = label
    pred_test_y.append(best_label)

print(test_y)   # ['apple', 'banana', 'kiwi', 'lime', 'orange', 'peach', 'pineapple']
print(pred_test_y)  # ['apple', 'banana', 'kiwi', 'lime', 'orange', 'peach', 'pineapple']