Open In Colab 在 Colab 上执行或在 GitHub 上查看/下载此 notebook

傅里叶变换与频谱图

在语音和音频处理中，时域信号通常会被转换为另一个域。但是，为什么我们需要转换音频信号呢？

信号的一些语音特征/模式（例如，基频 (pitch)，共振峰 (formants)）在时域下查看时可能不太明显。通过适当设计的变换，可能更容易从信号本身提取所需的信息。

最流行的变换是傅里叶变换 (Fourier Transform)，它将时域信号转换为频域 (frequency domain) 中的等效表示。在接下来的章节中，我们将描述傅里叶变换以及其他相关的变换，如短时傅里叶变换 (Short-Term Fourier Transform, STFT) 和频谱图 (spectrograms)。

1. 傅里叶变换

时间离散序列 \(f[n]={f[0],f[1],..f[N-1]}\) 的傅里叶变换称为离散傅里叶变换 (Discrete Fourier Transform, DFT)，其定义如下

\(F_{k} = \sum_{n=0}^{N-1} f_{n} e^{-j\frac{2\pi}{N}kn}\)

逆变换称为逆离散傅里叶变换 (Inverse Discrete Fourier Transform, IDFT)，它将频域信号 \(F_k\) 映射回时域信号 \(f_n\)

\(f_{n} = \sum_{k=0}^{N-1} F_{k} e^{j\frac{2\pi}{N}kn}\)

这两种表示是等价的，应用它们不会丢失信息。这只是表示同一信号的不同方式。

直觉是什么？

傅里叶变换背后的思想是将信号表示为一系列频率递增的复正弦波的加权和。例如，复指数项 \(e^{j\frac{2\pi}{N}kn}\) 决定了这种“复正弦波”的频率

\(e^{j\frac{2\pi}{N}kn} = cos(\frac{2\pi}{N}kn) +j sin(\frac{2\pi}{N}kn)\).

而项 \(F_{k}\) 则是另一个复数，它决定了频率分量的幅度和移位（相位）。可以证明，使用 N 个具有适当幅度 (amplitude) 和相位 (phase) 的复正弦波，我们可以建模任何信号。换句话说，复正弦波是构成您信号的基本砖块。如果您像搭乐高积木一样恰当地组合它们，您就可以创建任何想要的信号（周期性的和非周期性的都可以）。

这种变换的复杂度是 \(O(N^2)\)，因为对于频率表示 \(F_k\) 的每个元素 k，我们都必须遍历序列的所有 N 个元素。这使得计算长序列的 DFT 和 IDFT 变得不可能。

幸运的是，存在名为快速傅里叶变换 (Fast-Fourier Transform, FFT) 的算法，其计算复杂度为 \(O(Nlog(N))\)。FFT 将输入序列分割成小块，然后组合它们的 DFT。

“复正弦波”这个概念可能很难理解。不过，您可以在网上找到许多包含精美图形动画的优秀资料来帮助理解（参见参考文献中的教程）。目前，我们只需将傅里叶变换视为一种将实值序列映射为复值序列的线性变换。

在计算一些 DFT 之前，让我们先下载一些语音信号并安装 speechbrain

%%capture
!wget https://www.dropbox.com/s/u8qyvuyie2op286/spk1_snt1.wav

%%capture
# Installing SpeechBrain
BRANCH = 'develop'
!git clone https://github.com/speechbrain/speechbrain.git -b $BRANCH
%cd /content/speechbrain/
!python -m pip install .

import torch
import matplotlib.pyplot as plt
from speechbrain.dataio.dataio import read_audio

signal = read_audio('/content/spk1_snt1.wav')
print(signal.shape)

# fft computation
fft = torch.fft.fft(signal.squeeze(), dim=0)
print(fft)
print(fft.shape)

正如您所见，输入信号是实数（因此虚部填充为零）。DFT 是一个张量，其中包含变换的实部和虚部。

现在我们计算 DFT 的幅度和相位并绘制出来

# Real and Imaginary parts
real_fft = fft.real
img_fft = fft.imag

mag = torch.sqrt(torch.pow(real_fft,2) + torch.pow(img_fft,2))
phase = torch.arctan(img_fft/real_fft)

plt.subplot(211)
x_axis = torch.linspace(0, 16000, mag.shape[0])
plt.plot(x_axis, mag)

plt.subplot(212)
plt.plot(x_axis, phase)
plt.xlabel('Freq [Hz]')

从图中可以注意到一些有趣的事情

  • 幅度图是对称的。x 轴的最后一个元素对应于采样频率 \(f_s\)，在本例中为 16kHz。由于这种对称性，只需绘制从 0 到 \(fs/2\) 的幅度即可。这个频率称为奈奎斯特频率 (Nyquist frequency)。

  • 相位图非常嘈杂。这也是预料之中的。相位出了名的难以解释和估计。

我们现在绘制从 0 到奈奎斯特频率的幅度

half_point = mag[0:].shape[0]//2
x_axis = torch.linspace(0, 8000, half_point)
plt.plot(x_axis, mag[0:half_point])
plt.xlabel('Frequency')

我们可以看到，语音信号的大部分能量集中在频谱的较低部分。实际上，许多重要的音素，如元音，其大部分能量都在频谱的这一部分。

此外，我们可以注意到幅度谱中的一些峰值。让我们放大以便更清晰地看到它们

plt.plot(mag[0:4000])
plt.xlabel('Frequency')

这些峰值对应于基频 (pitch)（即声带振动的频率）和共振峰 (formants)（对应于声道谐振频率）。

现在让我们尝试回到时域

signal_rec = torch.fft.ifft(fft, dim=0)
signal_rec = signal_rec # real part
signal_orig = signal

# Plots
plt.subplot(211)
plt.plot(signal_orig)

plt.subplot(212)
plt.plot(signal_rec)
plt.xlabel('Time')

print(signal_orig[0:10])
print(signal_rec[0:10])

从图中可以看出，信号可以在时域中重构。由于一些数值舍入误差，两个信号非常相似但不完全相同（参见前 10 个样本的打印输出）。

2. 短时傅里叶变换 (STFT)

语音是一种随时间演变的“动态”信号。因此，引入一种混合时频表示来显示语音的频率成分如何随时间演变是合理的。这种表示称为短时傅里叶变换。

STFT 按以下方式计算

  1. 使用重叠的滑动窗口（例如，汉明窗 (hamming)，汉宁窗 (hanning)，布莱克曼窗 (blackman)）将时域信号分割成多个块。

  2. 对每个小块计算 DFT

  3. 将所有 DFT 组合成一个单一的表示

现在我们计算一个语音信号的 STFT

from speechbrain.processing.features import STFT

signal = read_audio('/content/spk1_snt1.wav').unsqueeze(0) # [batch, time]

compute_STFT = STFT(sample_rate=16000, win_length=25, hop_length=10, n_fft=400) # 25 ms, 10 ms
signal_STFT = compute_STFT(signal)

print(signal.shape)
print(signal_STFT.shape)

  • STFT 表示的第一个维度是批量轴 (batch axis)（SpeechBrain 期望这样做，因为它旨在并行处理多个信号）。

  • 第三个是频率分辨率 (frequency resolution)。它对应于 FFT 点数 (\(n_{fft}\)) 的一半，因为如前所述，FFT 是对称的。

  • 最后一个维度汇集了 STFT 表示的实部和虚部。

与傅里叶变换类似，STFT 也有一个逆变换，称为逆短时傅里叶变换 (Inverse Short-Term Fourier Transform, ISTFT)。使用适当设计的窗口，我们可以完美重构原始信号。

from speechbrain.processing.features import ISTFT

compute_ISTFT = ISTFT(sample_rate=16000, win_length=25, hop_length=10)
signal_rec = compute_ISTFT(signal_STFT)
signal_rec = signal_rec.squeeze() # remove batch axis for plotting

# Plots
plt.subplot(211)
plt.plot(signal_orig)

plt.subplot(212)
plt.plot(signal_rec)
plt.xlabel('Time')

3. 频谱图

如前所述，傅里叶变换的幅度比相位包含更多信息。因此，我们可以取 STFT 表示的幅度，得到所谓的频谱图。频谱图是最流行的语音表示之一。

让我们看看频谱图是什么样的

spectrogram = signal_STFT.pow(2).sum(-1) # power spectrogram
spectrogram = spectrogram.squeeze(0).transpose(0,1)

spectrogram_log = torch.log(spectrogram) # for graphical convenience

plt.imshow(spectrogram_log.squeeze(0), cmap='hot', interpolation='nearest', origin='lower')
plt.xlabel('Time')
plt.ylabel('Frequency')
plt.show()

频谱图是一种二维表示，可以绘制成图像（黄色区域对应于幅度较高的时间-频率点）。从频谱图中，您可以看到频率成分如何随时间演变。例如，您可以清晰地分辨元音（其频率模式的特点是对应于基频和共振峰的多条线）和擦音（其特点是存在连续的高频成分）。通常，我们绘制的是功率谱图，它对应于 STFT 幅度的平方。

频谱图的时间分辨率和频率分辨率取决于用于计算 STFT 的窗口长度。

例如，如果我们增加窗口的长度，我们可以获得更高的频率分辨率（但时间分辨率较低）

signal = read_audio('/content/spk1_snt1.wav').unsqueeze(0) # [batch, time]

compute_STFT = STFT(sample_rate=16000, win_length=50, hop_length=10, n_fft=800)
signal_STFT = compute_STFT(signal)

spectrogram = signal_STFT.pow(2).sum(-1)
spectrogram = spectrogram.squeeze(0).transpose(0,1)
spectrogram = torch.log(spectrogram)

plt.imshow(spectrogram.squeeze(0), cmap='hot', interpolation='nearest', origin='lower')
plt.xlabel('Time')
plt.ylabel('Frequency')
plt.show()

反之亦然，我们可以获得更高的时间分辨率，但代价是频率分辨率降低

signal = read_audio('/content/spk1_snt1.wav').unsqueeze(0) # [batch, time]

compute_STFT = STFT(sample_rate=16000, win_length=5, hop_length=5, n_fft=800)
signal_STFT = compute_STFT(signal)

spectrogram = signal_STFT.pow(2).sum(-1)
spectrogram = spectrogram.squeeze(0).transpose(0,1)
spectrogram = torch.log(spectrogram)

plt.imshow(spectrogram.squeeze(0), cmap='hot', interpolation='nearest', origin='lower')
plt.xlabel('Time')
plt.ylabel('Frequency')
plt.show()

尽管频谱图信息量很大，但它是不可逆的。事实上，在计算它时，我们只使用了 STFT 的幅度，而没有使用相位。

频谱图是计算一些流行的语音特征的起点，例如滤波器组 (Filter Banks, FBANKs) 和梅尔频率倒谱系数 (Mel-Frequency Cepstral Coefficients, MFCCs)，这些特征是另一个教程的主题。

参考文献

[1] L. R. Rabiner, Ronald W. Schafer, “Digital Processing of Speech Signals”, Prentice-Hall, 1978

[2] S. K. Mitra Digital Signal Processing: A Computer-Based Approach 幻灯片

[3] https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/

[4] https://sites.northwestern.edu/elannesscohn/2019/07/30/developing-an-intuition-for-fourier-transforms/

引用 SpeechBrain

如果您在研究或商业中使用 SpeechBrain，请使用以下 BibTeX 条目引用它

@misc{speechbrainV1,
  title={Open-Source Conversational AI with {SpeechBrain} 1.0},
  author={Mirco Ravanelli and Titouan Parcollet and Adel Moumen and Sylvain de Langen and Cem Subakan and Peter Plantinga and Yingzhi Wang and Pooneh Mousavi and Luca Della Libera and Artem Ploujnikov and Francesco Paissan and Davide Borra and Salah Zaiem and Zeyu Zhao and Shucong Zhang and Georgios Karakasidis and Sung-Lin Yeh and Pierre Champion and Aku Rouhe and Rudolf Braun and Florian Mai and Juan Zuluaga-Gomez and Seyed Mahed Mousavi and Andreas Nautsch and Xuechen Liu and Sangeet Sagar and Jarod Duret and Salima Mdhaffar and Gaelle Laperriere and Mickael Rouvier and Renato De Mori and Yannick Esteve},
  year={2024},
  eprint={2407.00463},
  archivePrefix={arXiv},
  primaryClass={cs.LG},
  url={https://arxiv.org/abs/2407.00463},
}
@misc{speechbrain,
  title={{SpeechBrain}: A General-Purpose Speech Toolkit},
  author={Mirco Ravanelli and Titouan Parcollet and Peter Plantinga and Aku Rouhe and Samuele Cornell and Loren Lugosch and Cem Subakan and Nauman Dawalatabad and Abdelwahab Heba and Jianyuan Zhong and Ju-Chieh Chou and Sung-Lin Yeh and Szu-Wei Fu and Chien-Feng Liao and Elena Rastorgueva and François Grondin and William Aris and Hwidong Na and Yan Gao and Renato De Mori and Yoshua Bengio},
  year={2021},
  eprint={2106.04624},
  archivePrefix={arXiv},
  primaryClass={eess.AS},
  note={arXiv:2106.04624}
}