采样率转换
前言
两年前，我简单研究了 Direct Sound 和 WASAPI 的使用方法，发现了 WASAPI 的 IAudioClient 接口的 GetMixFormat 可以获取系统内部混音器结构。但是在调用 Initialize 时如果设置的 wav 格式不符合系统内部混音器结构会导致调用失败。然而，在其 flags 参数添加 AUDCLNT_STREAMFLAGS_AUTOCONVERTPCM 就可以解决这个问题。后来发现还有一个 AUDCLNT_STREAMFLAGS_SRC_DEFAULT_QUALITY 说是可以获得更好的采样率转换质量。
因此我开始好奇采样率转换是如何进行的，并在当时进行了一些研究，阅读了一些代码，又照猫画虎对着文献写了些代码，能够跑起来但设计上存在问题，且质量非常一般。由于了解到存在知识盲区，因此就暂时搁置，直到最近我又重新研究了这个话题。
简介
采样率转换（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）领域，建议掌握一定的基础，以便于了解本文内容。推荐书籍有奥本海姆的《信号与系统》和《离散时间信号处理》，当然其他的书籍也是可以的。这些书籍的内容不必全都了解（我自己就没看太多），根据需要即可。下面列出重点内容：
- 傅里叶变换、时域和频域的关系
- 采样定理及相关内容
- 滤波器，特别是 FIR 滤波器
其中 1 点应该算是广为人知的了。2 和 3 点则是对于后面部分的重点，也是本文最重要的部分。
当然还需要一定代码能力，本文不涉及 MATLAB 相关内容，而是使用 Rust 语言。当然不会很复杂的，因为我自己也是属于 Rust 初学者，所以有很多地方都不了解，写这个也是练练手。而且总体代码量不大，也可以轻松的用 C/C++ 等其他语言实现。此外，还需要一定的 Python 基础，因为分析时需要使用。
采样率转换的方法
如果你看过上面说的基本书籍，通常会给出采样率转换相关的内容。一般来说，对于采样率从高到低的过程，称为减采样（或者下采样，downsampling）；反之则称为增采样（或者上采样，upsampling）。而具体的实施需要低通滤波器的配合，这样的过程则称为抽取（decimation）和内插（interpolation）。
上述过程都是针对整数倍数转换的，也就是说只能从 44100Hz 到 22050Hz，或者反过来。而要实现非整数倍的转换，则可以级联内插器和抽取器。比如从 44100Hz 到 48000Hz，可以先找出最大公因数 300，从而得出需要内插系数 L=160，抽取系数 M=147。
为了克服巨大的系数，就出现了多采样率信号处理的方法。比如对于前面的 L=160，可以分解成 L1=40 和 L2=4；而 M=147 可以分解为 M1=49 和 M2=3。当然具体的分解方法有详细的论证和严格的数学证明，这里只是举个例子。这样的好处是可以减少计算量，因为对于带限信号，44100Hz 采样率下的最大有效频率是 22050Hz，在 L=160 的情况下，滤波器为了能抑制镜像频率需要较高的阶数，需要很大的计算量。而先以 L1=40 进行内插可以允许一定的镜像频率，因为在下一级 L2=4 时内插时可以使用相对较高阶数的滤波器来去除这个镜像频率。这样下来总体需要的计算量是远远小于不分解的情况。同理在抽取的时候，第一级抽取可以允许一定的混叠，因为第二级可以进行过滤。
当然以上都是书中提到的方法，在实际的使用中并不一定是最合适的。比如按照上面的描述，我们在实现 44100Hz 和 48000Hz 采样率互相转换的时候，需要复杂的多级系统，实现起来也是相当的麻烦。
比如我们有时候只是为了实现简单，快速的转换，完全可以使用基本的线性插值加 IIR 滤波器。而对于较高要求的情况，则可以使用精度更插值算法，如 sinc 函数插值。文献 https://ccrma.stanford.edu/~jos/resample/resample.html 给出了一种基于 sinc 函数的带限插值的算法，网站 https://src.infinitewave.ca/ 提供了各种常见软件的各项技术指标。
正如文献所说，插值算法有拉格朗日插值、三次样条插值、贝塞尔样条插值等，但是这些算法对于图像的优秀效果并不适用于音频。对于音频来说可听性是最重要的，因此文献提出了一种公共领域的算法，用来实现带限信号的插值。该算法主要使用 sinc 函数插值，并针对软件和硬件设备进行了优化，有较高的性能且实现简单，使用灵活。本文也针对这种算法进行了研究，并分析了具体的细节。
如何判断转换的好坏
对于采样率转换来说，通常有很多的指标，通常是检测噪音的引入程度。细分一下包括对无关频率的影响、混叠和镜像频率的产生等。前面提供的网站给出了详细的技术指标，本文主要测试 1kHz 固定音调的正弦音频，以及频率在 $(0, F_S/2)Hz$ 的以二次函数变化的正弦音频（其中 $Fs$ 表示采样率）。前者通常称为 beep，后者则称为 sweep。
为了能够准确地反映转换的质量，需要使用特定的工具来显示音频文件的频谱（Spectrum）以及 Spectrogram（可译为语谱图，也通常叫做时频图或者频谱图）。反正名字挺乱的不同的人可能说的不一样，建议以英文为主。其中 Spectrum 通常是对于某一时间段（比如 1024 个样本）或者全部样本的傅里叶变换的结果，对于 beep 音频来说很合适；而 Spectrogram 则是关于时间、频率和功率的图，功率用颜色表示，通过短时傅里叶变换实现，非常适合 sweep 音频。
除此之外，对于 sinc 插值方法，还需要使用冲激函数测试其滤波器的性能。冲激函数 $\delta(t)$ 在数学上的表述是当 $t=0$ 时，$\delta(t)=\infty$；当 $t \ne 0$ 时，$\delta(t)=0$；且 $\int_{t=-\infty}^{\infty}\delta(t)dt=1$ 。不过实际上，我们只需要让 $x[0]=1$ ，让其他值都为 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：
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;
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;
}
}
}
}
}
说一下设计思路。首先不考虑其他状态，只考虑正常状态，忽略掉相等的情况，转换比率有两种情况——大于 1 和小于 1，但是不论对于哪一种，我们在每一次插值的时候都需要确保插值位置在 0 到 1 之间。对于比率大于 1 的情况，此时步长小于 1，那么同样两个样本可能会插值多次；而当比率小于 1，也就是步长大于 1 的时候，就存在跳过部分样本的情况。因此我们只需要构造这样一个结构，即每次迭代都更新位置，如果位置超出了，就读入一定的样本。
所以这一段代码就很好解释了：
while self.pos >= 1.0 {
self.pos -= 1.0; // 移动1个样本位置
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);
然后考虑第一次输入，需要专门关照一下，因为线性插值无论如何都是需要获取原始输入的第一个样本作为起始的。
if let Some(s) = iter.next() { // 尝试读取
self.last_in[1] = s;
self.pos = 1.0;
self.state = State::Normal; // 切换到正常状态
} else { // 读取失败返回None，此时依然是First状态
return None;
}
这里值得注意的是，在读取到第一个样本的时候，将位置设置成了 1.0，这样可以确保在接下来的代码中可以移进这个样本，而又不会导致计算出现问题。
从挂起状态恢复也是类似的，只有读取到了数据才会切换到正常状态：
if let Some(s) = iter.next() {
self.last_in[1] = s;
self.state = State::Normal;
} else {
return None;
}
不过显而易见的，只根据前后两个样本进行插值计算，显然有较大的误差，尤其是因为声音信号的波形通常都是是由一系列不同频率的周期信号组成，而两个样本不足以表达这些信号的周期。但是线性插值并非一无是处，因为其速度快，而且对于低频信号的影响并没有那么大，所以对于一些性能较差的设备可能还是存在一定的价值。同时，对于线性信号（三角波、锯齿波、方波等）来说，也是有不错的效果。
转换质量分析
这是使用线性插值算法将 44.1k 的 beep 文件转换到 48k 的结果：
对比前面的图，是不是差距还挺大的？其实仔细看可以发现，最大的噪音在 -60dB 以下，剩下的基本都在 -80dB 以下，其余多次混叠的噪音则都在 -100dB 以下，实际上造成的影响并没有那么大。另一张 48k 下采样到 44.1k 的图就不放了，基本上是差不多的。说实话，这样的噪音我反正是听不出来的，想听到得多大音量，怕不是耳朵都要聋了。
而 sweep 就有意思了，因为 sweep 有两个特点，一个是频率不断变化，另一个是频率高点达到奈奎斯特频率，这会导致出现强烈的频率镜像现象和频率混叠现象。
对照原图，可以看到在低频的时候噪音的强度还不突出，而到了高频部分噪音明显变多。看起来是不是很乱？好像每次到了边界就会反弹，而且是不断的反弹，直到信号衰减到检测不到。这其实是由多次的混叠和镜像共同造成的。这回一下就能听出来了，尤其是当声音来到高频的时候，那些低频的混叠就会非常明显。
这就是线性插值的局限性了，因为对于 44100 和 48000，我们之前已经讨论过二者的比例关系是 147:160。从频域的角度来说，线性插值就像是直接先上采样 160 倍，再下采样 147 倍，而上采样的过程就出现了镜像，而且频率高达原先的 160 倍，之后再下采样，又会产生混叠，所以就乱七八糟了。这就是非整数倍转换的一个很大的问题。
为了说明这种情况，这里给出 48000Hz 使用线性插值到 96000Hz 的图像和相反的图像。
总之可以看到，每次到达边界的时候，就会出现一次反弹。具体反弹的次数取决于转换的倍数，倍数越大就会越多。而反弹的位置则取决于是上采样还是下采样。很明显，上采样的时候频率是完全对称的，所以这种现象被叫做镜像。而下采样的时候，高于奈奎斯特频率的信号又会重新出现，这种现象叫做混叠，至于为什么叫混叠，这需要从频域的角度来看了，就不细说了。
但是值得注意的是，对于上采样的镜像来说，通常是听不到的，所以实际上问题不大，而下采样时的混叠会对声音造成明显的可听的影响，这是采样率转换的重点。为了降低干扰，一般情况线性插值还会加入一个 IIR 滤波器，比如经典的巴特沃斯（Butterworth）滤波器。不过 IIR 滤波器会导致非线性相位，在有的时候并不是最好的选择。
由于线性插值并不是本文的核心内容，所以就没有深入分析，只是简单做一个说明，因为其远不如接下来马上要探讨的 sinc 函数插值方法。一般来说，线性插值仅适用于需要快速处理，对质量要求不高的低端嵌入式设备场景。
带限 sinc 函数插值
sinc 函数，定义为 $sinc(x)=\frac{sin(\pi x)}{\pi x}$ ，其中当输入为 0 的时候，结果是 1 。这个函数为什么长这样，这主要还是傅里叶变换的功劳，一般教材都会有推导过程，了解即可。sinc 函数的一大好处就是它是理想的低通滤波器。怎么理解呢，就是在频域上竖着切掉一定频率的信号，只保留低频的信号，对其进行傅里叶逆变换，得到的图像就是 sinc 函数图像。但是为什么说是理想的滤波器呢，因为 sinc 函数只有在正负无穷的时候才会达到 0，显然是过于理想了。
sinc 插值是完美的信号重建方法。因为采样定理告诉我们，对于带限信号，使用低通滤波器就可以重建连续时间的函数（带限信号在频域上是无限的，所以需要低通，这也就是为什么会有混叠和镜像）。
实际上正如之前提到的，整数倍采样率转换实际上只需要低通滤波器就行了。回看上面的两张展示镜像和混叠的图，如果有一个低通滤波器，事情就可以得到解决。而这个低通滤波器，和 sinc 插值就像是提前安排好了一样，为什么不直接使用 sinc 插值的同时滤波呢？
sinc 插值的方法其实很简单，使用这个公式可以实现： $\hat{x}(t)=\sum_{n=-\infty}^{\infty}x(nT_s)h_s(t-nT_s)$ ，其中，$h_s(t)=sinc(F_st)$ 。仔细看可以发现这个公式其实是对整个时域信号做了一个卷积操作，我们知道时域的卷积对应了频域的乘积，所以相当于进行了滤波。但是这是理想的情况，实际上我们不可能有无穷多的点，而且也不可能用那么多的点只是为了计算一个值。
所以我们需要近似，最简单的方法就是把 sinc 函数截断，比如只取前后各 32 个点。又因为截断必然会引起误差（尤其是选择的点数较少的时候），所以还需要加窗。这样下来就是一个 64 阶的 FIR 低通滤波器，只不过这个滤波器的系数表不能提前计算（但是可以近似计算）。
其实和线性插值几乎一样，我们只需要先找到要插值的位置，然后找到这个地方前后一共 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 窗函数都有偶对称的特性，所以只需要保留其一半的系数就可以了，这样又能节省一半的空间。
这里先给出这些数学函数的代码：
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 去哪了？怎么是 sinc_c 了？实际上直接使用 sinc 函数是一个理想的情况，而对于 sinc 低通滤波器来说，其系数计算的公式是这样的：$h[n]=\frac{sin(\omega_c(n-M/2))}{\pi(n-M/2)}\cdot \hat w[n]$ ，其中 $\omega_c\in(0,\pi)$ ，且 $\hat w[n]$ 是窗函数。因此我们使用了上面的 sinc_c ，采用了归一化的截止频率。而零阶修正贝塞尔函数针对迭代情况，展开了公式里的阶乘，从而方便计算。Kaiser 窗也是进行了一定的修改，从而使得函数和 sinc_c 一样是偶对称的。
为了生成一个系数表，我们需要以下几个参数：
-
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 的点。
因为这些参数排列组合可以无限多，所以我们需要一定的理论来预设一些合理的参数，具体方法这里暂且略过，等到后面需要的时候再说。先来说一说如何利用前面这 5 个参数来构造一个插值滤波器系数。
let half = order as f64 * 0.5;
let h_len = (quan as f64 * half).ceil() as usize;
let mut filter = Vec::with_capacity(h_len);
for i in 0..h_len {
let pos = i as f64 / quan as f64;
let coef = sinc_c(pos, cutoff) * kaiser(pos, order, kaiser_beta);
filter.push(coef);
}
是不是很简单，就像之前说的那样。不过这里并没有出现 ratio 参数，实际上滤波器的 cutoff 计算是需要的，后面很快就会说明。下面说说如何利用滤波器系数进行插值。
由于需要保存插值滤波器表和 (M+1) 个点作为缓存，我们使用 VecDeque 作为缓冲区保存之前输入的样本。同时由于不要考虑第一次输入，因此只需要两个状态即可。其他的代码结构基本上是一样的，就不多赘述了。
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,
}
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.fract();
interp += s * h;
}
}
interp
}
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;
}
}
}
}
}```
完整代码会在稍后给出，因为还包含了有关初始化参数计算相关的代码。
### 初始化参数的计算
在教材上有使用窗函数法设计 FIR 滤波器的内容，我们可以按照教材的方法计算我们需要的参数。首先需要确定阻带衰减水平 A，然后据此可以计算 Kaiser 窗的 β 值。
```rust
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
}
}
通常来说我们选择的 A 都是大于 50 的，所以基本上就是按照 $\beta=0.1102\times(A-8.7)$ 计算的。
然后根据公式 $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$，可以看到计算量并没有增加多少，不过滤波器系数的存储空间增加是必然的。
那么量化数量 Q 应该怎么选取呢？前面我们看到当 A=100 时，Q=128 看上去没什么问题，那么将 A 提高到 120，是否会遇到问题呢？按照之前的参数设置，观察 44k 到 48k 情况下的阻带衰减情况。
此时 M=168，看起来依然能达到要求，但是如果将 A 设为 110，却发现效果并不好：
这也是一个很奇怪的问题，实际上 Q=256 才能满足 A=110，而按照这样子应该 Q=512 才能满足 A=120 才对。而实际上按照 Q=128 得出的 Spectrogram 图如下：
仔细对比图上的颜色可以发现有 -110dB 左右的噪音。使用 A=120，Q=512 的结果如下：
这样才符合 120dB 的衰减。怀疑产生这个情况的原因是对于冲激响应测试，A=120 和 Q=128 在数值上可能存在巧合，因此产生了一个极具迷惑性的结果。所以我们需要多种方法共同测试才能得出更加准确的结果。
那么继续提高 A 到 140 呢？为了减少图片数量，就不贴图了，实际上需要 Q=2048 才能满足 A=140 的要求了。
为了给 24bit 音频留出足够的余量，我们通常会使用更高的衰减来满足至少 144dB 的动态要求，因此选择 A=160，看看需要多大的 Q 才能满足。然而发现即使是 Q 来到了 4096 甚至更高，依然无法满足这个要求。
怀疑可能是 wav 文件用的 f32 存在精度瓶颈，因为 f32 实际有效精度为 24 位，因此将 f64 的计算结果转换到 f32 的过程中损失了一定的精度。尽管 scipy 支持 f64 的 wav 文件，但是由于 hound 库不支持，因此尝试直接将原始的 f64 数据写入到文件来分析结果。
由此发现直到 Q=8192 才达到瓶颈，不过此时对存储空间的占用显然是很夸张了，因为滤波器每一阶就需要 8k 个双精度浮点数。此时选择 Q=4096 配合 A=160 也不失为一种选择，只不过实际只能达到 155dB 左右的衰减。如果为了进一步降低 Q，可以选择 Q=2048 和 A=150，此时也能勉强达到要求。
以 A=150，Q=2048，测试一下这么高的标准下的 96k 到 44k 和 48k 的下采样转换。
首先是 96k 到 44k，此时为了 20kHz 的通带，需要高达 464 阶的滤波器才能满足。
当然实际上使用这个标准还是过于极端了，我们用 A=120，Q=512 来试试看 96k 的下采样，此时仍然需要 366 阶滤波器。
看起来其实效果并不差了，而且实际上基本不可能听出来。至于 96k 到 48k 的图，我就不放了，结果是差不多的，因为过渡带选的比较宽（4000Hz），所以就不需要很高的阶数了。
下面的表格总结了一下不同衰减下 96k 下采样需要的滤波器阶数情况，顺便也计算一下 192k 的情况：
|
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 尚且不够完善，我们选择 Criterion.rs 进行性能测试。不过由于网速问题，第一次运行似乎会下载 plotters 这个东西，看起来和卡住了一样，我是等了很久才能正常跑。
由于需要测试的项目繁多，我挑了几个有代表性的，比如 44k 和 48k 的互转，以及 96k 的下采样。内容主要是测试初始化滤波器表所需要的时间，以及生成一秒钟数据需要的时间消耗。为了提高测试压力，我选用了 A=150，Q=2048 下的情况。考虑到误差范围，选择多次测试且不保留任何小数，并且取稳定和较好的结果。使用的 Rust 工具链都通过 Rustup 更新到测试平台默认设置的最新版（目前版本号 1.80.1）。
主要测试代码如下：
use criterion::{black_box, criterion_group, criterion_main, Criterion};
use simple_src::{sinc, Convert};
fn src_bench_init(c: &mut Criterion) {
c.bench_function("src 44k to 48k a150 init", |b| {
b.iter(|| {
let manager = sinc::Manager::new(48000.0 / 44100.0, 150.0, 2048, 2050.0 / 22050.0);
black_box(manager);
})
});
//其他测试
}
fn src_bench_process_1s(c: &mut Criterion) {
let manager = sinc::Manager::new(48000.0 / 44100.0, 150.0, 2048, 2050.0 / 22050.0);
c.bench_function("src 44k to 48k a150 1s", |b| {
b.iter(|| {
let iter = (0..).map(|x| x as f64).into_iter();
for s in manager.converter().process(iter).take(48000) {
black_box(s);
}
})
});
//其他测试
}
criterion_group!(benches, src_bench_init, src_bench_process_1s);
criterion_main!(benches);
由于测试是单声道的情况，如果是立体声则计算量至少翻倍，需要的时间也翻倍。同时为了确保性能，数据都是在插电情况下得出，离电情况下性能通常会打个七八折，这些需要注意。
表格给出了各种测试条件下的滤波器阶数、初始计算量（也反映存储占用）和每秒计算量，方便与测试结果进行比较：
|
44.1->48 |
48->44.1 |
96->44.1 |
96->48 |
| 滤波器阶数 |
213 |
232 |
464 |
238 |
| 初始计算量 |
218112 |
237568 |
475136 |
243712 |
| 每秒计算量 |
10.224M |
10.2312M |
20.4624M |
11.424M |
测试结果如下：
| 平台 |
44.1->48 (init/1s) |
48->44.1 (init/1s) |
96->44.1 (init/1s) |
96->48 (init/1s) |
| 平台 1 |
23ms/40ms |
25ms/40ms |
50ms/88ms |
26ms/35ms |
| 平台 2 |
28ms/56ms |
31ms/56ms |
61ms/137ms |
31ms/53ms |
| 平台 3 |
17ms/138ms |
18ms/139ms |
36ms/367ms |
18ms/150ms |
| 平台 4 |
23ms/86ms |
25ms/87ms |
50ms/280ms |
25ms/80ms |
| 平台 5 |
23ms/79ms |
25ms/81ms |
50ms/240ms |
26ms/76ms |
| 平台 6 |
13ms/105ms |
14ms/105ms |
28ms/215ms |
14ms/114ms |
| 平台 7 |
19ms/77ms |
21ms/77ms |
41ms/152ms |
21ms/79ms |
| 平台 8 |
18ms/150ms |
20ms/150ms |
40ms/298ms |
20ms/172ms |
关于测试平台的详细说明（收集了这么多平台的信息也是不容易啊）：
- 平台 1：高通骁龙 8cx Gen 3 @ 2.69 GHz，Windows on ARM
- 平台 2：高通骁龙 865，红米 K30 Pro，Termux Ubuntu
- 平台 3：AMD Ryzen 5 4600H，Windows
- 平台 4：AMD Ryzen 5 4600H，WSL 2 Ubuntu
- 平台 5：AMD Ryzen 7 4800H，WSL 2 Ubuntu
- 平台 6：Intel Core Ultra 7 155H，Windows
- 平台 7：GitHub Actions，ubuntu-latest
- 平台 8：GitHub Actions，windows-latest
以下是一些补充说明：
- 平台 1 在 Windows 和 WSL 2 Ubuntu 性能几乎一样，故没有分别列出。
- 平台 3 和平台 4 是我自己的电脑，可以发现 Windows 和 Linux 明显差异。
- 平台 5 是找朋友搞的，由于 CPU 差距不大，可以和平台 4 进行对照。
- 平台 6 是另一个朋友搞的，算是较新的 Intel 平台，但是只搞了 Windows 的。
- 平台 7 和平台 8 是 GitHub Actions，可以看出 Windows 和 Linux 的显著差距。
测试发现使用 x86 的平台在计算插值时性能居然不如 arm64 的平台，电脑的性能甚至不如手机。而且 x86 下 Windows 和 Linux 平台也存在不少差异。目前排查发现除去查表造成的缓存问题，主要问题是在插值时使用的 fract 函数。在 arm64 这个函数耗时极短，而在 x86 这个函数就非常浪费，特别是在 Windows 平台，耗时比 Linux 平台多一倍。以后有机会再仔细分析这个函数的性能问题。
补充：最新版本代码性能有一定的提升，且不再使用这个 benchmark 方式了，上述测试仅供参考。
目前使用 divan 进行测试，选用 A=120 和 A=144 这两种情况，以及多种常见采样率之间的互相转换（通带均为 20kHz），一个典型的输出结果如下：
bench1 fastest │ slowest │ median │ mean │ samples │ iters
├─ init_a120 │ │ │ │ │
│ ├─ 44k to 48k 1.684 ms │ 3.035 ms │ 1.698 ms │ 1.76 ms │ 100 │ 100
│ ├─ 44k to 96k 1.684 ms │ 2.613 ms │ 1.704 ms │ 1.745 ms │ 100 │ 100
│ ├─ 48k to 44k 1.834 ms │ 2.157 ms │ 1.849 ms │ 1.868 ms │ 100 │ 100
│ ├─ 48k to 96k 942.7 µs │ 1.355 ms │ 948.7 µs │ 962.3 µs │ 100 │ 100
│ ├─ 96k to 44k 3.67 ms │ 5.186 ms │ 3.705 ms │ 3.758 ms │ 100 │ 100
│ ╰─ 96k to 48k 1.884 ms │ 2.448 ms │ 1.897 ms │ 1.927 ms │ 100 │ 100
├─ init_a144 │ │ │ │ │
│ ├─ 44k to 48k 9.163 ms │ 10.47 ms │ 9.279 ms │ 9.37 ms │ 100 │ 100
│ ├─ 44k to 96k 9.173 ms │ 10.32 ms │ 9.279 ms │ 9.377 ms │ 100 │ 100
│ ├─ 48k to 44k 9.979 ms │ 11.23 ms │ 10.13 ms │ 10.19 ms │ 100 │ 100
│ ├─ 48k to 96k 4.985 ms │ 6.397 ms │ 5.035 ms │ 5.111 ms │ 100 │ 100
│ ├─ 96k to 44k 19.95 ms │ 21.42 ms │ 20.2 ms │ 20.29 ms │ 100 │ 100
│ ╰─ 96k to 48k 10.23 ms │ 11.85 ms │ 10.37 ms │ 10.46 ms │ 100 │ 100
├─ proc_a120_10ms │ │ │ │ │
│ ├─ 44k to 48k 405.5 µs │ 1.054 ms │ 409.3 µs │ 423.4 µs │ 1000 │ 1000
│ ├─ 44k to 96k 815.5 µs │ 1.399 ms │ 827.2 µs │ 854.3 µs │ 1000 │ 1000
│ ├─ 48k to 44k 405.2 µs │ 1.131 ms │ 409.2 µs │ 429.7 µs │ 1000 │ 1000
│ ├─ 48k to 96k 461.1 µs │ 1.076 ms │ 474 µs │ 495.9 µs │ 1000 │ 1000
│ ├─ 96k to 44k 831 µs │ 1.794 ms │ 848.5 µs │ 877.3 µs │ 1000 │ 1000
│ ╰─ 96k to 48k 450.7 µs │ 1.227 ms │ 454.6 µs │ 475 µs │ 1000 │ 1000
╰─ proc_a144_10ms │ │ │ │ │
├─ 44k to 48k 594.4 µs │ 1.599 ms │ 617.9 µs │ 659.8 µs │ 1000 │ 1000
├─ 44k to 96k 1.182 ms │ 5.041 ms │ 1.224 ms │ 1.32 ms │ 1000 │ 1000
├─ 48k to 44k 583.2 µs │ 4.393 ms │ 604.6 µs │ 657.1 µs │ 1000 │ 1000
├─ 48k to 96k 616.2 µs │ 1.249 ms │ 638.1 µs │ 679.2 µs │ 1000 │ 1000
├─ 96k to 44k 1.209 ms │ 5.265 ms │ 1.4 ms │ 1.592 ms │ 1000 │ 1000
╰─ 96k to 48k 592 µs │ 1.472 ms │ 606.1 µs │ 625.8 µs │ 1000 │ 1000
结语
有关采样率转换的研究，我全部写在这里了。从知识的学习，到代码的编写，看了不少书，也重写了几次代码，学会了很多东西，又花了很长时间写完本文，时间跨度从巴黎奥运会开始前到现在黑神话都发售了好几天。这也算是我今年的年度大作了，记得去年也是研究了很久的 IEEE-754 浮点数标准和牛顿迭代法，并完成了软浮点的代码实现。实现这样的作品极其耗费时间和精力，但学习的过程和学到的知识却是非常重要的，这也是我完成这些内容的原因。
本文虽然到这结束了，但是知识的学习远没有结束。有关采样率转换还有很多我不知道的内容，正等待进一步的探索。因此本文肯定存在不少疏漏之处，还希望高手能够指出。
本文全部代码在此 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 工具（没用多少）
致谢
感谢我的朋友不厌其烦配合测试，给我提供宝贵的数据。
[复制链接]
发表于 2024-8-30 17:57
|
发表于 2024-9-2 20:56
