我有一个来自麦克风四面体阵列的四通道音频信号。我希望将它从 48 kHz 上采样到 240 kHz。
是否有首选的音频插值方法?对于音频的特定情况,三次插值(或任何其他)是否比线性插值有任何优势?
假设我使用三次插值,我是分别对每个通道进行插值,还是在所有四个通道上使用双三次插值有什么好处?
我有一个来自麦克风四面体阵列的四通道音频信号。我希望将它从 48 kHz 上采样到 240 kHz。
是否有首选的音频插值方法?对于音频的特定情况,三次插值(或任何其他)是否比线性插值有任何优势?
假设我使用三次插值,我是分别对每个通道进行插值,还是在所有四个通道上使用双三次插值有什么好处?
对于音频的特定情况,三次插值(或任何其他)是否比线性插值有任何优势?
您既不会使用音频。原因很简单:您通常为音频信号假设的信号模型非常“傅立叶”,也就是说,它们假设声音由加权谐波振荡组成,并且本质上是带限的。
线性插值和三次插值都不尊重这一点。
相反,您将使用带有抗成像/抗混叠滤波器的重采样器,这是一个很好的低通滤波器。
让我们退后一步:
当我们有一个在时间上离散的信号时,即在一个规则的时间点格上被采样,它的频谱是周期性的——它每(采样频率)重复一次。
现在,当然,我们很少这样看,因为我们知道我们的采样只能代表的带宽,我们通常只绘制从 0 到的频谱,例如:
S(f)
^
|---
| \
| \ ---
| --/ \
| \------\
+----------------------'---> f
0 f_s/2
现在,实际情况是,事实上,我们知道对于实值信号,频谱与对称:
S(f)
^
---|---
/ | \
--- / | \ ---
/ \-- | --/ \
/------/ | \------\
---'----------------------+----------------------'--->
-f_s2/2 0 f_s/2
但是,由于与“采样实例脉冲序列”相乘的某事物的频谱的周期性性质,该事物在两侧无限重复,但我们通常只“看到” 1. 奈奎斯特区(标记为:
)
: S(f) :
: ^ :
: ---|--- : -------
… : / | \ : / \ …
: --- / | \ --- : --- / \ ---
: / \-- | --/ \ : / \-- --/ \
: /------/ | \------\ : /------/ \------\
-------'----------------------+----------------------'---------------------------------------------'-->
-f_s/2 0 f_s/2 f_s
当我们增加采样率时,我们“只是”增加了观察宽度。只是一个随机的例子:
S(f)
^
---|--- :------
… / | \ /: \ …
--- / | \ --- --- / : \ ---
/ \-- | --/ \ / \-- : --/ \
/------/ | \------\ /------/ : \------\
-------'----------------------+----------------------'---------------------------------------------'-->
-f_s/2 0 f_s/2 new f_s/2 f_s
试试看!拿一个音频文件,让你喜欢的工具向你展示它的频谱。然后,只需在每个样本后插入一个,保存为新的音频文件(python 非常适合此类实验),并显示其频谱。您将在左侧看到原始音频(正半部分)频谱,在右侧看到其镜像!
现在,要摆脱这些图像,您只需对原始奈奎斯特带宽进行低通滤波。
这就是重采样器所做的一切:改变采样率,并确保输出信号中不会出现重复和折叠(别名)。
如果您通过整数因子(例如,48 kHz -> 192 kHz)进行上采样,那么您只需在每个输入样本后插入零,然后插入低通滤波器;真的就是这么简单。
在理想情况下,该滤波器将是一个矩形:让原始带宽保持不变,抑制所有不从那里。具有矩形光谱形状的滤波器在时域中具有(无限!)sinc 形状,所以这就是 sinc 插值(以及为什么它非常完美)。
由于该 sinc 是无限长的,而且您的信号不是,嗯,这不是真正可实现的。但是,您可以使用截断的 sinc 插值。
事实上,即使这样也太过分了:无论如何,您的原始音频具有低通特性!(仅仅是因为在对模拟音频源进行采样之前总是需要抗混叠滤波器;更不用说高频是听不见的,无论如何。)
因此,在插入这些零后,您只需使用“足够好”的低通滤波器。这样可以避免计算工作量,甚至可能比 sinc 的截断更好。
现在,如果您的问题绝对不是整数插值怎么办?例如,240000 / 44800 绝对不是整数。那么该怎么办?
在这种相对良性的情况下,我会选择一个合理的重采样器:首先,我们增加一个整数因子,因此得到的采样率是目标采样率的倍数。我们将按照上面的说明进行低通滤波,将结果信号限制在其原始 44.8 kHz/2 带宽,然后通过应用下采样,即抗混叠将其过滤到目标 240 kHz/2 带宽,然后扔掉的个样本。
真的就是这么简单!
其实我们可以进一步简化:由于抗镜像滤波器在22.4kHz处截止,而抗混叠滤波器只在120kHz之后,后者是多余的,可以去掉,这样一个合理的重采样器的整体结构变成:
上采样 -> 核心过滤器 -> 下采样
(事实上,我们甚至可以应用多速率处理并翻转顺序,大大减少了工作量,但这里太过分了。)
那么,你们这里的价格是多少?对于 44800 Hz 输入,240000 Hz 输出,最小公倍数是 3360000 Hz = 3360 kHz,增加 75 倍,低通滤波器,然后减少 14。所以,你需要一个 1/75 频带低通筛选。使用 python 或 octave 很容易设计一个!
音频专用模数转换器 (ADC) 通常具有内部或外部模拟低通滤波器,并以目标采样频率的倍数对模拟滤波信号进行采样。然后,该高速数字信号由数字抽取滤波器进行低通滤波,并抽取到最终的采样频率。如果我们采用DigiKey目前成本最低的24 位、48 kHz 采样频率 (fs) 音频 ADC,即 Asahi Kasei AK5720,其数据表指出:
AK5720 以 64fs 对模拟输入进行采样。数字滤波器抑制高于阻带的噪声,64fs 的倍数除外。AK5720 包含一个抗混叠滤波器(RC 滤波器),用于衰减 64fs 左右的噪声。
当 fs = 48 kHz 时,数字抽取滤波器(图 1)在其从 28.4 kHz 到 3.044 MHz 的第一个阻带中抑制噪声,不希望的过渡和通带以 3.072 MHz 的倍数为中心。在模数转换之前,这些频带中的噪声已被模拟电阻电容 (RC) 滤波器滤除。RC 滤波器滚降为 6 dB/倍频程。第一个阻带的边缘大约有 7 个八度音程。将 RC 滤波器截止频率设置为 40 kHz 将在 3.044 MHz 处产生约 36 dB 的衰减。由于衰减相对较低,因此需要通过系统设计确保在信号进入 ADC 时这些频段内没有太多噪声。声学信号可能会超过 20 kHz 音频频带,因此某些被拒绝的信号可能不是原本的噪声。在抽样中,
图 1. Asahi Kasei AK5720 24 位音频 ADC 的 48 kHz 采样频率的数字抽取滤波器规格概要。该 ADC 可以选择使用对称或短延迟抽取滤波器,牺牲群延迟平坦度以在大部分音频频带上实现低群延迟。这两个滤波器具有几乎相同的幅度频率响应规格。
通常,使用Nyquist-Shannon 采样定理给出的框架处理数字音频它可以通过样本的 sinc 插值来完美地重建连续时间信号。偏离理想框架会产生噪音,应记录在案。这有助于避免过度设计系统的附加部件。例如,如果 ADC 的抽取滤波器的阻带衰减约为 71 dB,则使用具有 120 dB 阻带衰减的插值滤波器对 48 kHz 信号进行上采样是没有意义的,因为较低质量的滤波器将具有较低的复杂性,并且不会显着增加整体噪声水平。如果处理音频信号不是为了收听,而只是为了分析,例如延迟估计,质量要求可能会更宽松。
我对@CedronDawg 的回答投了反对票,因为它歪曲了音频信号的采样方式,就好像连续时间信号没有带限一样。我赞成@MarcusMüller 的回答,因为它正确解释了在奈奎斯特-香农采样中通过整数因子对信号进行上采样。我不同意线性或三次插值不会用于音频的笼统说法。如果质量要求不是太高或信号带宽远小于采样频率的一半,则两者都可以使用,但是我不支持本应用的分段多项式插值。
因为你的上采样率是一个整数,所以“新频率”对可听频带没有混叠,可听频带之外的图像将听不到。然后人们可能会认为,只有与可听频带中平坦频率响应的不必要偏差才重要,并且由于其他原因,您还希望衰减可听频带之外的频谱图像。正如@MarcusMüller 在评论中所指出的,这些原因可能是为了减少狗的痛苦、节省放大器功率、符合某些规范或减少互相关计算中的错误。我不知道您的应用程序是否会从等波纹(图 3)或最小二乘误差滤波器中受益更多。两种类型都可以设计。在您的应用程序中,线性甚至分段三次(图 1)。
您的上采样因子为 240 kHz / (48 kHz) = 5。这是一个固定比率,这意味着分段线性或分段三次插值将等效于通过在之间添加四个新的零值样本来稀释输入信号每个原始连续样本对,将信号乘以等于上采样因子 5 的“上采样增益因子”以补偿由于信号稀释引起的基带衰减,并使用有限脉冲响应 (FIR) 对结果信号进行滤波筛选。这使得分段多项式插值与@MarcusMüller 的答案中描述的上采样框架兼容。
您可以通过使用分段线性或分段三次插值方法对单位脉冲信号进行插值来获得等效的 FIR 滤波器系数,例如通过此 Octave 脚本进行分段三次 Hermite 插值:
pkg load signal
function retval = hermite_upsample(y, R) # Piece-wise cubic Hermite upsample sequence y to R times its sampling frequency, with output endpoints matching the input endpoints. The cubic polynomial tangents at input samples y[k] and y[k+1] are centered differences (y[k+1]-y[k-1])/2 and (y[k+2]-y[k])/2. The input sequence is assumed zero beyond its endpoints.
retval = zeros(1, (length(y) - 1)*R + 1);
n = 1;
for k = 1:length(y)-1
ykm1 = 0;
ykp2 = 0;
if (k - 1 >= 1)
ykm1 = y(k-1);
endif
if (k + 2 <= length(y))
ykp2 = y(k+2);
endif
c0 = y(k);
c1 = 1/2.0*(y(k+1)-ykm1);
c2 = ykm1 - 5/2.0*y(k) + 2*y(k+1) - 1/2.0*ykp2;
c3 = 1/2.0*(ykp2-ykm1) + 3/2.0*(y(k)-y(k+1));
for x = [0:R-(k<length(y)-1)]/R
retval(n) = ((c3*x+c2)*x+c1)*x+c0;
n += 1;
endfor
endfor
endfunction
R = 240000/48000 # Upsampling ratio
b = hermite_upsample([0, 0, 1, 0, 0], R) # impulse response, equal to the equivalent FIR filter coefficients
freqz(b/R) # Plot frequency response excluding upsampling gain factor
plot(b, "x") # Plot impulse response including upsampling gain factor
脉冲响应b
包括上采样增益因子。得到的等效 FIR 滤波器的阶数相对较低,这意味着它在衰减光谱图像方面不是很有效(图 2)。有关光谱图像的解释,请参阅@MarcusMüller 的答案。
图 2. 分段三次 Hermite 插值在上采样到原始采样频率 5 倍时的质量性能。上图:不包括上采样增益因子 5 的 Hermite 插值的幅度频率响应。频率以目标采样频率表示。底部:包含上采样增益因子 5 的 Hermite 插值的脉冲响应。理想的上采样低通滤波器应在频率 π/5 处具有截止频率,并具有拉伸的sinc 函数脉冲响应(包括上采样增益因子)。
还有其他分段三次 Hermite 插值/样条(有时也称为 Catmull-Rom 样条)变体。此处使用的变体根据其相邻样本计算每个样本的正切,根据我的经验,如果我们仅限于分段三次插值方法,该方法在输入采样间隔上形成三次多项式,则它是音频上采样的不错选择,基于四个周围的输入样本。
在您的情况下,双三次插值将等效于三次插值,因为双三次插值通常是可分离的,并且您不会在原始通道“之间”形成新通道。双三次没有什么好处。
您可以通过使用更长的 FIR 滤波器来获得更好的质量性能(图 3),该滤波器可以使用标准的低通 FIR 滤波器设计方法进行设计,例如通过以下 Octave 脚本:
pkg load signal
N = 290; # Filter length - 1
fs_0 = 48000; # Source sampling frequency
fs_1 = 240000; # Target sampling frequency
R = fs_1/fs_0; # Upsampling ratio
f_max = 20000; # Maximum frequency of interest (Eigenmike em32 bandlimit per release notes v17.0)
weight_passband = 1; # Pass band error weight
weight_stopband = 200; # Stop band error weight
b = remez(N, [0, 2*f_max/fs_1, fs_0/fs_1, 1], [R, R, 0, 0], [weight_passband, weight_stopband]) # Stop band starts at fs_0/2 to prevent aliasing which might give artifacts in delay estimation
freqz(b/R) # Plot frequency response excluding upsampling gain factor
plot(b, "x") # Plot impulse response including upsampling gain factor
图 3. 上述 Octave 脚本的 FIR 滤波器在上采样到原始采样频率的 5 倍时的质量性能。顶部:由上述 Octave 脚本生成的 FIR 滤波器的幅度频率响应,不包括上采样增益因子 5。频率以目标采样频率表示。底部:由上述 Octave 脚本生成的 FIR 滤波器的脉冲响应,其中包括上采样增益因子 5。
FIR 滤波器的计算复杂度将通过在实现中考虑到大多数输入样本将是零值来降低。如果您需要标准低通滤波器设计方法无法保证的插值属性,请查看我对 FIR 滤波器设计的回答:Window vs Parks-McClellan and Least-Squares,尽管我不知道如何处理您的具体上采样率为 5。如果有人这样做,他们应该写一个答案:如何使用 Parks-McClellan 算法设计奈奎斯特插值滤波器?插值属性将允许在每 5 个样本处输出输入样本,从而降低计算复杂度。
如果需要低计算复杂度,请注意扩展多项式评估或霍纳评估与分段多项式插值的直接 FIR 滤波器实现相比,分段三次多项式插值的计算复杂度更高。分段多项式插值方法有效地即时计算直接 FIR 滤波器系数,然后通过使用这些系数对稀释的输入进行滤波来生成每个输出样本。这是低效的,因为对于每 5 个输出样本,使用相同的系数,并且这些系数会被重复重新计算。由于这个原因,具有固定系数的直接 FIR 滤波器方法将是首选。与分段多项式插值相比,它还具有更多可以单独优化的系数,因此在相同有效 FIR 滤波器长度的情况下,您应该能够使用直接 FIR 方法获得质量更好的滤波器。
为了进行公平的比较,我们需要承认在固定整数比上采样场景中,可以使用前向差分法进一步优化分段多项式插值。对于相同的有效滤波器长度,我不知道这是否会比直接 FIR 滤波器实现运行得更快或更慢。至少一个直接整数比上采样 FIR 滤波器将易于优化,并且非常适合并行化和单指令多数据 (SIMD) 架构,并且与分段多项式插值不同,可以轻松缩放到更高质量(更长的滤波器) . 出于这个原因,并且您可能需要分段多项式插值无法提供的高质量,我建议使用直接 FIR 滤波器方法。
对于 FIR 滤波器,可以通过采用多速率 FIR 滤波方法获得进一步的加速,例如先上采样 2 倍,然后再上采样 2.5 倍,对后一个滤波器的频率响应要求更宽松。有很多关于用于插值的多级 FIR 滤波的文献。也许你很幸运,有一篇论文以 5 的上采样率为例:Yong Ching Lim 和 Rui Yang,“ On the synthesis of very sharp decimators and interpolators using the frequency-response masking technology ”,在IEEE Transactions关于信号处理,卷。53,没有。4,第 1387-1397 页,2005 年 4 月。doi:10.1109/TSP.2005.843743. 还有无限脉冲响应 (IIR) 滤波器解决方案,特别是两路全通半带滤波器,但具有相位频率响应失真。如果您目前只是进行原型设计,这可能太多了。
大多数人出于某种原因上采样,不清楚你的目标是什么。
由于您提到数据来自阵列,我怀疑您要么将使用额外的时间粒度来为波束成形提供延迟,要么使用额外的样本来简化时间延迟的测量。
我的回答将涵盖波束成形。正如 Marcus Mueler 的回答所建议的那样,与 5 点插值相比,完全向上多速率上采样的延迟要低得多。如果您正在做的所有事情都是操纵光束,那么您实际上不需要(尽管它没有伤害)进行上采样。如果延迟不是问题,我赞同马库斯的回答。
正如 Cedron 的回答所指出的,一组插值滤波器也可能具有较低的复杂性,如果功耗是一个问题,这可能很重要。
就跨通道插值而言,它可能作为运动补偿方案的一部分工作,但 4 个通道并不能为您提供太多插值。
从根本上说,答案将取决于你为什么以及你有什么限制。它也不仅仅是线性和立方。
如果你能得到一份
Nielsen, Richard O. 声纳信号处理。阿泰克大厦公司,1991 年。
对时域波束形成的细节有很好的处理。
抱歉,MM,我同意 Havakok 的观点:实际上,时域插值解决方案应该也可以做到,并且在计算方面要便宜得多。(假设大多数频率内容低于奈奎斯特)。
我会使用三次插值,所以在原始样本点上没有任何“角落”,这当然是更高频率音调的构造(引入)。
通道绝对应该独立插值。
赛德
跟进马库斯:
我认为实际看到一些示例会有所帮助。
1) 线性插值 - 黑线
2) 三次插值 - 红线
3)傅里叶插值 - 绿线
(这不是 sinc 函数的 FIR 实现。相反,我采用 DFT,对其进行零填充,然后采用逆 DFT。)
首先是脉搏。
似乎是 sinc 函数的东西不是。它是狄利克雷核函数,又名 sinc。[请参阅我的博客文章https://www.dsprelated.com/showarticle/1038.php中从 (28) 开始的“随着 N 变大”部分,以了解它们之间的关系。 接下来是一个大正弦波。它们在这里都是很好的近似值。 这是一个相当平稳的信号。端点彼此靠近以使 DFT 公平。 这是一个相当粗略的信号。端点彼此相距很远,以显示 DFT 的环绕弱点。 那么,哪种插值方法实际上更好呢?显然不是线性的。否则,我猜取决于您的标准。
假设我有一段信号是纯抛物线。三次插值将为您提供精确的插值,而 DFT 方法将提供良好的近似值。假设另一部分在 DFT 帧中具有整数个周期的纯音,则反之亦然。
苹果和橙子。
由于四面体麦克风阵列波束成形,我认为 OP 想要上采样以提高延迟估计粒度。查看这些图表,我相信三次插值可以更好地匹配以彼此的分数时间延迟采样的相同信号,所以我坚持我的答案,但这是另一天的测试。
我也坚持使用它会减少计算量,并且正如 SP 指出的那样,延迟更低。
我在 Gambas 中编写了一个程序来生成这些图表。样本值由滚动条控制,因此非常易于使用。我已经在 Gambas 论坛的Interpolation Methods Comparison Project上发布了源代码。
如果没有 Gambas,则需要安装它。最新版本是 13.3.0。存储库参考是 PPA:gambas-team/gambas3
奥利,
是的,我指的是在点附近引入的涟漪,而不是环绕。我不同意你的观点,它们的位置很大程度上取决于粗网格间距,并且会阻碍延迟匹配方法。它们恰好处于粗采样的奈奎斯特频率(每个样本半个周期),因此将在精细插值采样中引入。
您似乎也忽略了我的抛物线形信号部分的反例,并将您的分析集中在正弦音调上。如果我在任何延迟距离处对抛物线进行粗略采样,我将在采样位置获得抛物线上的点。现在,当我进行三次插值时,插值点将与基础信号完全匹配,因此延迟计算也可以是精确的。(我很注重准确性。)
您都缺少的另一点是 sinc 函数与连续情况有关,它只是离散情况下的近似值。
管道,
是的,由于正在解决的问题,我只处理时域评估,“找到延迟”,本质上是一个时域问题。我的观点基于数学经验,在这种情况下尚未得到严格证实。实际上,我喜欢被证明是错误的(尤其是如果我自己做而不是在其中揉鼻子),因为它会导致学习新事物而不是确认我先前存在的偏见。
奥利、马库斯、罗伯特、派普、
足够诡辩的讨论可以在针头上跳舞的天使的数量,让我们拿一根针,一些天使,然后数一数。请提供您推荐的特定算法,包括任何 FIR 滤波器的大小和系数值。它需要在我的 16 点样本集上工作,但我可以根据需要归零。一个快速的代码示例将是理想的。然后我可以进行一些实际的数值测量并为我的“可忽略的谐波”言论辩护。
这是我的三次插值代码:
Paint.MoveTo(myDW, myDH + myBars[0].Value) 对于 n = 1 到 myCount - 3 p0 = myBars[n - 1].Value p1 = myBars[n].Value p2 = myBars[n + 1].Value p3 = myBars[n + 2].Value c1 = p2 - p0 c2 = 2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3 c3 = 3.0 * (p1 - p2) + p3 - p0 对于 m = 1 至 myDW - 1 v = m / myDW f = p1 + 0.5 * v * (c1 + v * (c2 + v * c3)) Paint.LineTo((n + 1 + v) * myDW, myDH + f) 下一个 Paint.LineTo((n + 2) * myDW, myDH + p2) 下一个 Paint.Stroke()
进步:
我没有 Octave(或 MATLAB),我不使用 SciLab,所以我无法对 Olli 的代码做任何事情。但是我看了一下图片,所以这就是我所做的:
'---- 建造一棵奥利冷杉 Dim o As Integer Dim a, f As Float f = Pi(0.2) '2 Pi / 10 我的奥利菲尔[100] = 1.0 对于 o = 1 到 100 a = f * o myOlliFir[100 + o] = Sin(a) / a myOlliFir[100 - o] = myOlliFir[100 + o] 下一个
公平地说,由于端点不为零,我人为地将它们扩展到完整的 FIR 宽度。请注意,我的计算是有效的,因为我实际上不需要将填充的零与 FIR 值相乘并将它们相加。不过,这种方法需要更多的计算才能实现。
'---- 奥利插值 Dim o, t As Integer 对于 o = 0 到 65 v = 0 s = 95 - o 对于 t = s - 5 到 0 步 -5 v += myCoarseSamples[0] * myOlliFir[t] 下一个 对于 c = 0 到 15 v += myCoarseSamples[c] * myOlliFir[s] s += 5 下一个 对于 t = s 至 200 步骤 5 v += myCoarseSamples[15] * myOlliFir[t] 下一个 myOlliValues[o] = v 下一个
我的样本信号是单颗牙齿。黑线代表真正的连续信号。红线是三次插值,绿线是 FIR 插值。采样是完美的,因此样本值是这些点的信号值。两种插值都使用同一组采样值,并且对基础信号视而不见。
那么,额外的计算会导致更好的拟合吗?
接下来是来自两个不同快照的延迟计算。额外的计算是否使这更准确?我对此表示高度怀疑。
我要延迟做延迟处理。我不确定它是否会为讨论增加很多内容,而且我还有其他更紧迫的事情要处理。
我已经在我发布原始代码的同一个论坛线程中发布了生成后一个图表的程序。
https://forum.gambas.one/viewtopic.php?f=4&t=702
它包含除牙齿之外的其他信号。你会很高兴知道 FIR 技术在纯正弦波上优于三次插值,但并不显着。抛物线形状则相反。那里没有惊喜。
在我看来,没有一个案例表明 FIR 技术所需的额外计算需要额外的工作来显着改善结果。还有许多情况(特别是齿和步)三次插值更接近于基础信号。
我非常鼓励 OP 安装 Gambas 并下载这个程序(假设 Linux 可用)。
这是我实现的第一个 sinc 滤波器,它确实有效。它并不总是比三次插值更好,但是当它这样做时并没有明显更好。但是计算成本要高得多。鉴于 Olli 的 290 长度达到 58 个粗点,每个输出点需要 58 次乘法和 58 次加法,而三次方需要 4 次乘法和 3 次加法(如果您包括计算系数而不是使用,则在这种情况下加上 0.8 次乘法和 1 次加法)查找数组)。
做 12 倍以上的工作是否值得?
我不这么认为,但这是OP的选择。我坚持我的开场白:“实际上,时域插值解决方案应该做得很好,并且在计算方面要便宜得多。”,但我学到了一点。