python做語音信號處理

作者：凌逆戰

時間：2019年11月1日

博客園地址：https://www.cnblogs.com/LXP-Never/p/10078200.html

---

音頻信號的讀寫、播放及錄音

標准的python已經支持WAV格式的書寫，而實時的聲音輸入輸出需要安裝pyAudio(http://people.csail.mit.edu/hubert/pyaudio)。最后我們還將使用pyMedia(http://pymedia.org)進行Mp3的解碼和播放。

　　音頻信號是模擬信號，我們需要將其保存為數字信號，才能對語音進行算法操作，WAV是Microsoft開發的一種聲音文件格式，通常被用來保存未壓縮的聲音數據。

語音信號有三個重要的參數：聲道數、取樣頻率和量化位數。

• 聲道數：可以是單聲道或者是雙聲道
• 采樣頻率：一秒內對聲音信號的采集次數，44100Hz采樣頻率意味着每秒鍾信號被分解成44100份。換句話說，每隔$\frac{1}{44100}$秒就會存儲一次，如果采樣率高，那么媒體播放音頻時會感覺信號是連續的。
• 量化位數：用多少bit表達一次采樣所采集的數據，通常有8bit、16bit、24bit和32bit等幾種

例如CD中所儲存的聲音信號是雙聲道、44.1kHz、16bit。

如果你需要自己錄制和編輯聲音文件，推薦使用Audacity(http://audacity.sourceforge.net)，它是一款開源的、跨平台、多聲道的錄音編輯軟件。在我的工作中經常使用Audacity進行聲音信號的錄制，然后再輸出成WAV文件供Python程序處理。

wave-讀wav文件

wava模塊為WAV聲音格式提供了方面的界面，他不支持壓縮/解壓，但支持單聲道/立體聲。

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

file通常為是字符串格式的文件名或者文件路徑

　　　　例如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表示原始字符串。

mode是缺省參數，可以不填，也可以是"rb":只讀模式；"wb":只寫模式。注意不支持讀/寫格式。

該open()函數可用於with聲明中。當with塊完成時，Wave_read.close()或Wave_write.close()方法被調用。Wave_read是讀取的文件流。

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.short)

將讀取的字符串數據轉換為一維short類型的數組。

通過np.fromstring函數將字符串轉換為數組，通過其參數dtype指定轉換后的數據格式(由於我們的聲音格式是以兩個字節表示一個取樣值，因此采用short數據類型轉換)

現在的wave_data是一個一維的short類型的數組，但是因為我們的聲音文件是雙聲道的，因此它由左右兩個聲道的取樣交替構成：LR

wave_data.shape = (-1, 2)　　# -1的意思就是沒有指定,根據另一個維度的數量進行分割，得到n行2列的數組。

getnchannels, getsampwidth, getframerate, getnframes等方法可以單獨返回WAV文件的特定的信息。

Wave_read.close（）　　關閉文件流wave

Wave_read.getnchannels（）　　返回音頻通道的數量（1對於單聲道，2對於立體聲）。

Wave_read.getsampwidth（）　　以字節為單位返回樣本寬度

Wave_read.getframerate（）　　返回采樣頻率。

Wave_read.getnframes（）　　　返回音頻幀數。

Wave_read.rewind（）　　　　　　將文件指針倒回到音頻流的開頭。

Wave_read.tell（）　　　　　　返回當前文件指針位置。 

讀取通道數為2的音頻信號

# -*- coding: utf-8 -*-
# 讀Wave文件並且繪制波形
import wave
import matplotlib.pyplot as plt
import numpy as np

# 打開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(3,1,1)
plt.plot(time, wave_data[:,0])
plt.xlabel("time (seconds)")
plt.ylabel("Amplitude")
plt.title("Left channel")
plt.grid()  # 標尺

plt.subplot(3,1,3)
# 右聲道波形
plt.plot(time, wave_data[:,1], c="g")
plt.xlabel("time (seconds)")
plt.ylabel("Amplitude")
plt.title("Left channel")
plt.title("right channel")
plt.grid()

plt.show()

效果圖：

第二種讀取文件的方式：

from scipy.io import wavfile

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

audio 是直接經過歸一化的數組

第三種讀取音頻的方式：

import librosa

y, sr = librosa.load(filename)

該函數是會改變聲音的采樣頻率的。如果 sr 缺省，librosa.load()會默認以22050的采樣率讀取音頻文件，高於該采樣率的音頻文件會被下采樣，低於該采樣率的文件會被上采樣。因此，如果希望以原始采樣率讀取音頻文件，sr 應當設為 None。具體做法為 y, sr = librosa(filename, sr=None)。

audio 是直接經過歸一化的數組

wave-寫wav音頻

在寫入第一幀數據時，先通過調用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（）　　返回文件中的當前位置

# -*- 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()

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()

第二種寫音頻文件的方法

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()

第三種寫音頻文件的方法

librosa.output.write_wav(path, y, sr, norm=False)

參數：

• path：str，保存輸出wav文件的路徑
• y：np.ndarry 音頻時間序列
• sr：y的采樣率
• norm：True/False，是否啟動幅值歸一化

合成有音調的音樂

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
}

音頻播放

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

錄音

以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，我們的研究就建立在每一幀的語音特征分析上。

　　提取的不同的語音特征參數對應着不同的語音信號分析方法：時域分析、頻域分析、倒譜域分析...由於語音信號最重要的感知特性反映在功率譜上，而相位變化只起到很小的作用，所有語音頻域分析更加重要。

信號加窗 

1、矩形窗

$$w(n)=\left\{\begin{matrix} 1&&0\leq n\leq L-1\\ 0&&其他 \end{matrix}\right.$$

2、漢明窗(Hamming)

$$w(n)=\left\{\begin{matrix} \frac{1}{2}(1-cos(\frac{2\pi n}{L-1}))&&0\leq n\leq L-1\\  0&&其他 \end{matrix}\right.$$

3、海寧窗(Hanning)

$$w(n)=\left\{\begin{matrix} 0.54-0.46cos(\frac{2\pi n}{L-1})&&0\leq n\leq L-1\\  0&&其他 \end{matrix}\right.$$

通常對信號截斷、分幀需要加窗，因為截斷都有頻域能量泄露，而窗函數可以減少截斷帶來的影響。

窗函數在scipy.signal信號處理工具箱中，如hanning窗：

import matplotlib.pyplot as plt
import scipy.signal as signal
plt.figure(figsize=(6,2))
plt.plot(signal.hanning(512))
plt.show()

信號分幀

　　在分幀中，相鄰兩幀之間會有一部分重疊，幀長(wlen) = 重疊(overlap)+幀移(inc)，如果相鄰兩幀之間不重疊，那么由於窗函數的形狀，截取到的語音幀邊緣會出現損失，所以要設置重疊部分。inc為幀移，表示后一幀第前一幀的偏移量，fs表示采樣率，fn表示一段語音信號的分幀數。

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

　　信號分幀的理論依據，其中x是語音信號，w是窗函數：

$$y(n)=\sum_{n=-N / 2+1}^{N / 2} x(m) w(n-m)$$

　　加窗截斷類似采樣，為了保證相鄰幀不至於差別過大，通常幀與幀之間有幀移，其實就是插值平滑的作用。

給出示意圖：

 這里主要用到numpy工具包，涉及的指令有：

• np.repeat：主要是直接重復

• np.tile：主要是周期性重復

對比一下：

向量情況：

矩陣情況：

對於數據：

repeat操作：

tile操作：

對應結果：

對應分幀的代碼實現：

 這是沒有加窗的示例：

import numpy as np
import wave
import os
#import math
 
def enframe(signal, nw, inc):
    '''將音頻信號轉化為幀。
    參數含義：
    signal:原始音頻型號
    nw:每一幀的長度(這里指采樣點的長度，即采樣頻率乘以時間間隔)
    inc:相鄰幀的間隔（同上定義）
    '''
    signal_length=len(signal) #信號總長度
    if signal_length<=nw: #若信號長度小於一個幀的長度，則幀數定義為1
        nf=1
    else: #否則，計算幀的總長度
        nf=int(np.ceil((1.0*signal_length-nw+inc)/inc))
    pad_length=int((nf-1)*inc+nw) #所有幀加起來總的鋪平后的長度
    zeros=np.zeros((pad_length-signal_length,)) #不夠的長度使用0填補，類似於FFT中的擴充數組操作
    pad_signal=np.concatenate((signal,zeros)) #填補后的信號記為pad_signal
    indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T  #相當於對所有幀的時間點進行抽取，得到nf*nw長度的矩陣
    indices=np.array(indices,dtype=np.int32) #將indices轉化為矩陣
    frames=pad_signal[indices] #得到幀信號
#    win=np.tile(winfunc(nw),(nf,1))  #window窗函數，這里默認取1
#    return frames*win   #返回幀信號矩陣
    return frames
def wavread(filename):
    f = wave.open(filename,'rb')
    params = f.getparams()
    nchannels, sampwidth, framerate, nframes = params[:4]
    strData = f.readframes(nframes)#讀取音頻，字符串格式
    waveData = np.fromstring(strData,dtype=np.int16)#將字符串轉化為int
    f.close()
    waveData = waveData*1.0/(max(abs(waveData)))#wave幅值歸一化
    waveData = np.reshape(waveData,[nframes,nchannels]).T
    return waveData
 
filepath = "./data/" #添加路徑
dirname= os.listdir(filepath) #得到文件夾下的所有文件名稱 
filename = filepath+dirname[0]
data = wavread(filename)
nw = 512
inc = 128
Frame = enframe(data[0], nw, inc) 

def enframe(signal, nw, inc, winfunc):
    '''將音頻信號轉化為幀。
    參數含義：
    signal:原始音頻型號
    nw:每一幀的長度(這里指采樣點的長度，即采樣頻率乘以時間間隔)
    inc:相鄰幀的間隔（同上定義）
    '''
    signal_length=len(signal) #信號總長度
    if signal_length<=nw: #若信號長度小於一個幀的長度，則幀數定義為1
        nf=1
    else: #否則，計算幀的總長度
        nf=int(np.ceil((1.0*signal_length-nw+inc)/inc))
    pad_length=int((nf-1)*inc+nw) #所有幀加起來總的鋪平后的長度
    zeros=np.zeros((pad_length-signal_length,)) #不夠的長度使用0填補，類似於FFT中的擴充數組操作
    pad_signal=np.concatenate((signal,zeros)) #填補后的信號記為pad_signal
    indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T  #相當於對所有幀的時間點進行抽取，得到nf*nw長度的矩陣
    indices=np.array(indices,dtype=np.int32) #將indices轉化為矩陣
    frames=pad_signal[indices] #得到幀信號
    win=np.tile(winfunc,(nf,1))  #window窗函數，這里默認取1
    return frames*win   #返回幀信號矩陣

語音信號的短時時域處理

短時能量和短時平均幅度

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

• 區分濁音和清音段，因為濁音的短時能量$E(i)$比清音大很多；
• 區分聲母和韻母的分界和無話段和有話段的分界

短時平均過零率

　　對於連續語音信號，過零率意味着時域波形通過時間軸，對於離散信號，如果相鄰的取樣值改變符號，則稱為過零。

作用：

發濁音時由於聲門波引起譜的高頻跌落，所以語音信號能量約集中在3kHz以下

發清音時多數能量集中在較高的頻率上，

因為高頻意味着高的短時平均過零率，低頻意味着低的短時平均過零率，所以濁音時具有較低的過零率，而清音時具有較高的過零率。

利用短時平均過零率可以從背景噪聲中找出語音信號，

2、可以用於判斷寂靜無話段與有話段的起點和終止位置。

3、在背景噪聲較小的時候，用平均能量識別較為有效，在背景噪聲較大的時候，用短時平均過零率識別較為有效。

短時自相關函數

短時自相關函數主要應用於端點檢測和基音的提取，在韻母基因頻率整數倍處將出現峰值特性，通常根據除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()

語譜圖

絕大部分信號都可以分解為若干不同頻率的正弦波。
這些正弦波中，頻率最低的稱為信號的基波，其余稱為信號的諧波。
基波只有一個，可以稱為一次諧波，諧波可以有很多個，每次諧波的頻率是基波頻率的整數倍。諧波的大小可能互不相同。
以諧波的頻率為橫坐標，幅值（大小）為縱坐標，繪制的系列條形圖，稱為頻譜。頻譜能夠准確反映信號的內部構造。

　　語譜圖綜合了時域和頻域的特點，明顯的顯示出來了語音頻率隨時間的變化情況，語譜圖的橫軸為時間，縱軸為頻率任意給定頻率成分在給定時刻的強弱用顏色深淺表示。顏色深表示頻譜值大，顏色淺表示頻譜值小，語譜圖上不同的黑白程度形成不同的紋路，稱為聲紋，不用講話者的聲紋是不一樣的，可以用做聲紋識別。

其實得到了分幀信號，頻域變換取幅值，就可以得到語譜圖，如果僅僅是觀察，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('“概率”語譜圖')

語音識別

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']

我對上面這段代碼專門寫了一篇博客來進一步講解和分析，想詳細了解的讀者可以移步https://www.cnblogs.com/LXP-Never/p/11415110.html，語音數據集在這里。

參考文獻

網址：用python做科學計算 http://old.sebug.net/paper/books/scipydoc/index.html#

python標准庫wave模塊https://docs.python.org/3.6/library/wave.html

《python機器學習經典案例》美Prateek Joshi著

傅里葉變換的介紹：http://www.thefouriertransform.com/

各種音階及其對應的頻率 http://pages.mtu.edu/~suits/notefreqs.html

這篇博客的代碼https://github.com/LXP-Neve/Speech-signal-processing

這個網站有很多numpy寫的語音信號處理代碼