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
import numpy as np
import scipy
import matplotlib.pyplot as plt

# 获取音频的时域信号序列
wav=wave.open('test.wav','rb')
print(wav.getframerate(),wav.getnchannels(),wav.getnframes(),wav.getsampwidth())
c=wav.readframes(int(0.025*wav.getframerate()))
npc=np.frombuffer(c,dtype=np.short)
npc=npc.reshape((-1,2)).T[0]
print(npc,npc.shape,int(0.025*wav.getframerate()))

n=int(0.025*wav.getframerate())
plt.plot(np.arange(n),npc)
plt.show()
fr=wav.getframerate()
npcf=scipy.fft.fft(npc)
abs_npc=np.abs(npcf) # 振幅
angle_npc=np.angle(npcf) # 相位
x=np.arange(n/2)*fr*2/n**2
x[0]=x[0]*2

plt.plot(x,abs_npc[:n//2])
plt.show()

这里重点介绍快速傅里叶变换，即scipy.fft.fft()函数的原理。

首先，我们手里有一段时域信号序列，这段序列的长度为n，音频的采样率为fr，对这个序列进行fft，返回复数数组，每个复数表示一个正弦波，含有振幅、相位和频率三个信息。npcf的长度也为n，且是个复数数组。npcf[1:n/2]表示频域的各个正弦波，npc中信号的最大周期T=n/fr，则获知的信号的基频（最小频率）为fr/n，根据奈奎斯特采样定理，还原的信号频率最大为fr/2，npcf[1:n/2]代表1*fr/n到n/2*fr/n的频率的正弦波，npcf[0]表示直流分量，算是偏置项吧。npcf[i]=real+j*imag，则正弦波的相位就是复数的幅角，arctan(imag/real)。振幅是复数的模，amplitude=\sqrt{real^2+imag^2}，但是fft的返回值的模是放大值，直流分量的振幅放大了N倍，弦波分量的振幅放大了N/2倍

数学原理：
对长度为N的时序序列x，可通过离散傅里叶变换得到长度为N的X，X_k=\sum_{n=0}^{N-1}x_n*e^{-i2πkn/N}=\sum_{n=0}^{N-1}x_n*[cos(2πkn/N)-i*sin(2πkn/N)]，k=0,1,2,3,\cdots,N-1

参考文献：Python机器学习（五十六）SciPy fftpack(傅里叶变换)

离散余弦变换(DCT)和离散傅里叶变换(DFT)的关系

创建于2023.3.17/18.1，修改于2023.3.18/10.27