【音頻特征提取】傅里葉變換算法源碼學習記錄

背景

最近用到了相關操作提取音頻信號特征，故在此總結一下幾個操作
他們分別是 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')