【音频特征提取】傅里叶变换算法源码学习记录

背景

最近用到了相关操作提取音频信号特征，故在此总结一下几个操作
他们分别是 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 1elif norm == "ortho":return sqrt(n)elif norm == "forward":return nraise 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_normif a.shape[axis] != n: # 调整输入数组的形状，长的截断 短的补0s = 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] = nz = zeros(s, a.dtype.char)z[tuple(index)] = aa = zif 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# functionsdef execute(*args, **kwargs): # real signature unknownpass# 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 nelif norm == "ortho":return sqrt(n)elif norm == "forward":return 1raise 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_normif 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] = nz = zeros(s, a.dtype.char)z[tuple(index)] = aa = zif 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 / 4win_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]) # 结束索引# 计算当前块的STFTstft_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')