Python对音频进行频谱分析

一、储备函数

1.1 计时装饰器、查找最接近的值、音高频率创建与读取

import time
import math
import numpy as np
import pandas as pd


# 定义时间装饰器
def timeit(func):
    def wrapper(*args, **kwargs):
        start_time = time.time()
        result = func(*args, **kwargs)
        using = (time.time() - start_time) * 1000
        print(f'运行时间：{using} 毫秒')
        return result
    return wrapper


# 普通未排序数组查找最接近的值
@timeit
def find_nearest(array, value):
    array = np.array(array)
    idx = (np.abs(array-value)).argmin()
    return array[idx]


# 已排序数组查找最接近的值
@timeit
def find_nearest_sorted(array, value):
    idx = np.searchsorted(array, value, side="left")
    if idx > 0 and (idx == len(array) or math.fabs(value - array[idx-1]) < math.fabs(value - array[idx])):
        return array[idx-1]
    else:
        return array[idx]


# 计算音高频率，参数一般为440、442或432
def calculate_frequency(standard_frequency):
    list_frequency = []
    for i in range(132):
        fre = round(standard_frequency / 32 * math.pow(2, (i-9.0) / 12),3)
        list_frequency.append(fre)
    return list_frequency[12:]


# 利用计算音高创建音高频率df
@timeit
def create_frequency():
    df_frequency = pd.DataFrame()
    df_frequency['letter_name'] = ['C0', 'C♯/D♭0', 'D0', 'D♯/E♭0', 'E0', 'F0', 'F♯/G♭0', 'G0', 'G♯/A♭0', 'A0', 'A♯/B♭0', 'B0', 'C1', 'C♯/D♭1', 'D1', 'D♯/E♭1', 'E1', 'F1', 'F♯/G♭1', 'G1', 'G♯/A♭1', 'A1', 'A♯/B♭1', 'B1', 'C2', 'C♯/D♭2', 'D2', 'D♯/E♭2', 'E2', 'F2', 'F♯/G♭2', 'G2', 'G♯/A♭2', 'A2', 'A♯/B♭2', 'B2', 'C3', 'C♯/D♭3', 'D3', 'D♯/E♭3', 'E3', 'F3', 'F♯/G♭3', 'G3', 'G♯/A♭3', 'A3', 'A♯/B♭3', 'B3', 'C4', 'C♯/D♭4', 'D4', 'D♯/E♭4', 'E4', 'F4', 'F♯/G♭4', 'G4', 'G♯/A♭4', 'A4', 'A♯/B♭4', 'B4', 'C5', 'C♯/D♭5', 'D5', 'D♯/E♭5', 'E5', 'F5', 'F♯/G♭5', 'G5', 'G♯/A♭5', 'A5', 'A♯/B♭5', 'B5', 'C6', 'C♯/D♭6', 'D6', 'D♯/E♭6', 'E6', 'F6', 'F♯/G♭6', 'G6', 'G♯/A♭6', 'A6', 'A♯/B♭6', 'B6', 'C7', 'C♯/D♭7', 'D7', 'D♯/E♭7', 'E7', 'F7', 'F♯/G♭7', 'G7', 'G♯/A♭7', 'A7', 'A♯/B♭7', 'B7', 'C8', 'C♯/D♭8', 'D8', 'D♯/E♭8', 'E8', 'F8', 'F♯/G♭8', 'G8', 'G♯/A♭8', 'A8', 'A♯/B♭8', 'B8', 'C9', 'C♯/D♭9', 'D9', 'D♯/E♭9', 'E9', 'F9', 'F♯/G♭9', 'G9', 'G♯/A♭9', 'A9', 'A♯/B♭9', 'B9']
    df_frequency['frequency_440'] = calculate_frequency(440)
    df_frequency['frequency_442'] = calculate_frequency(442)
    return df_frequency


# 利用读取音高创建音高频率df
@timeit
def read_frequency():
    df = pd.read_excel('f:/音高频率.xlsx', index_col='音名')
    letter_name = []
    for i in range(10):
        for j in df.index.to_list():
            letter_name.append(j.strip() + str(i))
    df_frequency = pd.DataFrame()
    df_frequency['frequency_440'] = df.stack().sort_values()
    df_frequency['letter_name'] = letter_name
    df_frequency.reset_index(inplace=True)
    df_frequency = df_frequency[['letter_name', 'frequency_440']]
    return df_frequency


df_frequency = create_frequency()   # 2.986431121826172 毫秒
# df_frequency = read_frequency()     # 378.781795501709 毫秒

df_frequency.to_excel('f:/frequency.xlsx')
print(df_frequency.loc[df_frequency['letter_name']=='A4'])

1.2 440音高频率表

音名	0	1	2	3	4	5	6	7	8	9
C	16.352	32.703	65.406	130.81	261.63	523.25	1046.5	2093	4186	8372
C♯/D♭	17.324	34.648	69.296	138.59	277.18	554.37	1108.7	2217.5	4434.9	8869.8
D	18.354	36.708	73.416	146.83	293.66	587.33	1174.7	2349.3	4698.6	9397.3
D♯/E♭	19.445	38.891	77.782	155.56	311.13	622.25	1244.5	2489	4978	9956.1
E	20.602	41.203	82.407	164.81	329.63	659.26	1318.5	2637	5274	10548
F	21.827	43.654	87.307	174.61	349.23	698.46	1396.9	2793.8	5587.7	11175
F♯/G♭	23.125	46.249	92.499	185	369.99	739.99	1480	2960	5919.9	11840
G	24.5	48.999	97.999	196	392	783.99	1568	3136	6271.9	12544
G♯/A♭	25.957	51.913	103.83	207.65	415.3	830.61	1661.2	3322.4	6644.9	13290
A	27.5	55	110	220	440	880	1760	3520	7040	14080
A♯/B♭	29.135	58.27	116.54	233.08	466.16	932.33	1864.7	3729.3	7458.6	14917
B	30.868	61.735	123.47	246.94	493.88	987.77	1975.5	3951.1	7902.1	15804

对采样率为44.1kHz的AAC音频进行解码时，一帧的解码时间须控制在23.22毫秒内。

二、Pyaudio

我们可以使用Python自带的标准库wave。
但实际我并不推荐使用wave，用soundfile, librosa和pydub之类的第三方库会方便和强大得多

初衷 语音识别领域对音频文件进行频谱分析是一项基本的数据处理过程，同时也为后续的特征分析准备数据。

前驱知识
Python需要使用的相关库

wave
https://docs.python.org/3/library/wave.html
pyaudio
http://people.csail.mit.edu/hubert/pyaudio/
numpy
https://www.runoob.com/numpy/numpy-tutorial.html
pylab
https://www.programcreek.com/python/example/2345/pylab.title
音频帧概率详解
1.采样率（Sample Rate）：每秒从连续信号中提取并组成离散信号的采样个数，它用赫兹（Hz）来表示。一般音乐CD的采样率是44100Hz，所以视频编码中的音频采样率保持在这个级别就完全足够了，通常视频转换器也将这个采样率作为默认设置。
2.帧率（Frame rate）：是用于测量显示帧数的量度。所谓的测量单位为每秒显示帧数(Frames per Second，简称：FPS）或“赫兹”（Hz）。
3.码率（Bit Rate）：指视频或音频文件在单位时间内使用的数据流量，该参数的单位通常是Kbps，也就是千比特每秒。通常2000kbps~3000kbps就已经足以将画质效果表现到极致了。码率参数与视频文件最终体积大小有直接性的关系
4.正常人听觉的频率范围大约在20Hz~20kHz之间，根据奈奎斯特采样理论，为了保证声音不失真，采样频率应该在40kHz左右。常用的音频采样频率有8kHz、11.025kHz、22.05kHz、16kHz、37.8kHz、44.1kHz、48kHz等，如果采用更高的采样频率，还可以达到DVD的音质

背景知识:
(一个AAC原始帧包含一段时间内1024个采样及相关数据)
分析：
1.AAC
音频帧的播放时间=一个AAC帧对应的采样样本的个数/采样频率(单位为s)
一帧 1024个 sample。采样率 Samplerate 44.1KHz，每秒44100个sample, 所以根据公式 音频帧的播放时间=一个AAC帧对应的采样样本的个数/采样频率
当前AAC一帧的播放时间是= 1024*1000/44100= 22.32ms(单位为ms)
2.MP3
mp3 每帧均为1152个字节， 则：
frame_duration = 1152 * 1000 / sample_rate
例如：sample_rate = 44100HZ时，计算出的时长为26.122ms，这就是经常听到的mp3每帧播放时间固定为26ms的由来。
3.H264
视频的播放时间跟帧率有关:
frame_duration = 1000/帧率(fps)
例如：fps = 25.00 ，计算出来的时常为40ms，这就是同行所说的40ms一帧视频数据。

# 打开提前准备的WAV文档，文件路径根据需要做修改

wf = wave.open("文档路径", "rb")

# 创建PyAudio对象
p = pyaudio.PyAudio()
stream = p.open(format=p.get_format_from_width(wf.getsampwidth()),
channels=wf.getnchannels(),
rate=wf.getframerate(),
output=True
nframes = wf.getnframes()
framerate = wf.getframerate()


# 读取完整的帧数据到str_data中，这是一个string类型的数据
str_data = wf.readframes(nframes)
wf.close()


# 将音频波形数据转换为数组
# A new1-D array initialized from raw binary ortext data in a string.
wave_data = numpy.fromstring(str_data, dtype=numpy.short)


# 将wave_data数组改为2列，行数自动匹配。在修改shape的属性时，需使得数组的总长度不变。
wave_data.shape = -1,2


# 将数组转置
wave_data = wave_data.T
#time 也是一个数组，与wave_data[0]或wave_data[1]配对形成系列点坐标
#time = numpy.arange(0,nframes)*(1.0/framerate)


# 绘制波形图
#pylab.plot(time, wave_data[0])
#pylab.subplot(212)
#pylab.plot(time, wave_data[1], c="g")
#pylab.xlabel("time (seconds)")
#pylab.show()

# 采样点数，修改采样点数和起始位置进行不同位置和长度的音频波形分析
N=44100
start=0 #开始采样位置
df = framerate/(N-1) # 分辨率
freq = [df*n for n in range(0,N)] #N个元素
wave_data2=wave_data[0][start:start+N]
c=numpy.fft.fft(wave_data2)*2/N


# 常规显示采样频率一半的频谱
d=int(len(c)/2)


# 仅显示频率在4000以下的频谱
while freq[d]>4000:
d=10
pylab.plot(freq[:d-1],abs(c[:d-1]),'r')
pylab.show()