目录
- 背景
- 快速理解
- FFT（快速傅里叶变换）
- IFFT（逆傅里叶变换）
- STFT（短时傅里叶变换）
- 代码实现
- FFT源代码
- IFFT源代码
- FFT、IFFT自己实验
- STFT源代码
- STFT自己实验
- 总结
背景
最近用到了相关操作提取音频信号特征，故在此总结一下几个操作
他们分别是 STFT，FFT，IFFT
快速理解
FFT（快速傅里叶变换）
一张全景照片，包含景区所有风景。
FFT 就像是这张全景照片的处理器，它能快速分析整张照片，告诉你主要的景点在哪里。适合对整个信号（全景照片）进行整体的频率分析。
IFFT（逆傅里叶变换）
一张P好的全景照片，包含景区所有风景。
IFFT 操作就像是P图，给原始图片上添加了各种贴纸，字体，滤镜，可能有很多个图层，当你最后点下了保存按钮，多个图层合成为一张图片的过程就是IFFT过程。
STFT（短时傅里叶变换）
一本相册，一系列连续的照片，每一张照片都展示了风景的一部分。
STFT 就像是分析这本相册，它把每一张照片（时间片段）单独分析，告诉你每一张里有哪些景点。这就能让你知道在不同的时间点，风景（频率成分）是如何变化的。
代码实现
FFT源代码
numpy库作为基准来研究，版本1.24.3，只贴出源码中与FFT操作相关的部分
# 根据不同的 norm 计算不同的归一化因子
def _get_forward_norm(n, norm):
if n < 1:
raise ValueError(f"Invalid number of FFT data points ({n}) specified.")
if norm is None or norm == "backward":
return 1
elif norm == "ortho":
return sqrt(n)
elif norm == "forward":
return n
raise ValueError(f'Invalid norm value {norm}; should be "backward",'
'"ortho" or "forward".')
# 计算之前的前处理函数
def _raw_fft(a, n, axis, is_real, is_forward, inv_norm):
# a: 输入的数组
# n: 变换轴长度，如果指定的 n 小于输入的长度，则会截取输入；如果大于输入的长度，则在末尾用零填充；默认和a一样长。
# axis: 指定计算FFT的轴
# is_real 是否是实数
# is_forward 是否正向FFT
# inv_norm 归一化因子
axis = normalize_axis_index(axis, a.ndim)
if n is None:
n = a.shape[axis]
fct = 1/inv_norm
if a.shape[axis] != n: # 调整输入数组的形状，长的截断 短的补0
s = list(a.shape)
index = [slice(None)]*len(s)
if s[axis] > n:
index[axis] = slice(0, n)
a = a[tuple(index)]
else:
index[axis] = slice(0, s[axis])
s[axis] = n
z = zeros(s, a.dtype.char)
z[tuple(index)] = a
a = z
if axis == a.ndim-1: # 判断变换轴是否为数组的最后一个轴
r = pfi.execute(a, is_real, is_forward, fct) # FFT 变换
else:
a = swapaxes(a, axis, -1)
r = pfi.execute(a, is_real, is_forward, fct)
r = swapaxes(r, axis, -1)
return r
@array_function_dispatch(_fft_dispatcher)
def fft(a, n=None, axis=-1, norm=None):
# a: 输入的数组
# n: 变换轴长度，如果指定的 n 小于输入的长度，则会截取输入；如果大于输入的长度，则在末尾用零填充；默认和a一样长。
# axis: 指定计算FFT的轴
# norm: 用于指定正向/反向变换的归一化模式。
a = asarray(a) # 将输入统一转换为 numpy 数组
if n is None:
n = a.shape[axis]
inv_norm = _get_forward_norm(n, norm) # 计算变换的归一化因子
output = _raw_fft(a, n, axis, False, True, inv_norm) # 计算FFT之前的前处理
return output
这段代码的入口在fft函数，一路看下去发现在我们可见的代码范围内，只是针对输入的数组做了一些预处理，并没有真正FFT计算的部分，真正的计算部分由r = pfi.execute(a, is_real, is_forward, fct)调用，我们跟进去之后代码如下，是由c/c++实现的底层库了：
# encoding: utf-8
# module numpy.fft._pocketfft_internal
# from C:\Users\admin\.conda\envs\pann-cpu-py38\lib\site-packages\numpy\fft\_pocketfft_internal.cp38-win_amd64.pyd
# by generator 1.147
# no doc
# no imports
# functions
def execute(*args, **kwargs): # real signature unknown
pass
# no classes
# variables with complex values
__loader__ = None # (!) real value is '<_frozen_importlib_external.ExtensionFileLoader object at 0x00000177A3437A00>'
__spec__ = None # (!) real value is "ModuleSpec(name='numpy.fft._pocketfft_internal', loader=<_frozen_importlib_external.ExtensionFileLoader object at 0x00000177A3437A00>, origin='C:\\\\Users\\\\admin\\\\.conda\\\\envs\\\\pann-cpu\\\\lib\\\\site-packages\\\\numpy\\\\fft\\\\_pocketfft_internal.cp38-win_amd64.pyd')"
IFFT源代码
这部分同样使用numpy库作为基准来研究，版本1.24.3
# 根据不同的 norm 计算不同的归一化因子
def _get_backward_norm(n, norm):
if n < 1:
raise ValueError(f"Invalid number of FFT data points ({n}) specified.")
if norm is None or norm == "backward":
return n
elif norm == "ortho":
return sqrt(n)
elif norm == "forward":
return 1
raise ValueError(f'Invalid norm value {norm}; should be "backward", '
'"ortho" or "forward".')
# 计算之前的前处理函数
def _raw_fft(a, n, axis, is_real, is_forward, inv_norm):
axis = normalize_axis_index(axis, a.ndim)
if n is None:
n = a.shape[axis]
fct = 1/inv_norm
if a.shape[axis] != n:
s = list(a.shape)
index = [slice(None)]*len(s)
if s[axis] > n:
index[axis] = slice(0, n)
a = a[tuple(index)]
else:
index[axis] = slice(0, s[axis])
s[axis] = n
z = zeros(s, a.dtype.char)
z[tuple(index)] = a
a = z
if axis == a.ndim-1:
r = pfi.execute(a, is_real, is_forward, fct)
else:
a = swapaxes(a, axis, -1)
r = pfi.execute(a, is_real, is_forward, fct)
r = swapaxes(r, axis, -1)
return r
@array_function_dispatch(_fft_dispatcher)
def ifft(a, n=None, axis=-1, norm=None):
# a: 输入的数组
# n: 变换轴长度，如果指定的 n 小于输入的长度，则会截取输入；如果大于输入的长度，则在末尾用零填充；默认和a一样长。
# axis: 指定计算FFT的轴
# norm: 用于指定正向/反向变换的归一化模式。
a = asarray(a) # 将输入统一转换为 numpy 数组
if n is None:
n = a.shape[axis]
inv_norm = _get_backward_norm(n, norm)
output = _raw_fft(a, n, axis, False, False, inv_norm)
return output
通过阅读源码，发现 FFT和 IFFT的代码实现几乎一致（仅有细节差别），都是通过改变传入_raw_fft的参数is_forward来确定做 FFT还是 IFFT。
FFT、IFFT自己实验
这里绘制的结果是三部分，分别是原始波形、FFT结果、IFFT结果
import librosa
import numpy as np
import matplotlib.pyplot as plt
# 加载音频文件
file_path = r'C:\Users\admin\Desktop\自己的分类数据\5.0\unbalanced\normal_audio\14-2024.02.25.10.48.37_(2.836402877697843, 7.836402877697843).wav'
y, sr = librosa.load(file_path, sr=None)
# 计算 FFT
fft_result = np.fft.fft(y)
# 计算 IFFT
ifft_result = np.fft.ifft(fft_result)
# 绘制原始音频信号
plt.figure(figsize=(14, 5))
plt.subplot(3, 1, 1)
plt.title('Original Audio Signal')
plt.plot(y)
plt.xlabel('Time')
plt.ylabel('Amplitude')
# 绘制 FFT 结果（只取前一半，因为FFT对称）
plt.subplot(3, 1, 2)
plt.title('FFT Result')
plt.plot(np.abs(fft_result)[:len(fft_result) // 2])
plt.xlabel('Frequency Bin')
plt.ylabel('Magnitude')
# 绘制 IFFT 结果
plt.subplot(3, 1, 3)
plt.title('Reconstructed Signal from IFFT')
plt.plot(ifft_result.real)
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.tight_layout()
# 保存图片
plt.savefig('fft_comparison.png')
STFT源代码
librosa库为基准来研究，版本0.9.1，仔细阅读学习了一下源代码，发现有很多写法还是十分精妙的，怪不得计算效率会高很多。
@deprecate_positional_args
@cache(level=20)
def stft(
y, # 待处理音频信号
*, # 分隔符，此分隔符之后的参数都要键值对输入
n_fft=2048, # * 在时间维度上的窗口大小，用于FFT计算
hop_length=None, # 每个帧之间的跳跃长度 默认n_fft / 4
win_length=None, # * 在音频频率维度上的窗口大小，用于窗口函数计算
window="hann", # 选用的窗口函数：对当前窗口信号进行加权处理，减少频谱泄漏和边界效应
center=True, # 控制每个帧是否在中心对齐，否则左边界对齐
dtype=None, # 输出数组的数据类型
pad_mode="constant", # 填充模式，用于确定如何处理信号边缘
):
# By default, use the entire frame 默认窗口长度 win_length 的设置
if win_length is None:
win_length = n_fft
# Set the default hop, if it's not already specified 默认跳跃长度 hop_length 的设置
if hop_length is None:
hop_length = int(win_length // 4)
# Check audio is valid 音频有效性检查
util.valid_audio(y, mono=False)
# 获取窗口函数 fftbins为True则会将窗口函数填充到 win_length 相同长度
fft_window = get_window(window, win_length, fftbins=True)
# Pad the window out to n_fft size 将窗口函数填充到 n_fft相同长度（为了处理win_length != n_fft的情形）
fft_window = util.pad_center(fft_window, size=n_fft)
# Reshape so that the window can be broadcast 重新构造fft_window的形状以便和 y 进行计算
# ndim： fft_window 扩展之后的总维度数，axes：不变的轴（将数据填充到-2轴以外的其他位置上，一般音频信号-2是时间轴）
# eg. y(2,4096) fft_window(1024,)，设置 ndim=3, axes=-2
# 因此扩展后的总维度为3，同时fft_window需要-2维度保持原有的值不变，其余位置进行广播，得到(2, 1024, 4096)
fft_window = util.expand_to(fft_window, ndim=1 + y.ndim, axes=-2)
# Pad the time series so that frames are centered 中心填充
if center:
if n_fft > y.shape[-1]:
warnings.warn(
"n_fft={} is too small for input signal of length={}".format(
n_fft, y.shape[-1]
),
stacklevel=2,
)
padding = [(0, 0) for _ in range(y.ndim)]
padding[-1] = (int(n_fft // 2), int(n_fft // 2))
y = np.pad(y, padding, mode=pad_mode)
elif n_fft > y.shape[-1]: # 如果 n_fft 大于输入信号的长度，报错
raise ParameterError(
"n_fft={} is too large for input signal of length={}".format(
n_fft, y.shape[-1]
)
)
# Window the time series. 按照n_fft和hop_length进行滑动窗口从时间维度切分y，y_frames(num_frames, frame_length)
# 补充说明：
# STFT 通常会生成一个三维矩阵，维度一般为 (batch_size, num_frames, num_bins)
# batch_size：批次大小，表示处理多个音频片段
# num_frames：时间帧数，表示音频片段被分割成多少个时间窗口
# num_bins(frame_length)：频率成分数，表示每个时间窗口内计算的频率分量数
y_frames = util.frame(y, frame_length=n_fft, hop_length=hop_length)
fft = get_fftlib() # 获取FFT库对象，便于后续操作
if dtype is None: # 设置数据类型
dtype = util.dtype_r2c(y.dtype)
# Pre-allocate the STFT matrix 预分配STFT矩阵
shape = list(y_frames.shape)
shape[-2] = 1 + n_fft // 2 # 因为FFT结果是镜像对称，所以存储空间减半，+1是因为FFT结果是复数
stft_matrix = np.empty(shape, dtype=dtype, order="F") # 预分配内存，empty函数创建的数组所有元素都没有默认值
# 1 np.prod(stft_matrix.shape[:-1]) * stft_matrix.itemsize
# 计算出所有时间窗口和所有批次中，一个频率成分所需的总内存（一列）
# stft_matrix.itemsize 返回 stft_matrix 中每个元素的字节大小
# stft_matrix.shape[:-1] 返回除了最后一个维度的剩余形状，eg.(1024,256)[:-1] 返回 1024
# np.prod(...) 求出...的所有维度的乘积（总元素个数）
# 2 计算出需要的列数 util.MAX_MEM_BLOCK // ...
# util.MAX_MEM_BLOCK 代表限定的最大可用内存大小
# 补充说明：
# 计算结果是列数的原因：在stft中，横轴x代表时间，纵轴y代表频率，因此每一列代表一个时间窗口的内存大小
n_columns = util.MAX_MEM_BLOCK // (
np.prod(stft_matrix.shape[:-1]) * stft_matrix.itemsize
)
n_columns = max(n_columns, 1)
# 分块计算
for bl_s in range(0, stft_matrix.shape[-1], n_columns):
bl_t = min(bl_s + n_columns, stft_matrix.shape[-1]) # 结束索引
# 计算当前块的STFT
stft_matrix[..., bl_s:bl_t] = fft.rfft(
fft_window * y_frames[..., bl_s:bl_t], axis=-2
)
return stft_matrix # 横轴是时间，纵轴是频率
STFT自己实验
了解这些之后自己写个例子试试看，stft_result 是得到的直接计算结果，D是全部求绝对值之后的结果，DB是为了更好的可视化做的后处理，我加上原始波形图，三者放在一起比较，保存了一张结果。
import librosa.display
import matplotlib.pyplot as plt
import numpy as np
# 加载音频文件
audio_file = r'C:\Users\admin\Desktop\test.wav'
y, sr = librosa.load(audio_file)
# 参数设置
n_fft = 2048 # FFT窗口大小
hop_length = 512 # 帧之间的跳跃长度
# 计算STFT
stft_result = librosa.stft(y, n_fft=n_fft, hop_length=hop_length)
# 获取幅度谱
D = np.abs(stft_result)
# 将幅度谱转换为分贝单位
DB = librosa.amplitude_to_db(D, ref=np.max)
# 绘制原始波形
plt.figure(figsize=(10, 8))
plt.subplot(3, 1, 1)
librosa.display.waveshow(y, sr=sr)
plt.title('Waveform')
plt.xlabel('Time')
# 绘制原始幅度谱 D
plt.subplot(3, 1, 2)
librosa.display.specshow(D, sr=sr, hop_length=hop_length, x_axis='time', y_axis='log')
plt.colorbar(format='%+2.0f')
plt.title('STFT Magnitude Spectrum (D)')
plt.xlabel('Time')
plt.ylabel('Frequency')
# 绘制分贝单位的幅度谱 DB
plt.subplot(3, 1, 3)
librosa.display.specshow(DB, sr=sr, hop_length=hop_length, x_axis='time', y_axis='log')
plt.colorbar(format='%+2.0f dB')
plt.title('STFT Magnitude Spectrum (DB)')
plt.xlabel('Time')
plt.ylabel('Frequency')
plt.tight_layout()
# 保存图片
plt.savefig('stft_comparison.png')
总结
FFT是时域信号——>频域信号，IFFT是频域信号——>时域信号，输入是一维音频信号的前提下，这两个操作得到的结果都是一维的特征。
一般 IFFT都会配合 FFT操作一起使用，先使用FFT，然后从频率角度处理音频信号，最后使用 IFFT操作重新将信号恢复。
STFT是时域信号——>时间-频率域信号，输入是一维音频信号的前提下，这个操作得到的结果都是二维的特征。
二者的使用场景根本区别在于是否需要了解频率随着时间的变化，如果需要就使用 STFT，否则使用 FFT + IFFT
