音频采样率转换
本文原始版本发布于 https://www.52pojie.cn/thread-1959816-1-1.html，此处进行了适当的精简，同时更新了一下代码（最新代码以 GitHub 仓库为准）。
前言
两年前，我研究了 WASAPI 播放音频的方法，详见https://www.cnblogs.com/PeaZomboss/p/17035785.html，挖了个坑，就是重采样的坑，具体就不细说了，前文后面提到了。因为现在学习了相关知识，所以现在来给它填上。
采样率转换（Sample Rate Conversion）通常也称作重采样（Resample），是对已有的数字信号采样进行处理的一种手段，在音频、图像处理等领域比较常见。本文将要讨论的内容，主要是在音频处理领域的常见采样率转换方法。
相信对于很多人来说，采样率转换似乎并不常见，但实际上你的手机、电脑每天都在运行相关的代码。在音频领域，只要手机或电脑在播放声音，那么系统极有可能会进行采样率转换的操作；在图像领域，对图片或视频的缩放就会对像素进行处理。
从音频角度来说，声卡驱动通常会提供多种不同的输出设置，如 44100Hz 16bit、48000Hz 24bit 等。但是系统通常只允许用户选择一种选项，可能是声卡不支持多采样率，也可能是系统不支持。因此当播放的音频文件采样率与系统设置的采样率不同时就涉及到采样率转换了。
在音乐行业，常见的 CD 音频标准是 44100Hz 16bit；而在影视行业，通常采用 48000Hz 的压缩格式（如 AAC）。比如使用音乐软件听歌，大部分情况下都是 44100Hz 的音频；而视频网站的视频内置的音频几乎都是 48000Hz 的。这种情况下，用户可以手动切换声卡设置，或者使用更专业的播放器，采用独占音频（一般是 WASAPI 或者是 ASIO）的方式使得声卡按照音频文件的采样率播放，从而达到无损的效果。而对于大部分普通板载声卡的驱动，可能只提供了一种采样率标准且甚至无法切换，因此必须使用采样率转换。
由于大部分声卡提供的输出采样率标准通常是 44100Hz 和 48000Hz，本文主要也着重于这二者之间的转换，但也会研究 96000Hz 的下采样问题（Hi-Res 通常要求至少 96k/24bit）。注意使用的算法其实是通用的，但是针对不同的采样率需要调整特定的参数才能达到最好的效果。
曾经的安卓系统就因为其内置的采样率转换而为人诟病，彼时的骁龙芯片更是让听歌体验变得极度糟糕。但是对于听歌来说，在音源已经达到 CD 标准的时候，好的播放和收听设备的影响显然是明显大于音源的。而且对于所谓的 Hi-Res，"母带级" 等音质，实际上更高的采样率对于听歌的来说毫无意义（也就是在音乐发行环节，参考此文：https://people.xiph.org/~xiphmont/demo/neil-young.html），但对于音乐的制作环节还是有价值的。
基本知识
首先要明确的是，本文所需的知识主要来自数字信号处理（Digital Signal Processing）领域，建议掌握一定的基础，以便于了解本文内容。推荐书籍有奥本海姆的《信号与系统》和《离散时间信号处理》，当然其他的书籍也是可以的。这些书籍的内容不必全都了解（我自己就没看太多），根据需要即可。
当然还需要一定代码能力，本文不涉及 MATLAB 相关内容，而是使用 Rust 语言。当然不会很复杂的，因为我自己也是属于 Rust 初学者，所以有很多地方都不了解。而且总体代码量不大，也可以轻松的用 C/C++ 等其他语言实现。此外，还需要一定的 Python 基础，因为分析时需要使用。
如果你看过上面说的基本书籍，通常会给出采样率转换相关的内容。一般来说，对于采样率从高到低的过程，称为减采样（或者下采样，downsampling）；反之则称为增采样（或者上采样，upsampling）。而具体的实施需要低通滤波器的配合，这样的过程则称为抽取（decimation）和内插（interpolation）。
上述过程都是针对整数倍数转换的，也就是说只能从 44100Hz 到 22050Hz，或者反过来。而要实现非整数倍的转换，则可以级联内插器和抽取器。比如从 44100Hz 到 48000Hz，可以先找出最大公因数 300，从而得出需要内插系数 L=160，抽取系数 M=147。
多采样率信号处理就是针对这个问题发展起来的。当然这样的实现也是相当复杂的，为了简单起见，文献 https://ccrma.stanford.edu/~jos/resample/resample.html 给出了一种基于 sinc 函数的带限插值的算法，网站 https://src.infinitewave.ca/ 提供了各种常见软件的各项技术指标。本文就按照这篇文献提出的方法，实现了具体的转换代码。
判断转换的质量
对于采样率转换来说，通常有很多的指标，通常是检测噪音的引入程度。细分一下包括对无关频率的影响、混叠和镜像频率的产生等。前面提供的网站给出了详细的技术指标，本文主要测试 1kHz 固定音调的正弦音频，以及频率在 \((0, F_S/2)Hz\) 的以二次函数变化的正弦音频（其中 \(Fs\) 表示采样率）。前者通常称为 beep，后者则称为 sweep。
为了能够准确地反映转换的质量，需要使用特定的工具来显示音频文件的频谱（Spectrum）以及 Spectrogram（可译为语谱图，也通常叫做时频图或者频谱图）。反正名字挺乱的不同的人可能说的不一样，建议以英文为主。其中 Spectrum 通常是对于某一时间段（比如 1024 个样本）或者全部样本的傅里叶变换的结果，对于 beep 音频来说很合适；而 Spectrogram 则是关于时间、频率和功率的图，功率用颜色表示，通过短时傅里叶变换实现，非常适合 sweep 音频。
除此之外，对于 sinc 插值方法，还需要使用冲激函数来获得滤波器的冲激响应，从而分析滤波器的性能和指标。冲激函数只需要在全是 0 的样本序列中出现一个值，无论多大就可以，因为最终都需要进行缩放。
生成测试文件
由于需要生成一个固定频率的正弦波和频率逐渐变化的正弦波，所以仅使用简单的 sin 函数显然无法满足要求。我们需要一个振荡器（Oscillator）来生成波形，这并不难实现。考虑正弦函数的特性，我们知道其相位在 \((0, 2\pi)\) 之间，幅度在 \((0, 1)\) 之间，频率则取决于周期 \(T\)。其生成公式为 \(x[n]=Asin(\omega n+\varphi)\)，其中 \(A\) 表示振幅，\(\omega=2\pi / T\)，\(\varphi\) 表示初始相位，\(n\) 是表示样本的整数。
由于数字音频的特点，实际的相位、频率和采样率相关。已知采样率 \(F_S\)，要求频率 \(F_A\)，显然周期 \(T=F_S/F_A\)，带入可得 \(\omega=2\pi F_A / F_S\) 。考虑到 sin 函数的周期性，我们可以在相位超过 \(2\pi\) 时减去 \(2\pi\)，一般我们使用 \(\tau\) 来表示 \(2\pi\) 。
因此代码可以这样来写：
use std::f64::consts::TAU; // 等于 2 * PI
struct Osc {
phase: f64, // 相位
omega: f64, // 角频率
freq: f64, // 目标频率
sample_rate: f64, // 采样率
}
impl Osc {
fn init(sample_rate: f64) -> Self {
Self {
phase: 0.0,
omega: 0.0,
freq: 0.0,
sample_rate,
}
}
fn set_freq(&mut self, freq: f64) {
self.freq = freq;
self.omega = TAU * freq / self.sample_rate;
}
fn next(&mut self) -> f64 {
let sample = self.phase.sin(); // 计算样本
self.phase += self.omega;
while self.phase >= TAU { // 用 while 是为了防止设定过高的频率
self.phase -= TAU;
}
sample
}
}
生成采样之后，还需要写入文件才行，这里我们使用了 Rust 库 hound，这个库提供了完整的 wav 文件的读写操作，支持整数和 32 位浮点数。为了减少量化带来的损失，同时方便操作，我们使用 32 位浮点数来保存文件。
下面是生成 beep 和 sweep 文件的代码：
use hound::{WavSpec, WavWriter};
fn gen_beep(sample_rate: u32) {
let filename = format!("beep_{}k.wav", sample_rate / 1000);
let spec = WavSpec {
channels: 1,
sample_rate,
bits_per_sample: 32,
sample_format: hound::SampleFormat::Float,
};
let mut writer = WavWriter::create(filename, spec).unwrap();
let mut osc = Osc::init(sample_rate as f64);
osc.set_freq(1000.0);
let sample_count = sample_rate * 5;
for _ in 0..sample_count {
let sample = osc.next() * 0.99;
writer.write_sample(sample as f32).unwrap();
}
writer.finalize().unwrap();
}
fn gen_sweep(sample_rate: u32) {
let filename = format!("sweep_{}k.wav", sample_rate / 1000);
let spec = WavSpec {
channels: 1,
sample_rate,
bits_per_sample: 32,
sample_format: hound::SampleFormat::Float,
};
let mut writer = WavWriter::create(filename, spec).unwrap();
let mut osc = Osc::init(sample_rate as f64);
let sample_count = sample_rate * 5;
let nyquist_freq = sample_rate as f64 / 2.0;
for i in 0..sample_count {
osc.set_freq(nyquist_freq * (i as f64 / sample_count as f64).powi(2));
let sample = osc.next() * 0.99;
writer.write_sample(sample as f32).unwrap();
}
writer.finalize().unwrap();
}
gen_beep 能生成指定采样率下一段长 5 秒的 1kHz 的正弦单声道音频，听起来就像嘟~的声音，为了避免将来转换时可能出现的振幅过大的问题，最后还将振幅乘了 0.99 以减少失真的可能性，这样实际最大音量约是 -0.1dBFS，因为通常我们假定 ±1 是一般系统能处理的最大值。gen_sweep 则是能生成指定采样率下同样规格的文件，区别是正弦的频率是不断变化的，从 0Hz 到奈奎斯特频率，振幅也同样进行了限制。
检测音频并绘图
由于生成的音频文件仅通过播放是几乎无法听出具体的转换结果是否符合预期，因此使用专业的工具进行分析是很重要的。一般来说这个时候应该请出 MATLAB 了，不过相比这个大家伙，我觉得还是 Python 更加灵活。我们使用 numpy、scipy 和 matplotlib.pyplot 来进行分析和绘制。
对于前面提到的 Spectrum，我们只需要对文件进行 FFT 就可以得到了；而 Spectrogram 则是需要使用 STFT 进行分析。除此之外还专门在 Spectrum 的基础上进行了一定的修改，以展示滤波器的冲激响应。具体代码如下：
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from scipy.signal import ShortTimeFFT
from scipy.signal.windows import kaiser
def _show_spectrum(fs, data, name, impulse: None | str = None):
passband = impulse == 'passband'
N = len(data)
half_N = N // 2
fft_data = abs(np.fft.fft(data))
fft_data = fft_data / half_N if impulse is None else fft_data / max(fft_data)
fft_dBFS = 20 * np.log10(fft_data)
freqs = np.fft.fftfreq(N, 1/fs)
plt.figure(figsize=(6, 4))
xticks = np.arange(0, fs // 2 + 1, 2000)
xticklabels = [f'{int(tick/1000)}' for tick in xticks]
ymin, ymax, ystep = (-3, 1, 0.5) if passband else (-200, 10, 20)
ax = plt.gca()
ax.set(xlabel='Frequency in kHz', ylabel='Magnitude in dBFS',
xlim=(0, fs//2), ylim=(ymin, ymax),
xticks=xticks, yticks=np.arange(ymin, ymax, ystep),
xticklabels=xticklabels, facecolor='black')
ax.plot(freqs[:half_N], fft_dBFS[:half_N], color='white')
ax.grid()
prefix = 'Passband of ' if passband else 'Spectrum of '
plt.title(prefix + name)
plt.show()
# 显示wav文件频谱
def show_wav_spectrum(filename):
fs, data = wavfile.read(filename)
_show_spectrum(fs, data, filename)
# 显示wav冲激响应频谱
def show_wav_impulse(filename, passband=False):
fs, data = wavfile.read(filename)
_show_spectrum(fs, data, filename, impulse='passband' if passband else '')
# 显示原始f64数据冲激响应频谱
def show_raw_impulse(filename, fs, passband=False):
data = np.fromfile(filename, np.float64)
_show_spectrum(fs, data, filename, impulse='passband' if passband else '')
# 显示wav文件spectrogram
def show_wav_spectrogram(filename):
fs, data = wavfile.read(filename)
N = len(data)
window_size = 2048
hop = window_size // 2
win = kaiser(window_size, 20)
SFT = ShortTimeFFT(win, hop, fs, scale_to='magnitude')
Sx = SFT.stft(data)
fig = plt.figure(figsize=(6, 4))
ax = plt.gca()
yticks = np.arange(0, fs // 2 + 1, 2000)
yticklabels = [f'{int(tick/1000)}' for tick in yticks]
ax.set(xlabel='Time in seconds', ylabel='Frequency in kHz',
yticks=yticks, yticklabels=yticklabels)
im = ax.imshow(20*np.log10(abs(Sx)), origin='lower', aspect='auto',
extent=SFT.extent(N), cmap='inferno', # 我觉得不错的配色
vmin=-200, vmax=1, # 最大最小值
interpolation='sinc') # sinc插值比较准确
fig.colorbar(im, label="Magnitude in dBFS", ticks=np.arange(-200,1,20))
plt.title(f'Spectrogram of {filename}')
plt.show()
这里有很多细节，包括对数据进行的各种处理，图像的坐标轴的处理，图像颜色的选择等。而且有关短时傅里叶变换的 ShortTimeFFT 在网上几乎没有资料，只能查官方文档。为了得到一张效果较为完美的图像，我也是查阅大量资料，前前后后多次修改代码。
看一下 44100Hz 采样率的 beep 文件的 Spectrum：
以及 48000Hz 采样率的：
能看出来因为 1kHz 和采样率的倍数关系，噪音的引入还是有所区别的。但是二者的噪音也都是在 -160dBFS 以下，属于几乎不可听的状态。
当然对于 44100Hz 采样率 sweep 的 Spectrogram 也放在这里供参考，48000Hz 的就不放了，因为图都是差不多的：
这些图是不是看上去都挺不错的，而经过采样率转换之后，它们还能保持这样的状态吗？接着往后看。
代码框架设计
实际上这一块是后来搞的，因为一开始根本没考虑这些，可以去代码仓库看看提交历史。不过确实狠狠学了不少 Rust 知识，只能说真复杂啊，和编译器斗智斗勇的。
首先是核心 Trait：
trait NextSample {
fn next_sample<F>(&mut self, f: &mut F) -> Option<f64>
where
F: FnMut() -> Option<f64>;
}
所有的和插值相关的代码都应该实现这个 Trait，这样为后续的包装提供了好处。注意这个 next_sample 是一个泛型函数，因此这个 Trait 不支持动态派发，这也是我单独拆分这个 Trait 的原因，因为后面还有一个 Convert 的 Trait，可以动态派发：
pub trait Convert {
fn process<I>(&mut self, iter: I) -> ConvertIter<'_, I, Self>
where
I: Iterator<Item = f64>,
Self: Sized;
}
impl<T: NextSample> Convert for T {
#[inline]
fn process<I>(&mut self, iter: I) -> ConvertIter<'_, I, Self> {
ConvertIter::new(iter, self)
}
}
process 接受一个迭代器参数，然后返回一个 ConvertIter 结构，该结构实现了一个迭代器，返回转换后的输出，而直接返回 impl Iterator<Item = f64> 也可以，但就是不支持动态派发了。由于使用了 next_sample 而不是直接使用迭代器的输入，所以当输入迭代器出现 None 时，再次输入有效数据还能接着输出剩余的内容，这样就可以避免数据无法通过一个连续的迭代器提供的情况。
ConvertIter 结构和迭代器实现如下：
pub struct ConvertIter<'a, I, C> {
iter: I,
cvtr: &'a mut C,
}
impl<'a, I, C> ConvertIter<'a, I, C> {
#[inline]
pub fn new(iter: I, cvtr: &'a mut C) -> Self {
Self { iter, cvtr }
}
}
impl<I, C> Iterator for ConvertIter<'_, I, C>
where
I: Iterator<Item = f64>,
C: NextSample,
{
type Item = f64;
#[inline]
fn next(&mut self) -> Option<Self::Item> {
self.cvtr.next_sample(&mut || self.iter.next())
}
}
看上去好像还整上生命周期了，好像很麻烦的样子，其实一点也不。这是为了告诉编译器这个迭代器的存活时间最多和实现了 NextSample 的对象相同。实际上就是简单调用了 next_sample 这个函数而已。因此我们接下来的关键就是如何为不同的插值方法都实现这个函数。而且看起来这些复杂的包装代码，其实不会造成任何损失，因为这些都是静态的，也就是编译时确定的，这就是所谓的零成本抽象。
因此对于任意一个实现了 Convert Trait 的转换器，都可以使用类似这样的代码实现转换：
let samples = vec![1.0, 2.0];
let mut cvtr = get_converter();
for s in cvtr.process(samples.into_iter()) {
println!("sample = {s}");
}
线性插值
首先出场的是线性插值，这应该是最容易想到的方法了，也是最容易实现的方法。
这个很好理解，比如从采样率 5 降低到 4，比率是 0.8，步长就变成了 1.25。那么对于时长为一秒的数据 \(x[0]=1, x[1]=3, x[2]=5, x[3]=7, x[4]=9\) ，我们需要 \(y[0]=x(0), y[1]=x(1.25), y[2]=x(2.5), y[3]=x(3.75)\) ，但是 \(x[n]\) 的索引不能是小数，所以就用线性插值公式计算。比如 \(y[1]=x(1.25)=x[1]+0.25 \times (x[2]-x[1])=3+0.25 \times (5-3)=3.5\) ，同样的可以算出每一个 \(y[m]\) 出来。
代码非常简单：
enum State {
First, // 第一次状态
Normal, // 正常状态
Suspend, // 挂起状态，也就是输入存在None
}
pub struct Converter {
step: f64, // 插值位置步长
pos: f64, // 当前插值位置，0到1
last_in: [f64; 2], // 插值点前后的两个样本
state: State, // 当前状态
}
fn next_sample<F>(&mut self, f: &mut F) -> Option<f64>
where
F: FnMut() -> Option<f64>,
{
loop {
match self.state {
State::First => {
if let Some(s) = f() { // 尝试读取
self.last_in[1] = s; // 缓存
self.pos = 1.0; // 此时置1.0是为了能在正常状态时接着读取，因为至少两个样本才能插值
self.state = State::Normal; // 切换到正常状态
} else {
return None; // 此时仍然是初始状态
}
}
State::Normal => {
while self.pos >= 1.0 {
self.pos -= 1.0; // 移动一个样本
self.last_in[0] = self.last_in[1]; // 移进缓存
if let Some(s) = f() {
self.last_in[1] = s; // 读入新的值
} else {
self.state = State::Suspend; // 切换到挂起状态
return None; // 返回空
}
}
// 执行插值的代码
let interp = self.last_in[0] + (self.last_in[1] - self.last_in[0]) * self.pos;
self.pos += self.step; // 增加当前的位置
return Some(interp);
}
State::Suspend => {
if let Some(s) = f() {
self.last_in[1] = s;
self.state = State::Normal; // 读取成功切换到正常状态
} else {
return None; // 否则继续返回None
}
}
}
}
}
不过显而易见的，只根据前后两个样本进行插值计算，显然有较大的误差，尤其是因为声音信号的波形通常都是是由一系列不同频率的周期信号组成，而两个样本不足以表达这些信号的周期。但是线性插值并非一无是处，因为其速度快，而且对于低频信号的影响并没有那么大，所以对于一些性能较差的设备可能还是存在一定的价值。同时，对于线性信号（三角波、锯齿波、方波等）来说，也是有不错的效果。
这是使用线性插值算法将 44.1k 的 beep 文件转换到 48k 的结果：
对比前面的图，是不是差距还挺大的？其实仔细看可以发现，最大的噪音在 -60dB 以下，剩下的基本都在 -80dB 以下，其余多次混叠的噪音则都在 -100dB 以下，实际上造成的影响并没有那么大。另一张 48k 下采样到 44.1k 的图就不放了，基本上是差不多的。说实话，这样的噪音我反正是听不出来的，想听到得多大音量，怕不是耳朵都要聋了。
而 sweep 就有意思了，因为 sweep 有两个特点，一个是频率不断变化，另一个是频率高点达到奈奎斯特频率，这会导致出现强烈的频率镜像现象和频率混叠现象。
对照原图，可以看到在低频的时候噪音的强度还不突出，而到了高频部分噪音明显变多。看起来是不是很乱？好像每次到了边界就会反弹，而且是不断的反弹，直到信号衰减到检测不到。这其实是由多次的混叠和镜像共同造成的。这回一下就能听出来了，尤其是当声音来到高频的时候，那些低频的混叠就会非常明显。
所以一般为了降低干扰，线性插值还会加入一个 IIR 滤波器，比如经典的巴特沃斯（Butterworth）滤波器。不过 IIR 滤波器会导致非线性相位，在有的时候并不是最好的选择。由于线性插值并不是本文的核心内容，所以就没有深入分析，只是简单做一个说明，因为其远不如接下来马上要探讨的 sinc 函数插值方法。一般来说，线性插值仅适用于需要快速处理，对质量要求不高的低端嵌入式设备场景。
sinc 插值
采样定理告诉我们，对于带限信号，使用低通滤波器就可以重建连续时间的函数。而带限条件下的 sinc 插值是完美的信号重建方法，因为 sinc 函数本身就是一个理想的低通滤波器。
sinc 函数，定义为 \(sinc(x)=\frac{sin(\pi x)}{\pi x}\) ，其中当输入为 0 的时候，结果是 1 。这个函数为什么长这样，这主要还是傅里叶变换的功劳，一般教材都会有推导过程，了解即可。因为 sinc 函数只有在正负无穷的时候才会达到 0，所以需要近似，最简单的方法就是把 sinc 函数截断，而实际使用中通常会加窗处理。这样就是一个 FIR 滤波器了，但是作为插值滤波器，和一般的滤波器的区别在于不能提前计算系数（可以近似计算）。
sinc 插值的方法其实很简单，使用这个公式可以实现： \(\hat{x}(t)=\sum_{n=-\infty}^{\infty}x(nT_s)h_s(t-nT_s)\) ，其中，\(h_s(t)=sinc(F_st)\) 。仔细看可以发现这个公式其实是对整个时域信号做了一个卷积操作，我们知道时域的卷积对应了频域的乘积，所以相当于进行了滤波。但是这是理想的情况，实际上我们不可能有无穷多的点，而且也不可能用那么多的点只是为了计算一个值。所以按照前面说的，使用窗函数法进行近似就是具体的方法了。
其实和线性插值几乎一样，我们只需要先找到要插值的位置，然后找到这个地方前后一共 M 个点（M 表示阶数，N 表示滤波器长度，M=N-1，也有地方直接用 N 表示阶数，用 L 表示长度），再一顿卷积就可以了。伪代码就像这样：
double sum = 0.0;
for (int k = -M; k <= M; k++) {
sum += x[k] * sinc(pos-k) * kaiser(pos-k, M, beta);
}
return sum;
其中 pos 表示插值的位置，这和线性插值是一样的。需要注意的是，Kaiser 窗提供了一个 beta 参数，可以控制对旁瓣的抑制。但是 beta 不是越大越好，这是和阶数有关系的，一般来说，阶数高了再使用更大的 beta，具体的组合需要实测才能得出。当然实际使用的时候，还需要加一个 cutoff 参数，用来控制滤波器的截止频率。
由于 Kaiser 窗函数计算复杂，如果像这样每一个点就要调用 M 次，显然性能是不够的。考虑到插值的位置是一个 0 到 1 的小数，我们完全可以量化这个插值系数，比如提前将其划分为 Q 份，并根据实际的系数对量化后的系数再进行一次线性插值（没想到吧，线性插值表示我还在呢）。这样我们只需要保存 N*Q 个点的滤波器系数表。通常来说这个 Q 值的选取取决于精度和存储空间，比如对于嵌入式设备来说，查找表大小受 Flash 大小限制。实际上，由于 sinc 函数和 Kaiser 窗函数都有偶对称的特性，所以只需要保留其一半的系数就可以了，这样又能节省一半的空间。
sinc 函数和 Kaiser 窗函数的代码如下：
use std::f64::consts::PI;
fn sinc_c(x: f64, cutoff: f64) -> f64 {
if x != 0.0 {
(PI * x * cutoff).sin() / (PI * x)
} else {
cutoff
}
}
fn bessel_i0(x: f64) -> f64 {
let mut y = 1.0;
let mut t = 1.0;
for k in 1..32 {
t *= (x / (2.0 * k as f64)).powi(2);
y += t;
if t < 1e-10 {
break;
}
}
y
}
fn kaiser(x: f64, order: u32, beta: f64) -> f64 {
let half = order as f64 * 0.5;
if (x < -half) || (x > half) {
return 0.0;
}
bessel_i0(beta * (1.0 - (x / half).powi(2)).sqrt()) / bessel_i0(beta)
}
当然实际上我们并没有直接使用 sinc 的定义，而是加入了 cutoff 参数，这样子主要是为了方便之后的使用。而零阶贝塞尔函数是 Kaiser 窗需要的，这里进行了一点点优化，展开了阶乘并转换成了迭代式，同时限制了最大迭代次数提供性能（此时误差很小）。而这个 kaiser 函数实际上并没有使用，而是直接用了贝塞尔函数，目的还是为了减少计算量（对于同一个 beta，分母的值是固定的）。
为了生成一个系数表，我们需要以下几个参数：
-
ratio：采样率转换比率。目标采样率除以原始采样率，除了用于步长的计算，还用于滤波器理想截止频率的计算。比如 44100Hz 转换到 48000Hz 时，为了抑制镜像频率，需要从 22050Hz 频率处进行过滤，由于原始文件采样率是 44100Hz，所以理想截止频率相对于奈奎斯特频率的比率就是 1；而从 48000Hz 到 44100Hz 的时候，为了抑制混叠，理想截止频率也是 22050Hz，其相对于奈奎斯特频率的比率是 22050/24000，而这恰好等于了 ratio。
-
order：滤波器阶数。阶数越大滤波过渡带越窄，滤波效果也越好，但也会增加计算量。一般来说，阶数为偶数的滤波器就是 Ⅰ 型 FIR 滤波器，奇数阶数的是 Ⅱ 型 FIR 滤波器，对于低通滤波器来说二者效果一致。滤波器延迟的样本数计算方法是 M / 2 * ratio 。
-
quan：量化数量。数量越多精度越高，但也占用更多存储空间。一般来说，量化数量不能过低，否则会影响插值的准确性，从而产生较多的噪音。具体取值和样本的精度有关，在前面提到的文献中有所论证。
-
kaiser_beta：Kaiser 窗函数 beta 值。越大对旁瓣的衰减更多，但在相同滤波器阶数下对过渡带的要求也会更高，因此需要和滤波器的阶数配合使用，不能盲目选择更大的值。一般可以通过经验公式计算目标衰减下的 beta 值。
-
cutoff：滤波器截止频率系数。实际的截止频率通常低于理想的情况，用来抑制镜像和混叠，具体取值和阶数有关。阶数越高，过渡带越窄，则可以越接近理想的截止频率。理想情况下，上采样时取 1，下采样时取转换比率。由于 FIR 滤波器的特性，截止频率是通带和阻带频率的中心点，也就是衰减 6dB 的点。
使用如下代码就可以根据这些参数计算插值滤波器的系数表了：
fn generate_filter_table(quan: u32, order: u32, beta: f64, cutoff: f64) -> Vec<f64> {
let len = order * quan / 2;
let i0_beta = bessel_i0(beta);
let half_order = order as f64 * 0.5;
let mut filter = Vec::with_capacity(len as usize + 1);
for i in 0..len {
let pos = i as f64 / quan as f64;
let i0_1 = bessel_i0(beta * (1.0 - (pos / half_order).powi(2)).sqrt());
let coef = sinc_c(pos, cutoff) * (i0_1 / i0_beta);
filter.push(coef);
}
filter.push(0.0);
filter
}
然后就可以利用这个表进行插值计算了：
fn interpolate(&self) -> f64 {
let coef = self.pos;
let mut interp = 0.0;
let pos_max = self.filter.len() - 1;
for (i, s) in self.buf.iter().enumerate() {
let index = i as f64 - self.half_order;
// 计算在插值滤波器系数表的位置
let pos = (coef - index).abs() * self.quan;
let pos_n = pos as usize;
// 如果位置超出最大范围就不计算了，因为此时系数是0
if pos_n < pos_max {
let h1 = self.filter[pos_n];
let h2 = self.filter[pos_n + 1];
// 对系数进行一次线性插值
let h = h1 + (h2 - h1) * (pos - pos_n as f64);
interp += s * h;
}
}
interp
}
这里代码不是很完整，只是贴出了主要的插值方法，总体的逻辑和线性插值很相似，甚至还更简单一点，完整代码会在初始化参数计算后面给出。
对于初始化参数的计算，我们可以按照教材上使用窗函数法设计 FIR 滤波器的部分内容。
首先需要确定阻带衰减水平 A，然后据此可以计算 Kaiser 窗的 β 值。通常来说我们选择的 A 都是大于 50 的，所以基本上就是按照 \(\beta=0.1102\times(A-8.7)\) 计算的。
fn calc_kaiser_beta(atten: f64) -> f64 {
if atten > 50.0 {
0.1102 * (atten - 8.7)
} else if atten >= 21.0 {
0.5842 * (atten - 21.0).powf(0.4) + 0.07886 * (atten - 21.0)
} else {
0.0
}
}
然后根据公式 \(M=\frac{A-8}{2.285\varDelta\omega}\)，可以根据阶数 M 或者过渡带宽度 \(\varDelta\omega\) 来计算。这里有一个重要的细节，就是对于采样率转换的情况，过渡带的宽度是需要根据上采用还是下采样进行缩放的，就像之前提到的截止频率 \(\omega_c\) 在下采样需要乘转换系数一样。比如从 44100 到 96000 的上采样，我们以 20000Hz 作为通带，22050Hz 作为阻带，那么截止频率为 21025Hz，对应到角频率就是通带 \(\omega_p\approx0.907\pi\)，阻带 \(\omega_s=\pi\)，过渡带宽度 \(\varDelta\omega\approx0.093\pi\)，截止频率 \(\omega_c\approx0.9535\pi\) 。但是如果其他条件不变，变成了下采样，此时最大频率来到了 48000Hz，所以通带 \(\omega_p\approx0.4167\pi\)，阻带 \(\omega_s=0.459375\pi\)，过渡带宽度 \(\varDelta\omega\approx0.0427\pi\)，截止频率 \(\omega_c\approx0.438\pi\)。此时发现过渡带宽度为原来的 \(\frac{44100}{96000}\)，刚好就是转换比率 R 的倒数。因此上述公式修改为 \(M=\frac{A-8}{2.285\varDelta\omega\cdot\min(R,1)}\) 即可计算我们需要的结果了。代码中不要忘了 \(\pi\) 值，这样得到的就是归一化的结果了。
fn calc_trans_width(ratio: f64, atten: f64, order: u32) -> f64 {
(atten - 8.0) / (2.285 * order as f64 * PI * ratio.min(1.0))
}
fn calc_order(ratio: f64, atten: f64, trans_width: f64) -> u32 {
f64::ceil((atten - 8.0) / (2.285 * trans_width * PI * ratio.min(1.0))) as u32
}
之后就可以通过输入其中一个参数来计算另一个参数了。利用这些方法，我们可以根据阻带衰减、阶数或过渡带宽度来计算所需的全部参数了。
// 如果输入的是过渡带宽度
let order = calc_order(ratio, atten, trans_width);
// 如果输入的是阶数
let trans_width = calc_trans_width(ratio, atten, order);
// 之后就是一样的计算了
let kaiser_beta = calc_kaiser_beta(atten);
let cutoff = ratio.min(1.0) * (1.0 - 0.5 * trans_width);
这里给出完整的代码：
use std::collections::VecDeque;
use std::f64::consts::PI;
use std::sync::Arc;
use super::NextSample;
#[inline]
fn sinc_c(x: f64, cutoff: f64) -> f64 {
if x != 0.0 {
(PI * x * cutoff).sin() / (PI * x)
} else {
cutoff
}
}
#[inline]
fn bessel_i0(x: f64) -> f64 {
let mut y = 1.0;
let mut t = 1.0;
for k in 1..32 {
t *= (x / (2.0 * k as f64)).powi(2);
y += t;
if t < 1e-10 {
break;
}
}
y
}
#[inline]
#[allow(dead_code)]
fn kaiser(x: f64, order: u32, beta: f64) -> f64 {
let half = order as f64 * 0.5;
if (x < -half) || (x > half) {
return 0.0;
}
bessel_i0(beta * (1.0 - (x / half).powi(2)).sqrt()) / bessel_i0(beta)
}
#[inline]
fn generate_filter_table(quan: u32, order: u32, beta: f64, cutoff: f64) -> Vec<f64> {
let len = order * quan / 2;
let i0_beta = bessel_i0(beta);
let half_order = order as f64 * 0.5;
let mut filter = Vec::with_capacity(len as usize + 1);
for i in 0..len {
let pos = i as f64 / quan as f64;
let i0_1 = bessel_i0(beta * (1.0 - (pos / half_order).powi(2)).sqrt());
let coef = sinc_c(pos, cutoff) * (i0_1 / i0_beta);
filter.push(coef);
}
filter.push(0.0);
filter
}
#[inline]
fn calc_kaiser_beta(atten: f64) -> f64 {
if atten > 50.0 {
0.1102 * (atten - 8.7)
} else if atten >= 21.0 {
0.5842 * (atten - 21.0).powf(0.4) + 0.07886 * (atten - 21.0)
} else {
0.0
}
}
#[inline]
fn calc_trans_width(ratio: f64, atten: f64, order: u32) -> f64 {
(atten - 8.0) / (2.285 * order as f64 * PI * ratio.min(1.0))
}
#[inline]
fn calc_order(ratio: f64, atten: f64, trans_width: f64) -> u32 {
f64::ceil((atten - 8.0) / (2.285 * trans_width * PI * ratio.min(1.0))) as u32
}
enum State {
Normal,
Suspend,
}
pub struct Converter {
step: f64,
pos: f64,
half_order: f64,
quan: f64,
filter: Arc<Vec<f64>>,
buf: VecDeque<f64>,
state: State,
}
impl Converter {
#[inline]
pub fn new(step: f64, order: u32, quan: u32, filter: Arc<Vec<f64>>) -> Self {
let taps = (order + 1) as usize;
let mut buf = VecDeque::with_capacity(taps);
buf.extend(std::iter::repeat(0.0).take(taps));
Self {
step,
pos: 0.0,
half_order: 0.5 * order as f64,
quan: quan as f64,
filter,
buf,
state: State::Normal,
}
}
#[inline]
fn interpolate(&self) -> f64 {
let coef = self.pos;
let mut interp = 0.0;
let pos_max = self.filter.len() - 1;
for (i, s) in self.buf.iter().enumerate() {
let index = i as f64 - self.half_order;
let pos = (coef - index).abs() * self.quan;
let pos_n = pos as usize;
if pos_n < pos_max {
let h1 = self.filter[pos_n];
let h2 = self.filter[pos_n + 1];
let h = h1 + (h2 - h1) * (pos - pos_n as f64);
interp += s * h;
}
}
interp
}
}
impl NextSample for Converter {
#[inline]
fn next_sample<F>(&mut self, f: &mut F) -> Option<f64>
where
F: FnMut() -> Option<f64>,
{
loop {
match self.state {
State::Normal => {
while self.pos >= 1.0 {
self.pos -= 1.0;
if let Some(s) = f() {
self.buf.pop_front();
self.buf.push_back(s);
} else {
self.state = State::Suspend;
return None;
}
}
self.pos += self.step;
let interp = self.interpolate();
return Some(interp);
}
State::Suspend => {
if let Some(s) = f() {
self.buf.pop_front();
self.buf.push_back(s);
self.state = State::Normal;
} else {
return None;
}
}
}
}
}
}
pub struct Manager {
ratio: f64,
order: u32,
quan: u32,
latency: usize,
filter: Arc<Vec<f64>>,
}
impl Manager {
pub fn with_raw(ratio: f64, quan: u32, order: u32, kaiser_beta: f64, cutoff: f64) -> Self {
let filter = generate_filter_table(quan, order, kaiser_beta, cutoff);
let latency = (ratio * order as f64 * 0.5).round() as usize;
Self {
ratio,
order,
quan,
latency,
filter: Arc::new(filter),
}
}
#[inline]
pub fn new(ratio: f64, atten: f64, quan: u32, trans_width: f64) -> Self {
let kaiser_beta = calc_kaiser_beta(atten);
let order = calc_order(ratio, atten, trans_width);
let cutoff = ratio.min(1.0) * (1.0 - 0.5 * trans_width);
Self::with_raw(ratio, quan, order, kaiser_beta, cutoff)
}
#[inline]
pub fn with_order(ratio: f64, atten: f64, quan: u32, order: u32) -> Self {
let kaiser_beta = calc_kaiser_beta(atten);
let trans_width = calc_trans_width(ratio, atten, order);
let cutoff = ratio.min(1.0) * (1.0 - 0.5 * trans_width);
Self::with_raw(ratio, quan, order, kaiser_beta, cutoff)
}
#[inline]
pub fn converter(&self) -> Converter {
Converter::new(
self.ratio.recip(),
self.order,
self.quan,
self.filter.clone(),
)
}
#[inline]
pub fn latency(&self) -> usize {
self.latency
}
#[inline]
pub fn order(&self) -> u32 {
self.order
}
}
转换质量分析
这里给出 A=100，Q=128，M=128 参数时 44k 到 48k 的具体情况：
可以看到基本满足要求，通带也刚好卡到了 20kHz 左右。如果以同样的参数，作用于 48k 到 44k 的下采样会怎么样呢？
其他两张图差不太多就不放了，关键就是这张通带图，可以看到通带已经不足 20kHz 了，唯一的解决办法就是加大阶数。但是加大阶数会增加计算量吗？我们接着讨论这个问题。
还是一样的配置，唯独把过渡带宽度设置为 \(\frac{2050}{22050}\)，这样就可以确保通带一定是满足 20kHz 的要求的。具体结果就不放了，因为衰减还是 100dB，所以效果其实和之前差不多，唯一的区别就是通带的频率都保证一样了。此时对于 44k 到 48k，M=138，而对于 48k 到 44k，M=151，由于计算量取决于输出点的数量乘阶数，所以对于上采样，每秒的计算量是 \(Calc=48000\times138=6624000\)，而下采样时计算量是 \(Calc=44100\times151=6659100\)，可以看到计算量并没有增加多少，不过滤波器系数的存储空间增加是必然的。
如果保持其他参数不变，将 A 增大到 120，就会出现量化精度不足的情况。但是光看通带图并不能看出问题来，需要配合扫频图才能得出。
看起来没什么问题，但是对比下图就知道问题挺大：
而只有将量化数量 Q 增大到 512 才能达到预期要求：
更大的衰减
如果需要测定更大的衰减下的参数选择，则需要使用 f64 浮点数来进行，这是因为 f32 的 wav 文件对于 144dB 以上的衰减是存在误差的，因为 f32 浮点数的有效精度只有 24 位。由于 hound 库不支持 f64（尽管 SciPy 支持），因此只能直接将 f64 的原始数据写入文件进行测试了。
这里说一下大概的情况：A=140 大概需要 Q=2048，A=150 则大约是 Q=4096，A=160 需要 Q=8192（这里的量化都只考虑 2 的幂次方）。当然也不是不可以降低一些要求，比如强行 A=160 和 Q=4096 也能有大概 155dB 的衰减，同理 A=150 和 Q=2048 也能达到 145dB 左右。
下图就是 A=150 和 Q=2048 下 96k 到 44.1k 的下采样结果（通带 20kHz，此时滤波器阶数达到了 464）：
同样条件下将 A 下调到 120，Q 下调到 512，则结果如图（此时依然需要 366 阶）：
下表给出了几种常见下采样的滤波器阶数（通带均为 20kHz）：
| 96->44.1 | 96->48 | 192->44.1 | 192->48 | |
|---|---|---|---|---|
| A=120 | M=366 | M=188 | M=731 | M=375 |
| A=150 | M=464 | M=238 | M=927 | M=475 |
性能测试
从前面的表格可以看到计算压力应该不小，于是我们尝试进行一个性能测试。由于 Rust 官方的 benchmark 尚且不够完善且需要 nightly 版本编译器，因此我们选择第三方性能测试框架。一开始我是选择了 Criterion.rs ，后来改用了更轻量级的 divan。
具体测试内容包括 A=120 和 A=150 两种衰减参数下 44.1k，48k 和 96k 互相转换的耗时。这个时间分成两部分：初始化滤波器系数耗时和处理 10 毫秒数据的耗时。
下表是使用骁龙 8cx Gen3（大核性能相当于骁龙 888 超频版本）测试得出的数据，采用的都是平均耗时：
| Init a120 | Init a150 | Proc a120 10ms | Proc a150 10ms | |
|---|---|---|---|---|
| 44k to 48k | 2.506 ms | 13.71 ms | 292.2 µs | 500.9 µs |
| 44k to 96k | 2.337 ms | 14.03 ms | 543.3 µs | 973.4 µs |
| 48k to 44k | 2.714 ms | 14.44 ms | 279.5 µs | 432.8 µs |
| 48k to 96k | 1.322 ms | 7.648 ms | 344.6 µs | 412.3 µs |
| 96k to 44k | 5.071 ms | 29.77 ms | 534.6 µs | 1.016 ms |
| 96k to 48k | 2.724 ms | 15.6 ms | 302.5 µs | 387 µs |
这个结果比手机端的骁龙 865 强了不少，相比 x86 平台的 AMD 锐龙 4800H 也更好，部分结果好于 GitHub Action 的机器（https://github.com/DrPeaboss/simple_src/actions/runs/10862155141/job/30144857691 ），但逊色于去年 Intel 推出的 Ultra 7 155H。基本上这些 CPU 性能都绰绰有余，可以同时跑好几路的实时转换。
注意到 divan 还会提供最好、最坏、中位数等数据，我们除了关心平均耗时，还需要注意最坏情况。如果处理 10 毫秒数据的耗时超过了 10 毫秒，则会出现爆音的情况，这个通常是线程无法获得足够的 CPU 资源引起的。因此通常需要提高线程优先级保证系统能调度足够的 CPU 资源。
当然一般在实时音频处理的时候不建议采用计算性能要求过高的参数，通常来说 A=120 就可以了，如果还是不行可以进一步降低到 A=100，这个时候绝大部分多年前的 CPU 的性能应该都吃得消了。
结语
有关此次采样率转换的研究，我全部写在这里了。从知识的学习，到代码的编写，看了不少书，也重写了几次代码，学会了很多东西，又花了很长时间写完本文，时间跨度从巴黎奥运会开始前到黑神话悟空发售后，历时一月有余。本文虽然到这结束了，但是知识的学习远没有结束。有关采样率转换还有很多我不知道的内容，正等待进一步的探索。同时本文肯定存在不少疏漏之处，还希望高手能够指出。
本文全部代码在此 GitHub 仓库可以找到：https://github.com/DrPeaboss/simple_src ，如果无法访问，也可以在此 Gitee 仓库找到：https://gitee.com/peazomboss/simple_src 。目前已经上传第一个版本到 crates.io，链接：https://crates.io/crates/simple_src ，使用命令 cargo add simple_src 即可直接添加到项目。
当然，如果需要成熟的采样率转换方案，那么 r8brain、libsamplerate、libsoxr 都是不错的选择。或者，通过 GitHub，一起完善这个简单的采样率转换器吧！
参考
- 奥本海姆《信号与系统》《离散时间信号处理》
- https://ccrma.stanford.edu/~jos/resample/resample.html 重采样算法
- https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.ShortTimeFFT.html SciPy 短时傅里叶变换文档
- https://matplotlib.org/stable/api/pyplot_summary.html#
matplotlib.pyplot库的 API 参考 - https://src.infinitewave.ca/ 各项采样率转换测试指标和各种软件数据
- https://kaisery.github.io/trpl-zh-cn/ Rust Book 中文翻译版
- https://doc.rust-lang.org/cargo/index.html 官方 Cargo 教程
- https://people.xiph.org/~xiphmont/demo/neil-young.html 为什么 192k/24bit 毫无意义，对所谓“无损”音乐有兴趣推荐阅读
- 一些喜欢瞎扯淡的不怎么好用的生成式 AI 工具（没用多少）