音频采样率转换的研究与代码实现
音频采样率转换
本文原始版本发布于 https://www.52pojie.cn/thread-1959816-1-1.html，此处进行了适当的精简，同时更新了一下代码（最新代码以 GitHub 仓库为准，或者 Gitee 仓库）。如需直接使用，请查看 crates.io 上的页面。
前言
两年前，我研究了 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 _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 spectrum(filename):
fs, data = wavfile.read(filename)
_spectrum(fs, data, filename)
# 显示wav冲激响应频谱
def impulse(filename, passband=False):
fs, data = wavfile.read(filename)
_spectrum(fs, data, filename, impulse='passband' if passband else '')
# 显示原始f64数据冲激响应频谱
def raw_impulse(filename, fs, passband=False):
data = np.fromfile(filename, np.float64)
_spectrum(fs, data, filename, impulse='passband' if passband else '')
# 显示wav文件spectrogram
def 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：
pub trait Convert {
fn next_sample<I>(&mut self, iter: &mut I) -> Option<f64>
where
I: Iterator<Item = f64>,
Self: Sized;
fn process<I>(&mut self, iter: I) -> ConvertIter<'_, I, Self>
where
I: Iterator<Item = f64>,
Self: Sized,
{
ConvertIter::new(iter, self)
}
}
所有的和插值相关的代码都应该实现这个 Trait，这样为后续的包装提供了好处。注意 next_sample 和 process 都是泛型函数，不兼容动态分发，为了让其支持动态分发，我发现约束 Self 为 Sized 就可以了，因为这样使用 Box<dyn Convert> 编译器可以确定要分配内存的大小。
process 包装了 next_sample ，其接受一个迭代器参数，然后返回一个 ConvertIter 结构，该结构实现了一个迭代器，可以获得转换后的输出。在设计上，当输入迭代器出现 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: Convert,
{
type Item = f64;
#[inline]
fn next(&mut self) -> Option<Self::Item> {
self.cvtr.next_sample(&mut self.iter)
}
}
看上去好像还整上生命周期了，好像很麻烦的样子，其实这是 Rust 的迭代器的基本用法。实际上内部就是简单调用了 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<I>(&mut self, iter: &mut I) -> Option<f64>
where
I: Iterator<Item = f64>,
{
loop {
match self.state {
State::First => {
if let Some(s) = iter.next() { // 尝试读取
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) = iter.next() {
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) = iter.next() {
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::Convert;
#[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 Convert for Converter {
#[inline]
fn next_sample<I>(&mut self, iter: &mut I) -> Option<f64>
where
I: Iterator<Item = f64>,
{
loop {
match self.state {
State::Normal => {
while self.pos >= 1.0 {
self.pos -= 1.0;
if let Some(s) = iter.next() {
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) = iter.next() {
self.buf.pop_front();
self.buf.push_back(s);
self.state = State::Normal;
} else {
return None;
}
}
}
}
}
}
当然还有用来包装这个转换器的代码，就不贴出来了。
转换质量分析
这里给出 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 |
上述都是之前的一些比较简单的实验，实际上按照文献给出的论证，衰减（A）和量化（Q）的关系大致是 \(A=12+12log_2Q\)，\(Q=2^{A/12-1}\)，不过实测会存在一些偏差，这并不奇怪。
因此我给出几个不错的参数搭配选择：
- A=96，Q=128：这个参数恰好达到了 16bit 音频的理论 96dB 动态上限（实际不止），且阶数低，存储占用低，适合快速的转换。
- A=108，Q=256 或 A=120，Q=512：这两组基本上对于有质量需求的 16bit 音频来说是完美符合的了，同时兼顾计算量和存储空间。
- A=144，Q=2048 或 A=156、Q=4096：极致的标准，可以满足 24bit 音频的转换，达到了所谓 Hi-Res 的标准，需要极大的计算量和存储空间。
当然略微偏大的 A 也不会造成问题，只会使计算出的阶数稍微高一点点。
性能测试
从前面的表格可以看到计算压力应该不小，于是我们尝试进行一个性能测试。由于 Rust 官方的 benchmark 尚且不够完善且需要 nightly 版本编译器，因此我们选择第三方性能测试框架。一开始我是选择了 Criterion.rs ，后来改用了更轻量级的 divan。
具体测试内容包括线性插值 1 秒数据的耗时，和 Sinc 插值 A=96、A=120 和 A=144 衰减参数下 44.1k，48k 和 96k 互相转换的初始化滤波器系数耗时和处理 10 毫秒数据的耗时。
一个典型的 benchmark 输出结果如下：
Timer precision: 100 ns
bench1 fastest │ slowest │ median │ mean │ samples │ iters
├─ 0. linear 1s │ │ │ │ │
│ ├─ 44k to 48k 61.29 µs │ 365.1 µs │ 61.39 µs │ 65.22 µs │ 1000 │ 1000
│ ├─ 44k to 96k 124.6 µs │ 488 µs │ 124.7 µs │ 131.5 µs │ 1000 │ 1000
│ ├─ 48k to 44k 58.69 µs │ 274.4 µs │ 58.79 µs │ 62.51 µs │ 1000 │ 1000
│ ├─ 48k to 96k 108.1 µs │ 508.1 µs │ 120.2 µs │ 125.3 µs │ 1000 │ 1000
│ ├─ 96k to 44k 74.29 µs │ 237.1 µs │ 95.69 µs │ 100.5 µs │ 1000 │ 1000
│ ╰─ 96k to 48k 76.19 µs │ 286 µs │ 76.29 µs │ 88.84 µs │ 1000 │ 1000
├─ 1. init a96 │ │ │ │ │
│ ├─ 44k to 48k 302.9 µs │ 730.7 µs │ 305.9 µs │ 322.2 µs │ 100 │ 100
│ ├─ 44k to 96k 302.9 µs │ 674 µs │ 306.2 µs │ 313.8 µs │ 100 │ 100
│ ├─ 48k to 44k 330.1 µs │ 746.5 µs │ 333.8 µs │ 347.1 µs │ 100 │ 100
│ ├─ 48k to 96k 169.7 µs │ 657.6 µs │ 181.9 µs │ 199.6 µs │ 100 │ 100
│ ├─ 96k to 44k 658.1 µs │ 1.385 ms │ 663.8 µs │ 698.3 µs │ 100 │ 100
│ ╰─ 96k to 48k 338.3 µs │ 400.4 µs │ 338.5 µs │ 342.4 µs │ 100 │ 100
├─ 1. proc a96 10ms │ │ │ │ │
│ ├─ 44k to 48k 259.5 µs │ 745.1 µs │ 261.7 µs │ 271.7 µs │ 1000 │ 1000
│ ├─ 44k to 96k 519.2 µs │ 1.99 ms │ 524.2 µs │ 571.1 µs │ 1000 │ 1000
│ ├─ 48k to 44k 259.7 µs │ 782.7 µs │ 260.3 µs │ 279.6 µs │ 1000 │ 1000
│ ├─ 48k to 96k 291.5 µs │ 939 µs │ 293 µs │ 313.8 µs │ 1000 │ 1000
│ ├─ 96k to 44k 520.4 µs │ 1.229 ms │ 524.6 µs │ 563.2 µs │ 1000 │ 1000
│ ╰─ 96k to 48k 288.8 µs │ 851.8 µs │ 291.6 µs │ 309.9 µs │ 1000 │ 1000
├─ 2. init a120 │ │ │ │ │
│ ├─ 44k to 48k 1.684 ms │ 2.722 ms │ 1.708 ms │ 1.783 ms │ 100 │ 100
│ ├─ 44k to 96k 1.686 ms │ 3.151 ms │ 1.73 ms │ 1.795 ms │ 100 │ 100
│ ├─ 48k to 44k 1.837 ms │ 2.578 ms │ 1.873 ms │ 1.924 ms │ 100 │ 100
│ ├─ 48k to 96k 942.2 µs │ 1.57 ms │ 1.005 ms │ 1.02 ms │ 100 │ 100
│ ├─ 96k to 44k 3.67 ms │ 5.418 ms │ 3.768 ms │ 3.907 ms │ 100 │ 100
│ ╰─ 96k to 48k 1.883 ms │ 2.953 ms │ 1.935 ms │ 1.998 ms │ 100 │ 100
├─ 2. proc a120 10ms │ │ │ │ │
│ ├─ 44k to 48k 329.4 µs │ 1.078 ms │ 332.4 µs │ 347.4 µs │ 1000 │ 1000
│ ├─ 44k to 96k 660.1 µs │ 2.017 ms │ 670 µs │ 700.5 µs │ 1000 │ 1000
│ ├─ 48k to 44k 334.9 µs │ 896.9 µs │ 340.6 µs │ 354.8 µs │ 1000 │ 1000
│ ├─ 48k to 96k 368.7 µs │ 1.007 ms │ 372.3 µs │ 390.6 µs │ 1000 │ 1000
│ ├─ 96k to 44k 676.4 µs │ 1.802 ms │ 694.2 µs │ 735.4 µs │ 1000 │ 1000
│ ╰─ 96k to 48k 370 µs │ 1.297 ms │ 376.1 µs │ 396 µs │ 1000 │ 1000
├─ 3. init a144 │ │ │ │ │
│ ├─ 44k to 48k 9.065 ms │ 10.6 ms │ 9.176 ms │ 9.332 ms │ 100 │ 100
│ ├─ 44k to 96k 9.055 ms │ 10.82 ms │ 9.188 ms │ 9.323 ms │ 100 │ 100
│ ├─ 48k to 44k 9.845 ms │ 11.7 ms │ 9.986 ms │ 10.09 ms │ 100 │ 100
│ ├─ 48k to 96k 4.978 ms │ 6.392 ms │ 5.049 ms │ 5.171 ms │ 100 │ 100
│ ├─ 96k to 44k 19.71 ms │ 22.43 ms │ 19.96 ms │ 20.14 ms │ 100 │ 100
│ ╰─ 96k to 48k 10.13 ms │ 11.64 ms │ 10.26 ms │ 10.36 ms │ 100 │ 100
╰─ 3. proc a144 10ms │ │ │ │ │
├─ 44k to 48k 542.1 µs │ 2.096 ms │ 554.1 µs │ 587.5 µs │ 1000 │ 1000
├─ 44k to 96k 1.074 ms │ 2.904 ms │ 1.104 ms │ 1.174 ms │ 1000 │ 1000
├─ 48k to 44k 536.1 µs │ 2.087 ms │ 548.9 µs │ 594.7 µs │ 1000 │ 1000
├─ 48k to 96k 578 µs │ 1.641 ms │ 607 µs │ 627.1 µs │ 1000 │ 1000
├─ 96k to 44k 1.069 ms │ 7.259 ms │ 1.222 ms │ 1.639 ms │ 1000 │ 1000
╰─ 96k to 48k 566 µs │ 1.427 ms │ 600.5 µs │ 625.9 µs │ 1000 │ 1000
注意到 divan 提供最好、最坏、中位数、平均值等数据，我们除了关心平均耗时，还需要注意最坏情况。如果处理 10 毫秒数据的耗时超过了 10 毫秒，则会出现爆音的情况，这个通常是线程无法获得足够的 CPU 资源引起的。因此通常需要提高线程优先级保证系统能调度足够的 CPU 资源。
当然一般在实时音频处理的时候不建议采用计算性能要求过高的参数，通常来说 A=120 就可以了，如果还是不行可以进一步降低到 A=108 或者 96，这个时候绝大部分多年前的 CPU 的性能应该都吃得消了。
代码优化
回到之前线性插值的部分代码（Sinc 插值也一样）：
while self.pos >= 1.0 {
self.pos -= 1.0;
self.last_in[0] = self.last_in[1];
if let Some(s) = iter.next() {
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);
这里其实存在一个小小的问题，就是浮点数累计误差问题。因为我们的转换比率是一个浮点数，而实际上浮点数在加加减减的过程中会逐渐产生累计误差，当这个次数不是很多的时候，误差可能在一个极小的量级，而如果次数越多，则误差越大。尽管 f64 的尾数高达 53 位，足以应对长时间的使用，但我们还是希望在软件层面做到零误差，而把误差的情况交给硬件。
因此可以使用定点数来进行优化，这样可以彻底消除长时间使用造成的累计误差。由于绝大多数情况下，我们使用的采样率都是整数，因此使用整数来确定插值的位置是绝对精准的。例如从 44100Hz 到 48000Hz 的转换，其倍率是 147:160，那么我们可以确定的是每输入 147 个样本，必然输出 160 个样本，那么我们完全可以确定插值的位置肯定在 \(\frac{1}{160}\) 和 \(\frac{159}{160}\) 之间，步长则是 \(\frac{147}{160}\)，所以还可以使用一个长度为 160 的表来存储这个比率，避免每次插值重复计算。
为了避免用户输入过于奇葩的比率关系（比如 1234:1233），导致查找表过大，也可以限制这个比率关系，比如不允许约分后出现较大的分子。同时为了限制较大的转化率，也可以限制比值过大或过小。同时考虑到之前留下的函数接口已经采用了浮点数作为比率，因此需要一个将浮点数转换为分数的方案。
综合前面这些要求，搜索后发现 num_rational 库完美满足。
使用如下代码来判断是否是支持的转换比率：
use num_rational::Rational64;
fn supported_ratio(ratio: Rational64) -> bool {
ratio > Rational64::default()
&& *ratio.numer() <= 1024
&& ratio.ceil().to_integer() <= 16
&& ratio.recip().ceil().to_integer() <= 16
}
然后就可以改造之前的一系列代码了。
// 这里传入的是比率的倒数，且已经通过检测并且约分
fn new(step: Rational64) -> Self {
let numer = *step.numer() as usize;
let denom = *step.denom() as usize;
let mut coefs = Vec::with_capacity(denom); // 插值系数表
for i in 0..denom {
coefs.push(i as f64 / denom as f64);
}
Self {
numer, // 分子，每次的增量
denom, // 分母，每次的减量
pos: 0, // 插值的位置
coefs,
last_in: [0.0; 2],
state: State::First,
}
}
这样再改造之前的插值过程：
while self.pos >= self.denom {
self.pos -= self.denom;
self.last_in[0] = self.last_in[1];
if let Some(s) = iter.next() {
self.last_in[1] = s;
} else {
self.state = State::Suspend;
return None;
}
}
let coef = self.coefs[self.pos];
let interp = self.last_in[0] + (self.last_in[1] - self.last_in[0]) * coef;
self.pos += self.numer;
return Some(interp);
本质和原先一样，都是确保插值的位置在 0 和 1 之间，但是采用精确的分数消除累加误差，同时理论上整数的加减相比浮点数应该更快，不过实际上的区别微乎其微。
还有一处比较重要的优化，就是 Sinc 插值时的性能问题。由于 FIR 滤波器的特点，每一个样本都需要进行一次卷积，采样率确定的情况下计算量取决于滤波器阶数，因此提高卷积的速度则是关键。考虑到卷积的特性和卷积核的对称特点，完全可以并行处理，这样可以将循环展开，一次处理多组数据。不过简单展开循环并不一定有效果，因为编译器可能优化期间已经完成了，因此我们需要考虑更多因素。
由于采样率转换使用的是插值滤波器，系数需要根据位置查表，而不是固定的，因此提高查表的效率就是优化的重点了。按照原先的迭代方法，查表的位置是先从大再到小，又变大的过程。此时对 CPU 的缓存不利。如果能在迭代的过程中让查表的位置从小到大，而不再变小，就有利于 CPU 的缓存。因此我们可以按照中心向两侧扩展的顺序迭代，从而提高缓存的命中率，减少因为查表造成的延迟。
而由于奇数阶的滤波器和偶数阶的滤波器存在一定的区别，我们可以单独处理这个问题，即单独计算偶数阶的中心系数。结合之前的分数比率改进，下面是优化后的插值代码：
fn interpolate(&self) -> f64 {
let coef = self.coefs[self.pos];
let mut interp = 0.0;
let pos_max = self.filter.len() - 1;
let taps = self.buf.len();
let iter_count = taps / 2;
let mut left; // 左指针，向左
let mut right; // 右指针，向右
if taps % 2 == 1 { // 偶数阶，奇数长度
let pos = coef * self.quan;
let posu = pos as usize;
let h1 = self.filter[posu];
let h2 = self.filter[posu + 1];
let h = h1 + (h2 - h1) * (pos - posu as f64);
interp += self.buf[iter_count] * h;
left = iter_count - 1;
right = iter_count + 1;
} else {
left = iter_count - 1;
right = iter_count;
}
let coef = coef + self.half_order;
for _ in 0..iter_count {
let pos1 = (coef - left as f64).abs() * self.quan;
let pos2 = (coef - right as f64).abs() * self.quan;
let pos1u = pos1 as usize;
let pos2u = pos2 as usize;
if pos1u < pos_max {
let h1 = self.filter[pos1u];
let h2 = self.filter[pos1u + 1];
let h = h1 + (h2 - h1) * (pos1 - pos1u as f64);
interp += self.buf[left] * h;
}
if pos2u < pos_max {
let h1 = self.filter[pos2u];
let h2 = self.filter[pos2u + 1];
let h = h1 + (h2 - h1) * (pos2 - pos2u as f64);
interp += self.buf[right] * h;
}
left = left.wrapping_sub(1);
right = right.wrapping_add(1);
}
interp
}
这样优化之后，通过 benchmark 测试发现，单次的处理时间大约减少了 10%，说明还是有一定的提升的。当然这里肯定还存在更好的优化方案，需要进一步的研究了。
进一步优化
上述的优化其实存在一些问题，比如有时候可能确实需要支持非分数比例的转换比率，或者对于例如 1234:1235 这样的比例超出了代码中的限制（为了更小的查找表），所以有必要同时保留这几种情况。目前这一部分的功能仅初步实现，高质量且完整的实现需要进一步重构。
实际上，按照文献中给出的方法，想要进一步得到更好的性能提升是很难的了，因为从一个巨大的系数表中去查询参数并非一项快速的操作。而且对于整数倍数的转换，显然绝大多数的系数都是浪费的。比如从 48000Hz 到 96000Hz，插值的位置不是 0 就是 0.5，而系数表却按照量化数量完整生成了，造成了空间的极大浪费，这也导致查表性能下降，造成时间的极大浪费。
即使是 44100Hz 到 48000Hz 这样的转换，实际上也只有 160 种插值位置，生成的系数表的大小仅为 160 倍的阶数，且精度是几乎没有损失的。而按照之前的操作，为了达到足够的精度，量化需要很高，远远高于 160。因此对于较小的分数比例的转换，有必要采用直接查表的方法，为有限个插值位置提供高速且高精度的系数。
这只需要构造一个二维数组即可实现，具体来说，以 1:2 的上采样为例，插值位置仅有 0 和 0.5，那么我们只需要构造一个 2 行 M 列的数组；对于 2:1 的下采样，则只需要 1 行 M 列。此时的空间占用极小且运行速度极快。生成这样的表的方法如下：
fn generate_fast_lut(len: usize, order: u32, beta: f64, cutoff: f64) -> Vec<Vec<f64>> {
let mut lut = Vec::with_capacity(len);
let i0_beta = bessel_i0(beta);
let half_order = order as f64 * 0.5;
let taps = order + 1;
for i in 0..len {
let pos = i as f64 / len as f64; // 每一个插值的位置
let mut coef_pos = Vec::with_capacity(taps as usize);
for j in (0..taps).rev() { // 计算系数
let pos = pos + j as f64 - half_order;
let coef = if (-half_order..=half_order).contains(&pos) {
let i0_1 = bessel_i0(beta * (1.0 - (pos / half_order).powi(2)).sqrt());
sinc_c(pos, cutoff) * (i0_1 / i0_beta)
} else {
0.0
};
coef_pos.push(coef);
}
lut.push(coef_pos);
}
lut
}
使用方法更是简单，因为完全不需要考虑插值的细节，直接对着滤波器表一通卷积就完事了：
fn interpolate(&self) -> f64 {
self.lut[self.pos]
.iter()
.zip(self.buf.iter())
.map(|(h, s)| h * s)
.sum()
}
大道至简，最快的方法永远是最简洁的方法。而且这适用于绝大多数的现实情况，因为目前主流的采样率基本都是 48000 的倍数，唯一的例外是 44100，但是好在二者的比例关系并不夸张，因此这个方法依然可用，只不过系数表的大小是几十倍，但也在可接受的范围内。
至于性能，看一下 benchmark 就知道了，和前面的对比发现有数倍的提升！
Timer precision: 100 ns
bench1 fastest │ slowest │ median │ mean │ samples │ iters
├─ 0. linear 1s │ │ │ │ │
│ ├─ 44k to 48k 65.69 µs │ 265.9 µs │ 65.79 µs │ 73.17 µs │ 1000 │ 1000
│ ├─ 44k to 96k 109 µs │ 427.2 µs │ 110 µs │ 116.1 µs │ 1000 │ 1000
│ ├─ 48k to 44k 64.69 µs │ 280.6 µs │ 64.69 µs │ 68.78 µs │ 1000 │ 1000
│ ├─ 48k to 96k 108.2 µs │ 394.3 µs │ 109.7 µs │ 122.9 µs │ 1000 │ 1000
│ ├─ 96k to 44k 82.79 µs │ 631.7 µs │ 85.29 µs │ 105.7 µs │ 1000 │ 1000
│ ╰─ 96k to 48k 84.19 µs │ 574.3 µs │ 109.9 µs │ 117.5 µs │ 1000 │ 1000
├─ 1. init a96 │ │ │ │ │
│ ├─ 44k to 48k 824.4 µs │ 1.444 ms │ 934.8 µs │ 957.8 µs │ 100 │ 100
│ ├─ 44k to 96k 1.55 ms │ 2.831 ms │ 1.651 ms │ 1.708 ms │ 100 │ 100
│ ├─ 48k to 44k 798.5 µs │ 3.655 ms │ 822.1 µs │ 909.3 µs │ 100 │ 100
│ ├─ 48k to 96k 5.399 µs │ 14.99 µs │ 6.099 µs │ 6.066 µs │ 100 │ 100
│ ├─ 96k to 44k 1.618 ms │ 2.573 ms │ 1.665 ms │ 1.75 ms │ 100 │ 100
│ ╰─ 96k to 48k 5.399 µs │ 12.19 µs │ 5.549 µs │ 5.851 µs │ 100 │ 100
├─ 1. proc a96 10ms │ │ │ │ │
│ ├─ 44k to 48k 51.49 µs │ 302.3 µs │ 52.19 µs │ 56.48 µs │ 1000 │ 1000
│ ├─ 44k to 96k 99.69 µs │ 315.6 µs │ 100.2 µs │ 108 µs │ 1000 │ 1000
│ ├─ 48k to 44k 49.99 µs │ 172.1 µs │ 50.09 µs │ 51.83 µs │ 1000 │ 1000
│ ├─ 48k to 96k 54.49 µs │ 245.6 µs │ 54.79 µs │ 56.79 µs │ 1000 │ 1000
│ ├─ 96k to 44k 100.2 µs │ 532.7 µs │ 100.8 µs │ 105.8 µs │ 1000 │ 1000
│ ╰─ 96k to 48k 56.29 µs │ 334.1 µs │ 56.49 µs │ 60.94 µs │ 1000 │ 1000
├─ 2. init a120 │ │ │ │ │
│ ├─ 44k to 48k 1.082 ms │ 4.76 ms │ 1.167 ms │ 1.284 ms │ 100 │ 100
│ ├─ 44k to 96k 2.171 ms │ 3.429 ms │ 2.264 ms │ 2.354 ms │ 100 │ 100
│ ├─ 48k to 44k 1.07 ms │ 2.084 ms │ 1.143 ms │ 1.19 ms │ 100 │ 100
│ ├─ 48k to 96k 7.499 µs │ 16.89 µs │ 7.499 µs │ 7.708 µs │ 100 │ 100
│ ├─ 96k to 44k 2.252 ms │ 5.177 ms │ 2.319 ms │ 2.455 ms │ 100 │ 100
│ ╰─ 96k to 48k 7.499 µs │ 10.39 µs │ 7.599 µs │ 7.704 µs │ 100 │ 100
├─ 2. proc a120 10ms │ │ │ │ │
│ ├─ 44k to 48k 62.89 µs │ 266.2 µs │ 64.29 µs │ 71.55 µs │ 1000 │ 1000
│ ├─ 44k to 96k 127.6 µs │ 514.6 µs │ 128 µs │ 133.4 µs │ 1000 │ 1000
│ ├─ 48k to 44k 63.09 µs │ 281.1 µs │ 63.29 µs │ 66.06 µs │ 1000 │ 1000
│ ├─ 48k to 96k 68.89 µs │ 128.4 µs │ 69.09 µs │ 70.38 µs │ 1000 │ 1000
│ ├─ 96k to 44k 126.7 µs │ 494.4 µs │ 127 µs │ 134.3 µs │ 1000 │ 1000
│ ╰─ 96k to 48k 71.79 µs │ 244.8 µs │ 71.99 µs │ 75.8 µs │ 1000 │ 1000
├─ 3. init a144 │ │ │ │ │
│ ├─ 44k to 48k 1.418 ms │ 2.255 ms │ 1.49 ms │ 1.536 ms │ 100 │ 100
│ ├─ 44k to 96k 2.804 ms │ 4.229 ms │ 2.879 ms │ 2.943 ms │ 100 │ 100
│ ├─ 48k to 44k 1.388 ms │ 2.727 ms │ 1.47 ms │ 1.495 ms │ 100 │ 100
│ ├─ 48k to 96k 9.699 µs │ 35.69 µs │ 9.799 µs │ 11.98 µs │ 100 │ 100
│ ├─ 96k to 44k 2.824 ms │ 3.494 ms │ 3.013 ms │ 3.032 ms │ 100 │ 100
│ ╰─ 96k to 48k 9.799 µs │ 15.29 µs │ 9.899 µs │ 10 µs │ 100 │ 100
╰─ 3. proc a144 10ms │ │ │ │ │
├─ 44k to 48k 77.59 µs │ 285.6 µs │ 78.09 µs │ 84.63 µs │ 1000 │ 1000
├─ 44k to 96k 159.6 µs │ 853.1 µs │ 162.2 µs │ 169.8 µs │ 1000 │ 1000
├─ 48k to 44k 77.39 µs │ 221.8 µs │ 77.59 µs │ 82.38 µs │ 1000 │ 1000
├─ 48k to 96k 83.49 µs │ 236.9 µs │ 83.59 µs │ 86.56 µs │ 1000 │ 1000
├─ 96k to 44k 153.3 µs │ 759.8 µs │ 154.1 µs │ 160.4 µs │ 1000 │ 1000
╰─ 96k to 48k 86.19 µs │ 430.8 µs │ 86.29 µs │ 90.97 µs │ 1000 │ 1000
接下来的优化目标可以从滤波器的角度减少计算量，例如使用最优化 FIR 滤波器设计方案，或者是最小相位滤波器等。当然肯定也还有其他方向的方法等待更多的发掘。
结语
有关此次采样率转换的研究，我全部写在这里了。从知识的学习，到代码的编写，看了不少书，也重写了几次代码，学会了很多东西，又花了很长时间写完本文，时间跨度从巴黎奥运会开始前到黑神话悟空发售后，历时一月有余，后续也不断更新优化。本文虽然到这结束了，但是知识的学习远没有结束。有关采样率转换以及数字信号处理还有很多我不知道的知识，正等待进一步的探索。同时本文肯定存在不少疏漏之处，还希望高手能够指出。
本文的全部代码都可在开头提高的仓库找到，如需使用可直接在项目使用 cargo add simple_src 或者手动编辑 Cargo.toml 文件。当然，如果需要成熟的采样率转换方案，那么 r8brain、libsamplerate、libsoxr 都是不错的选择。
尽管这是一项业余的研究，但我认为目前的转换质量、转换速度已经达到一个很不错的水平了。
参考
- 奥本海姆《信号与系统》《离散时间信号处理》
- 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 工具（没用多少）
版本历史
- 2024/08/30：完成初版
- 2024/09：进行了多次修改
- 2024/10/06：完成简化后的版本
- 2024/10/15：更新代码，修改部分内容
- 2024/10/19：添加代码优化方案
- 2024/10/28：更新 benchmark，以及一种分数比例的更优化方案
