记录学习音频短时傅里叶变换的过程，注意，这里并不会讲述傅里叶变换的原理
想了解原理可以自行搜索或查看文中的参考文章
运行使用Pycharm，环境为python3.10
1.导入软件包，设置快速傅里叶变换的参数
Win_size是对音频滑窗处理的窗体大小，Step_size是滑窗步长，Data_logs用于存放傅里叶变化后的值
（Win_size和Step_size的设置没有标准，但窗体越大，步长越短，结果所包含的频域信息就越多，时域信息就越少，参考文章对音频信号作短时傅里叶变换）
import waveimport numpy as npimport pandas as pdfrom matplotlib import pyplot as pltfrom scipy.io import wavfilefrom scipy.fftpack import fftWin_size = 16384Step_size = 8192Data_logs = []
2.导入wav
利用scipy库的wavfile读入wav，Fs为采样频率，Data为振幅数据
此处的testro文件是作者从耳聆网下载的一段人和Siri对话的声音。文件下载后是MP3文件，使用这个网站可以免费转换为wav文件
wavFile = "testro.wav"Fs, Data = wavfile.read(wavFile)
以下代码用于获取音频的统计信息，可以删去（代码参考了文章/jt8Ux）
f = wave.open(wavFile)Channels = f.getnchannels()SampleRate = f.getframerate()bit_type = f.getsampwidth() * 8frames = f.getnframes()# Duration 也就是音频时长 = 采样点数/采样率Duration = wav_time = frames / float(SampleRate) # 单位为sprint("通道数(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.62497916666667Data的统计信息如下：count 2.477999e+06mean 1.589993e-03std 7.394472e+02min -1.074700e+0425% -1.160000e+0250% 0.000000e+0075% 1.180000e+02max 9.835000e+03
3.快速傅里叶变换
对音频数据不断滑窗并进行傅里叶变换
由于FFT后获取的结果是一个a+bi形式的向量，所以需要对其取模（abs）以获取幅频特性，参考文章/p/78043f8f8306
lens = Data.shape[0]i = 0while 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).TData_logs = Data_logs / Win_size * 2Data_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 waveimport numpy as npimport pandas as pdfrom matplotlib import pyplot as pltfrom scipy.io import wavfilefrom scipy.fftpack import fftWin_size = 16384Step_size = 8192Data_logs = []wavFile = "testro.wav"Fs, Data = wavfile.read(wavFile)# 统计数据信息f = wave.open(wavFile)Channels = f.getnchannels()SampleRate = f.getframerate()bit_type = f.getsampwidth() * 8frames = f.getnframes()# Duration 也就是音频时长 = 采样点数/采样率Duration = wav_time = frames / float(SampleRate) # 单位为sprint("通道数(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 = 0while i + Win_size < lens:win = Data[i:i + Win_size]FFT = fft(win)Data_logs.append(abs(FFT))i = i + Step_sizeData_logs = np.array(Data_logs).TData_logs = Data_logs / Win_size * 2Data_logs[0] /= 2xlim = 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()
