信号与系统大作业之Python音频信号处理v1.0使用Python解析Wav文件

发表于 | 分类于 | 阅读次数:
字数统计: 2,882 字 | 阅读时长 ≈ 12 分钟

使用Python解析Wav文件

音频文件及代码已经放在GitHub上了

Wav文件格式

Waveform Audio File Format（WAVE，又或者是因为扩展名而被大众所知的WAV），是微软与IBM公司所开发在个人电脑存储音频流的编码格式，在Windows平台的应用软件受到广泛的支持，地位上类似于麦金塔电脑里的AIFF。此格式属于资源交换档案格式(RIFF)的应用之一，通常会将采用脉冲编码调制的音频资存储在区块中。也是其音乐发烧友中常用的指定规格之一。由于此音频格式未经过压缩，所以在音质方面不会出现失真的情况，但档案的体积因而在众多音频格式中较为大。通常使用三个参数来表示声音，量化位数，取样频率和采样点振幅。量化位数分为8位，16位，24位三种，声道有单声道和立体声之分，单声道振幅数据为n*1矩阵点，立体声为n*2矩阵点，取样频率一般有11025Hz(11kHz) ，22050Hz(22kHz)和44100Hz(44kHz) 三种。

WAV文档遵守资源交换档案格式之规则，在文件的前44(或46)字节放置标头(header)，使播放器或编辑器能够简单掌握文件的基本信息，其内容以区块(chunk)为最小单位，每一区块长度为4字节，而区块之上则由子区块包裹，每一子区块长度不拘，但须在前头先宣告标签及长度(字节)。标头的前3个区块日志案格式及长度；接着第一个子区块包含8个区块，记录声道数量、采样率等信息；接着第二个子区块才是真正的音频数据，长度则视音频长度而定。内容如下表所示。须注意的是，每个区块的端序不尽相同，而音频内容本身则是采用小端序。

DATA块内的信息是根据format chunk内的信息而决定的.由量化位数／声道数／采样率共同决定，下图为四种情况下DATA区域存储信息的格式：

Python读取WAV文件的信息

  • Python自带的Wave库中集成了许多函数，可以打开WAV文件并读出所需要的数据。
    • 可以使用wave.open(file,[mode])打开WAV文件，'r','rb'是只读格式，'w','wb'是只写模式。
    • 使用Wave库中自带的 Wave_read.getparams() 方法可以读取到WAV文件的参数，返回的参数是一个元组(nchannels, sampwidth, framerate, nframes, comptype, compname)
  • 也可以使用scipy.io.wavfile中的read("filename",mmap)返回的值是rate: int , data:numpy array较上一个方法比较简便的是不用再使用numpy.fromstring,numpy.shape,numpy.T等方法对数据进行操作转换成绘图所需要的格式

使用Numpy和Matplotlib来进行波形图的绘制

步骤如下

  • 通过wav库获得WAV文件的头文件中的信息，如采样率／声道数等等.
  • 提取出DATA区域的信息，用Numpy将string格式数据转化为数组
  • 通过判定声道数将DATA区域数据进行处理（对数组矩阵进行转换）
  • 得到每个绘制点的时间（x坐标）
  • 用Matplotlib库的Pyplot API 提供的方法绘制出波形图

代码如下

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
import wave as we
import numpy as np
import matplotlib.pyplot as plt
def wavread(path):
    wavfile = we.open(path,"rb")                        #打开WAV文件
    params = wavfile.getparams()                        #获取WAV头文件中的参数
    framerate,frame_wav= params[2],params[3]            #获取采样率和帧数
    datawav = wavfile.readframes(frame_wav)             #读取每一帧数据
    wavfile.close()                                     #关闭文件
    datause = np.fromstring(datawav,dtype = np.short)   #将String格式转换成int型数据
    datause.shape = -1,2                                #将横向矩阵转为纵向矩阵
    datause = datause.T	                                #矩阵转置
    time = np.arange(0, frame_wav) * (1.0/framerate)    #返回一个时间长度的list
    return datause,time

def main():
    path = input("The Path is:")                        #输入WAV文件路径
    wavdata,wavtime = wavread(path)
    plt.subplot(211)                                    #绘制两行一列图表的第一行
    plt.plot(wavtime, wavdata[0],color = 'green')       #绘制左声道的波形
    plt.subplot(212)                                    #绘制两行一列图表的第二行
    plt.plot(wavtime, wavdata[1])                       #绘制右声道的波形
    plt.show()

main()

使用网上下载的test.wav文件并绘制波形

可以发现这个WAV文件的左右声道波形完全不一样,说明这是个立体声双声道音频文件。

通过查阅scipy的文档，发现了一个specgram函数，可以直接绘制WAV文件的语音频谱图，语音频谱图的横坐标是时间，纵坐标是频率，坐标点值为语音数据能量，因而语音频谱图很好的表达了语音的音色随时间变化的趋势。有些经验丰富的人能够通过看语谱图而知道对应的语音信号的内容， 这种技术成为Spectrogram Reading。

绘制的波形图和语音频谱图如下：

使用自己手机录制的WAV文件进行绘图

图中的波形是录制的“啊”的声音波形。

使用Numpy、Scipy和matplotlib来绘制频谱图

代码如下：

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
import wave as we
import scipy
from scipy.io.wavfile import read
import numpy as np
from scipy.signal import hann
from scipy.fftpack import rfft
import matplotlib.pyplot as plt                     #引入所需要的包

path = input("The Path is:")
wavfile = we.open(path,"rb")                        #打开文件
params = wavfile.getparams()
framerate = params[2]
framerate,frameswav= params[2],params[3]            #读取所需要的参数
datawav = wavfile.readframes(frameswav)
wavfile.close()
time = np.arange(0, frameswav) * (1.0/framerate)
datause = np.fromstring(datawav,dtype = np.short)
window = hann(1024)                                 #使用窗函数，引入汉宁窗
datause.shape = -1,2
datause = datause.T
left = datause[0]
right = datause[1]
left_use = left[0:1024] * window                    #取信号前1024帧乘以窗函数
left_mag = 20 * scipy.log10(abs(rfft(left_use)))    #对左声道的信号做快速傅里叶变换并取对数
left_mag -= max(left_mag)
right_use = right[0:1024] * window
right_mag = 20 * scipy.log10(abs(rfft(right_use)))  #对右声道的信号做快速傅里叶变换并取对数
right_mag -= max(right_mag)
freqs = np.linspace(0,22050,1024)
plt.subplot(211)                                    #绘图
plt.plot(freqs,left_mag,color = 'green')
plt.title("Flute Spectrum")
plt.ylabel("Magnitude(dB)")
plt.xlabel("Frequency Bin")
plt.subplot(212)
plt.plot(freqs,right_mag)
plt.title("Flute Spectrum")
plt.ylabel("Magnitude(dB)")
plt.xlabel("Frequency Bin")
plt.show()

使用test.wav文件绘制的频谱图如下：

使用自己录制的WAV文件绘制的频谱图如下：

使信号通过数字滤波器

数字滤波器

数字滤波器是对数字信号进行滤波处理以得到期望的响应特性的离散时间系统。作为一种电子滤波器，数字滤波器与完全工作在模拟信号域的模拟滤波器不同。数字滤波器工作在数字信号域，它处理的对象是经由采样器件将模拟信号转换而得到的数字信号。数字滤波器具有比模拟滤波器更高的精度，甚至能够实现后者在理论上也无法达到的性能。例如，对于数字滤波器来说很容易就能够做到一个1000Hz的低通滤波器允许999Hz信号通过并且完全阻止1001Hz的信号，模拟滤波器无法区分如此接近的信号。由于奈奎斯特采样定理（Nyquist sampling theorem），数字滤波器的处理能力受到系统采样频率的限制。如果输入信号的频率分量包含超过滤波器1/2采样频率的分量时，数字滤波器因为数字系统的“混叠”而不能正常工作。如果超出1/2采样频率的频率分量不占主要地位，通常的解决办法是在模数转换电路之前放置一个低通滤波器（即抗混叠滤波器）将超过的高频成分滤除。否则就必须用模拟滤波器实现要求的功能。

设计数字滤波器

可以使用scipy.signal包中的firwin(),firwin2()函数来进行FIR滤波器的设计。

可以参考Scipy.Signal.Firwin()的官方文档来设计滤波器，官方文档上也给出了很多滤波器的例子，如下图所示：

使用Hann窗滤波器

Hann窗函数低通滤波器的频率响应图和相位图如下所示：

将WAV的数据文件转换成Numpy中的Array后和滤波器进行卷积，numpy.convolve()方法可以实现信号和滤波器的卷积，经过滤波后的结果如图

与未经过滤波器的波形相对比，可以发现高频部分有一定程度的衰减。

使用Kaiser窗滤波器

Kaiser窗高通滤波器的频率响应和相位图如下图所示：

将WAV文件的数据信号与滤波器进行卷积，经过滤波后的结果如图所示范：

可以很明显的看到，低频部分的信号被极大程度的衰减，可以说是成功完成了滤波。

代码如下

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
import wave as we
import scipy
import numpy as np
from scipy.signal import *
import matplotlib.pyplot as plt

# 读取WAV文件数据
def wavread(path):
    wavfile = we.open(path,'rb')
    params = wavfile.getparams()
    framerate,framewave = params[2],params[3]
    datawav = wavfile.readframes(framewave)
    wavfile.close()
    audio = np.fromstring(datawav,dtype=np.short)
    audio.shape = -1,2
    audio = audio.T
    time = np.arange(0,framewave,1.0/framerate)
    return audio,time,framerate

# 处理WAV文件数据
def fft_filter(signal):
    sampling_rate = 44100
    fft_size = 2200
    win = hann(2200)
    f1 = firwin(101,0.5,window=('kaiser',8),pass_zero = False)
    freqs = np.linspace(0,sampling_rate/2,fft_size/2+1)
    xs = signal[5800:8000] * win
    xc = np.convolve(xs,f1,mode='same')
    xf = np.fft.rfft(xc)/fft_size
    xn = np.fft.rfft(xs)/fft_size
    xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))
    xfn = 20*np.log10(np.clip(np.abs(xn), 1e-20, 1e100))
    return xfp,freqs,xc,xs,xfn

# 绘图
def main():
    path = input("The path is:")
    audio,time,framerate = wavread(path)
    audio_left = audio[0]
    audio_right = audio[1]
    y1,x1,z1,k1,n1 = fft_filter(audio_left)
    y2,x2,z2,k2,n2 = fft_filter(audio_right)
    plt.subplot(421)
    plt.plot(time[5800:8000],k1,color = 'green')
    plt.subplot(422)
    plt.plot(time[5800:8000],k2)
    plt.subplot(423)
    plt.plot(time[5800:8000],z1,color='green')
    plt.title("Filter Output")
    plt.xlabel("Frequency(Hz)")
    plt.subplot(424)
    plt.plot(time[5800:8000],z2)
    plt.title("Filter Output")
    plt.xlabel("Frequency(Hz)")
    plt.subplot(425)
    plt.plot(x1,n1,color='green')
    plt.title("Flute Spectrum Before Filter")
    plt.ylabel("Magnitude(dB)")
    plt.xlabel("Frequency(Hz)")
    plt.subplot(426)
    plt.plot(x2,n2)
    plt.title("Flute Spectrum Before Filter")
    plt.ylabel("Magnitude(dB)")
    plt.xlabel("Frequency(Hz)")
    plt.subplot(427)
    plt.plot(x1,y1,color='green')
    plt.title("Flute Spectrum")
    plt.ylabel("Magnitude(dB)")
    plt.xlabel("Frequency(Hz)")
    plt.subplot(428)
    plt.plot(x2,y2)
    plt.title("Flute Spectrum")
    plt.ylabel("Magnitude(dB)")
    plt.xlabel("Frequency(Hz)")
    plt.show()

main()

绘图

左侧为左声道信号，右侧为右声道信号

  • 第一行所示图片为原始信号
  • 第二行所示图片为经过滤波器后得到的信号
  • 第三行所示图片为原始信号的频域图
  • 第四行所示图片为经过滤波器后信号的频域图

参考