短时傅里叶变换（STFT）实例

记录学习音频短时傅里叶变换的过程，注意，这里并不会讲述傅里叶变换的原理
想了解原理可以自行搜索或查看文中的参考文章

运行使用Pycharm，环境为python3.10

1.导入软件包，设置快速傅里叶变换的参数

Win_size是对音频滑窗处理的窗体大小，Step_size是滑窗步长，Data_logs用于存放傅里叶变化后的值

（Win_size和Step_size的设置没有标准，但窗体越大，步长越短，结果所包含的频域信息就越多，时域信息就越少，参考文章对音频信号作短时傅里叶变换）

import wave
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.io import wavfile
from scipy.fftpack import fft

Win_size = 16384
Step_size = 8192
Data_logs = []

2.导入wav

利用scipy库的wavfile读入wav，Fs为采样频率，Data为振幅数据

此处的testro文件是作者从耳聆网下载的一段人和Siri对话的声音。文件下载后是MP3文件，使用这个网站可以免费转换为wav文件

wavFile = "testro.wav"
Fs, Data = wavfile.read(wavFile)

以下代码用于获取音频的统计信息，可以删去（代码参考了文章http://t.csdn.cn/jt8Ux）

f = wave.open(wavFile)
Channels = f.getnchannels()
SampleRate = f.getframerate()
bit_type = f.getsampwidth() * 8
frames = f.getnframes()
# Duration 也就是音频时长 = 采样点数/采样率
Duration = wav_time = frames / float(SampleRate)  # 单位为s

print("通道数(Channels)：", Channels)
print("采样率(SampleRate)：", SampleRate)
print("比特(Precision)：", bit_type)
print("采样点数(frames)：", frames)
print("帧数或者时间(Duration)：", Duration)
print("Data的统计信息如下：\n", pd.DataFrame(Data).describe())  # 统计数据信息

统计信息输出如下：

音频头参数： _wave_params(nchannels=1, sampwidth=2, framerate=48000, nframes=2477999, comptype='NONE', compname='not compressed')
通道数(Channels)： 1
采样率(SampleRate)： 48000
比特(Precision)： 16
采样点数(frames)： 2477999
帧数或者时间(Duration)： 51.62497916666667

Data的统计信息如下：
count  2.477999e+06
mean   1.589993e-03
std    7.394472e+02
min   -1.074700e+04
25%   -1.160000e+02
50%    0.000000e+00
75%    1.180000e+02
max    9.835000e+03

3.快速傅里叶变换

对音频数据不断滑窗并进行傅里叶变换

由于FFT后获取的结果是一个a+bi形式的向量，所以需要对其取模（abs）以获取幅频特性，参考文章https://www.jianshu.com/p/78043f8f8306

lens = Data.shape[0]
i = 0
while i + Win_size < lens:  # 防止列表出界
    win = Data[i:i + Win_size]
    FFT = fft(win)  # 快速傅里叶变换
    Data_logs.append(abs(FFT))  # 取模
    i = i + Step_size  # 滑向下一个窗口

4.变换后数据处理

将列表形式的数据转化为ndarray，对数据转置以便画图（因为一般来说频率是作为纵轴绘制的）

由于DFT的频谱是用谱密度定义的，即它的幅值表示的是单位带宽的幅值，所以FFT后幅值要除以N/2，参考文章FFT结果的物理意义

Data_logs = np.array(Data_logs).T

Data_logs = Data_logs / Win_size * 2
Data_logs[0] /= 2

然后截取一半的数据（因为FFT的结果是左右对称的，参考文章FFT的详细解释）
（顺便取个对数）

xlim = int(Data_logs.shape[0] / 2)
Data_logs = 20 * np.log10(Data_logs[0:xlim, :])

5.绘图

到这里就已经获得了音频的傅里叶结果了，最后使用plt绘制语谱图

imshow的参数解释：
extent：缩放横纵坐标值以符合实际，各个数字代表的含义为 [横轴最左边的值，横轴最右边的值，纵轴最下面的值，纵轴最上面的值]（参考文章imshow的origin和extent）
origin：原图的频率是0在最上方，从上往下的，不符合常规，设置“lower”可以上下翻转图片
aspect：自动控制轴的纵横比
至于为什么纵轴最上面的值是Fs//2（//的意思是取除后的整数部分），参考知乎链接

plt.imshow(Data_logs, extent=[0, len(Data) / Fs, 0, Fs // 2], origin='lower', aspect='auto')
plt.colorbar()
plt.xlabel('time(s)')
plt.ylabel('frequence(Hz)')
plt.title('STFT')
plt.show()

6.成图

如下

7.总结

从结果图像来看，低频段拥有更高的分贝，这是符合常理的。但频率在15500左右出现了断层，这是值得分析的一个点。

另外在查找资料的时候发现似乎使用汉明窗（Hamming）进行滑窗可以更好的获取频率信息，这是一个改进的方向，参考Hamming(汉明)窗的原理介绍及实例解析

写文章比较仓促，如果发现了任何错误欢迎交流！

附整体代码：

import wave
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.io import wavfile
from scipy.fftpack import fft

Win_size = 16384
Step_size = 8192
Data_logs = []

wavFile = "testro.wav"
Fs, Data = wavfile.read(wavFile)

# 统计数据信息
f = wave.open(wavFile)
Channels = f.getnchannels()
SampleRate = f.getframerate()
bit_type = f.getsampwidth() * 8
frames = f.getnframes()
# Duration 也就是音频时长 = 采样点数/采样率
Duration = wav_time = frames / float(SampleRate)  # 单位为s
print("通道数(Channels)：", Channels)
print("采样率(SampleRate)：", SampleRate)
print("比特(Precision)：", bit_type)
print("采样点数(frames)：", frames)
print("帧数或者时间(Duration)：", Duration)
print("Data的统计信息如下：\n", pd.DataFrame(Data).describe())

lens = Data.shape[0]
i = 0
while i + Win_size < lens:
    win = Data[i:i + Win_size]
    FFT = fft(win)
    Data_logs.append(abs(FFT))
    i = i + Step_size

Data_logs = np.array(Data_logs).T

Data_logs = Data_logs / Win_size * 2
Data_logs[0] /= 2

xlim = int(Data_logs.shape[0] / 2)
Data_logs = 20 * np.log10(Data_logs[0:xlim, :])

plt.imshow(Data_logs, extent=[0, len(Data) / Fs, 0, Fs // 2], origin='lower', aspect='auto')
plt.colorbar()
plt.xlabel('time(s)')
plt.ylabel('frequence(Hz)')
plt.title('STFT')
plt.show()