语音信号处理二十八——采样率转换
文章目录
前言一、采样率转换方法1.上采样与下采样2.转换方法1）模拟办法2）数字办法
3.信息丢失
二、整数倍采样率转换1.下采样1）下采样过程2）频谱变化3）下采样后两者DTFT关系4）抗混叠滤波器
2.上采样1）上采样点过程2）频谱变化3）上采样后两者DTFT关系4）滤波去除
三、非整数倍采样率转换1.原理2.转换思路3.滤波器的设计
总结
前言
做过音视频相关工作的同学都知道，采样率转换在不论是音频还是视频编码或者处理中，都经常用到。而采样我们都知道，是单位时间内对于信号等间隔进行量化编码的过程。
我们将单采样率系统称为单速率系统，多采样率系统称为多采样率系统，例如手机、wifi既有语音信号也有图像信号。这也是为什么要进行采样率转换的原因。
本篇文章中是针对于音频采样率转换原理做一个深入介绍。
|版本声明：山河君，未经博主允许，禁止转载
一、采样率转换方法
1.上采样与下采样
当对于一个模拟信号进行采样时：
x
[
n
]
=
x
c
[
n
T
]
=
>
x
d
[
n
]
=
x
c
[
n
T
d
]
x[n]=x_c[nT]=>x_d[n]=x_c[nT_d]
x[n]=xc[nT]=>xd[n]=xc[nTd] 通过不同的采样周期
T
<
T
d
T T f t > f d f_t>f_d ft>fd，此时 将数字信号从 f t f_t ft变为 f d f_d fd的过程，也就是高采样率到低采样率，称为下采样将数字信号从 f d f_d fd变为 f t f_t ft的过程，也就是低采样率到高采样率，称为上采样 2.转换方法 通过对于数字序列进行上采样或者下采样的办法有两种，分别是模拟办法和数字办法，本篇文章中着重介绍的是数字办法。 1）模拟办法 这种转换方法很好理解： 将 x [ n ] x[n] x[n]恢复为模拟信号 x c ( t ) x_c(t) xc(t)对于 x c ( t ) x_c(t) xc(t)进行采样率为 f d f_d fd的采样 2）数字办法 从名字上也很好理解： 上采样：通过插值办法下采样：通过抽样办法 3.信息丢失 在这个过程值得注意的是 操作信息变化备注上采样❌ 不增加信息插值 → 是“估计”，不是“测量”下采样✅ 会丢失信息频率受限，必须滤波，丢掉高频非整数变换❌ 不创造新信息整体来看，本质仍是插值 + 抽取组合 二、整数倍采样率转换 1.下采样 1）下采样过程 对于从高采样率转为低采样率的下采样的方法是抽样，而整数倍采样率转换对应的即是整数倍抽样就是在原序列的基础上进行整数倍的抽取，如： x d [ n ] = x [ D n ] , f d = 1 D f x_d[n]=x[Dn],\quad f_d=\frac{1}{D}f xd[n]=x[Dn],fd=D1f 如上图中给出 D = 2 D=2 D=2，整个时域图像看起来是压缩了，而对于 x d [ 1 ] x_d[1] xd[1]，则实为原序列的 x [ 2 ] x[2] x[2] 2）频谱变化 对于连续时间信号采样后的频谱 X ( e j ω ) X(e^{j\omega}) X(ejω)，是原始模拟信号的频谱 x c ( j Ω ) x_c(j\Omega) xc(jΩ)的无限重复，其幅度是原始模拟信号的频谱的 1 / T s 1/T_s 1/Ts： X ( e j ω ) = 1 T s ∑ k = − ∞ ∞ X c ( j ω T s − j k 2 π T s ) X(e^{j\omega})=\frac{1}{T_s}\sum_{k=-\infty}^{\infty}X_c(j\frac{\omega}{T_s}-jk\frac{2\pi}{T_s}) X(ejω)=Ts1k=−∞∑∞Xc(jTsω−jkTs2π) X ( e j ω ) X(e^{j\omega}) X(ejω)：离散时间傅里叶变换（DTFT），即采样信号的频谱 X c ( j Ω ) X_c(j\Omega) Xc(jΩ)：连续时间傅里叶变换（CTFT），原始模拟信号的频谱 ω \omega ω：角频率，通过对于 Ω \Omega Ω的归一化处理 T s T_s Ts即采样周期 k k k：频谱的周期性重复次数 将抽取的整数倍 D D D带入，上式即变成 X d ( e j ω ) = 1 D T s ∑ r = − ∞ ∞ X c ( j ω D T s − j r 2 π D T s ) , r = D k X_d(e^{j\omega})=\frac{1}{DT_s}\sum_{r=-\infty}^{\infty}X_c(j\frac{\omega}{DT_s}-jr\frac{2\pi}{DT_s}), \quad r=Dk Xd(ejω)=DTs1r=−∞∑∞Xc(jDTsω−jrDTs2π),r=Dk 其中，幅度缩小了 1 / D 1/D 1/D，而采样频率同样缩小了 1 / D 1/D 1/D，所以重复的次数变成了 D k Dk Dk，如同下图 在上文中给出 D = 2 D=2 D=2，此时明显看出在 [ 0 , 2 π / T s ] [0,2\pi/T_s] [0,2π/Ts]区间，下采样后的频谱明显多了一个，也就是 r = D k r=Dk r=Dk。 3）下采样后两者DTFT关系 我们知道时域上关系为 x d [ n ] = x [ n D ] x_d[n]=x[nD] xd[n]=x[nD]，现在使用离散时间傅里叶变化来看一下原序列的频谱和下采样后的频谱之间的关系： X ( e j ω ) = ∑ n = − ∞ ∞ x [ n ] e − j ω n X d ( e j ω ) = ∑ n = − ∞ ∞ x [ n D ] e − j ω n X(e^{j\omega})=\sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n} \\ X_d(e^{j\omega})=\sum_{n=-\infty}^{\infty}x[nD]e^{-j\omega n} X(ejω)=n=−∞∑∞x[n]e−jωnXd(ejω)=n=−∞∑∞x[nD]e−jωn 令 m = n D , n = m / D m=nD,n=m/D m=nD,n=m/D，即 m m m必须是 D D D的整倍数，我们使用 δ D \delta_D δD来对于 m m m进行过滤，即： δ D [ m ] = ∑ k = − ∞ ∞ δ [ m − k D ] \delta_D[m]=\sum_{k=-\infty}^{\infty}\delta[m-kD] δD[m]=k=−∞∑∞δ[m−kD] 因为求和的变化是 i i i，即上式只在 m = 0 , ± D , ± 2 D , . . . , ± k D m=0,\pm D,\pm 2D,...,\pm kD m=0,±D,±2D,...,±kD处为1，使用傅里叶级数展开 ∑ k = − ∞ ∞ δ [ m − k D ] = 1 D ∑ i = 0 D − 1 e j 2 π i m D \sum_{k=-\infty}^{\infty}\delta[m-kD]=\frac{1}{D}\sum_{i=0}^{D-1}e^{j\frac{2\pi im}{D}} k=−∞∑∞δ[m−kD]=D1i=0∑D−1ejD2πim 将整个公式使用 m m m代替 x d ( e j ω ) = ∑ m = − ∞ ∞ x [ m ] e − j ω m D × δ D [ m ] = ∑ m = − ∞ ∞ x [ m ] e − j ω m D 1 D ∑ i = 0 D − 1 e j 2 π i m D x_d(e^{j\omega})=\sum_{m=-\infty}^{\infty}x[m]e^{-j\omega \frac{m}{D}}\times \delta_D[m] \\ = \sum_{m=-\infty}^{\infty}x[m]e^{-j\omega \frac{m}{D}} \frac{1}{D}\sum_{i=0}^{D-1}e^{j\frac{2\pi im}{D}} xd(ejω)=m=−∞∑∞x[m]e−jωDm×δD[m]=m=−∞∑∞x[m]e−jωDmD1i=0∑D−1ejD2πim 整理可得 x d ( e j ω ) = 1 D ∑ i = 0 D − 1 ∑ m = − ∞ ∞ e − j ω m D e j 2 π i m D = 1 D ∑ i = 0 D − 1 ∑ m = − ∞ ∞ x [ m ] e − j m ω − 2 π i D x_d(e^{j\omega})=\frac{1}{D}\sum_{i=0}^{D-1}\sum_{m=-\infty}^{\infty}e^{-j\omega \frac{m}{D}}e^{j\frac{2\pi im}{D}} \\ = \frac{1}{D}\sum_{i=0}^{D-1}\sum_{m=-\infty}^{\infty}x[m]e^{-jm\frac{\omega-2\pi i}{D}} xd(ejω)=D1i=0∑D−1m=−∞∑∞e−jωDmejD2πim=D1i=0∑D−1m=−∞∑∞x[m]e−jmDω−2πi 我们知道对于 ∑ m = − ∞ ∞ x [ m ] e − j m ω − 2 π i D \sum_{m=-\infty}^{\infty}x[m]e^{-jm\frac{\omega-2\pi i}{D}} ∑m=−∞∞x[m]e−jmDω−2πi实际上就是 X ( e j ω − 2 π i D ) X(e^{j\frac{\omega-2\pi i}{D}}) X(ejDω−2πi)的DTFT，所以上式最终为： X d ( e j ω ) = 1 D ∑ i = 0 D − 1 X ( e j ω − 2 π i D ) X_d(e^{j\omega})=\frac{1}{D}\sum_{i=0}^{D-1}X(e^{j\frac{\omega-2\pi i}{D}}) Xd(ejω)=D1i=0∑D−1X(ejDω−2πi) 我们知道对于下采样后的 m m m和原序列的 n n n是一样的，都是 [ − ∞ , ∞ ] [-\infty,\infty] [−∞,∞]，此时我们对于 X d ( e j ω ) X_d(e^{j\omega}) Xd(ejω)进行分析: X ( e j ω − 2 π i D ) X(e^{j\frac{\omega-2\pi i}{D}}) X(ejDω−2πi)：表示每个副本是频谱 X ( e j ω ) X(e^{j\omega}) X(ejω)经过频率缩放和平移得到的 ∑ i = 0 D − 1 \sum_{i=0}^{D-1} ∑i=0D−1：表示 D D D个混叠副本相加 1 / D 1/D 1/D：表示频谱被压缩 D D D倍 最终看起来过程如下图 4）抗混叠滤波器 由上文可知，对于原始序列进行抽样后，原始的分量分量 ω \omega ω变成了 ω / D − 2 π i / D \omega/D-2\pi i/D ω/D−2πi/D ω / D \omega/D ω/D：相当于原始频率在坐标轴上进行了拉伸 − 2 π i / D -2\pi i/D −2πi/D：对于频率分量进行了移位 此时由于原始频率在坐标轴上进行了拉伸，所以可能产生不同频段开始互相重叠，如下图是连续时间信号的CTFT，以及原序列的DTFT和抽样后的DTFT： 对于原始频率在坐标轴上的拉伸，也可以理解为频谱内容压缩 D D D倍，这里可以使用采样定理来解释： 原采样率为： f s = 16 k H z f_s=16kHz fs=16kHz原采样率的Nyquist频率为： f N = 8 k H z f_N=8kHz fN=8kHz下采样因子 D = 2 D=2 D=2下采样后的采样率为 f s ′ = 8 k H z f_s'=8kHz fs′=8kHz下采样率的Nyquist频率为： f N = 4 k H z f_N=4kHz fN=4kHz 此时频率中的高频分量实际上就会叠加到低频分量中，这里因为原始序列中理应只保留 [ 0 k H z , 8 k H z ] [0kHz,8kHz] [0kHz,8kHz]的频率信息，所以我们需要过滤掉 [ 4 k H z , 8 k H z ] [4kHz,8kHz] [4kHz,8kHz]的频率，使用低通滤波器直接过滤掉 4 k H z 4kHz 4kHz以上的频率分量即可，这种滤波器叫做抗混叠滤波器 2.上采样 对于上采样是从低采样率到高采样率，我们使用的方法是插值，插值很好理解，就是在合适的时机插入序列，此时令插值因子是 I I I 1）上采样点过程 对比抽样，插值的过程正好是反过来的，时域序列和采样率直接关系如下： x I [ I n ] = x [ n ] , f I = I f x_I[In]=x[n],\quad f_I=If xI[In]=x[n],fI=If 如上图是和上文中的抽样反过来的，插值因子 I = 2 I=2 I=2，此时对于 x I [ n ] x_I[n] xI[n]和原序列 x [ n ] x[n] x[n]的值为： x I [ 0 ] = x [ 0 ] , x I [ 1 ] = 0 , x I [ 2 ] = x [ 1 ] , . . . , x I [ I n ] = x [ n ] x_I[0]=x[0],x_I[1]=0,x_I[2]=x[1],...,x_I[In]=x[n] xI[0]=x[0],xI[1]=0,xI[2]=x[1],...,xI[In]=x[n] 2）频谱变化 由上文可知，对于连续时间信号的CTFT和采样后离散时间信号的DTFT关系有： X ( e j ω ) = 1 T s ∑ k = − ∞ ∞ X c ( j ω T s − j k 2 π T s ) X(e^{j\omega})=\frac{1}{T_s}\sum_{k=-\infty}^{\infty}X_c(j\frac{\omega}{T_s}-jk\frac{2\pi}{T_s}) X(ejω)=Ts1k=−∞∑∞Xc(jTsω−jkTs2π) 将插值因子带入 X I ( e j ω ) = I T s ∑ k = − ∞ ∞ X c ( j I ω T s − j k 2 π I T s ) X_I(e^{j\omega})=\frac{I}{T_s}\sum_{k=-\infty}^{\infty}X_c(j\frac{I\omega}{T_s}-jk\frac{2\pi I}{T_s}) XI(ejω)=TsIk=−∞∑∞Xc(jTsIω−jkTs2πI) 此时我们对比上采样和下采样 下采样：高采样率到低采样率，过程中幅度缩小 1 / D 1/D 1/D，叠加 D D D次，坐标拉伸 ω / D \omega/D ω/D上采样：低采样率到低采样率，过程中幅度放大 I I I，坐标收缩 I / ω I/\omega I/ω，此时需要去除叠加，实际上就是滤波器的过程 3）上采样后两者DTFT关系 上文中说了对于时域序列 x I [ I n ] = x [ n ] x_I[In]=x[n] xI[In]=x[n]，所以对于原序列和上采样后的序列DTFT分别为： X ( e j ω ) = ∑ n = − ∞ ∞ x [ n ] e − j ω n X I ( e j ω ) = ∑ n = − ∞ ∞ x I [ n ] e − j ω n X(e^{j\omega})=\sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n} \\ X_I(e^{j\omega})=\sum_{n=-\infty}^{\infty}x_I[n]e^{-j\omega n} X(ejω)=n=−∞∑∞x[n]e−jωnXI(ejω)=n=−∞∑∞xI[n]e−jωn 而又由于 x I [ I n ] = x [ n ] x_I[In]=x[n] xI[In]=x[n]存在 对于非 I I I倍的 n n n点， x I [ n ] x_I[n] xI[n]的取值 x I [ n ] = 0 x_I[n]=0 xI[n]=0对于 I I I倍的 n n n点， x I [ n ] x_I[n] xI[n]的取值为 x I [ n ] = x [ n ] x_I[n]=x[n] xI[n]=x[n] 所以上式可以变为 X I ( e j ω ) = ∑ n = − ∞ ∞ x I [ n ] e − j ω n = ∑ n = − ∞ ∞ x I [ I n ] e − j ω I n = ∑ n = − ∞ ∞ x [ n ] e − j ω I n X_I(e^{j\omega})=\sum_{n=-\infty}^{\infty}x_I[n]e^{-j\omega n} \\ = \sum_{n=-\infty}^{\infty}x_I[In]e^{-j\omega In} \\ = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega In} XI(ejω)=n=−∞∑∞xI[n]e−jωn=n=−∞∑∞xI[In]e−jωIn=n=−∞∑∞x[n]e−jωIn 因此上采样后的 X I ( e j ω ) X_I(e^{j\omega}) XI(ejω)中的 x [ n ] e − j ω I n x[n]e^{-j\omega In} x[n]e−jωIn，对比原序列的 X ( e j ω ) X(e^{j\omega}) X(ejω)中的 e − j ω n e^{-j\omega n} e−jωn，实际上就是频率分量压缩了 I I I倍 4）滤波去除 从上文可知，对于插值实际上将频率分量压缩了 I I I倍，也就是对于原始序列 x [ n ] x[n] x[n]的频率分量从 [ − π , π ] [-\pi, \pi] [−π,π]变成了 [ − π / I , π / I ] [-\pi/I,\pi/I] [−π/I,π/I]，而想要去除上文说的叠加，我们知道它是存在 [ − π , − π / I ] [-\pi,-\pi/I] [−π,−π/I]和 [ π / I , π ] [\pi/I,\pi] [π/I,π]之间的，例如当 I = 3 I=3 I=3时： 频率轴（单位：rad/sample）: |-----|-----|-----|-----|-----|-----|-----| -π -2π/3 -π/3 0 π/3 2π/3 π ↑ ↑ ↑ ↑ ↑ 镜像 镜像 主谱 镜像 镜像 而我们实际上需要的只有主谱，这是一个滤波的过程，这个低通滤波器为： H ( e j ω ) = { 1 , ∣ ω ∣ ≤ π I 0 , o t h e r s H(e^{j\omega})=\begin{cases} 1, \quad |\omega|\leq \frac{\pi}{I} \\ 0, \quad others\end{cases} H(ejω)={1,∣ω∣≤Iπ0,others 于是滤波后的频谱为 Y ( e j ω ) = X I ( e j ω ) H ( e j ω ) = X ( e j I ω ) H ( e j ω ) Y(e^{j\omega})=X_I(e^{j\omega})H(e^{j\omega})=X(e^{jI\omega})H(e^{j\omega}) Y(ejω)=XI(ejω)H(ejω)=X(ejIω)H(ejω) 整个过程看起来如下： 三、非整数倍采样率转换 非整数倍采样率转换实际上是转换成整数倍采样率转换来进行的 1.原理 把非整数倍的问题变成“两个整数之间的关系”来处理： 若目标采样率是原始采样率的 I D \frac{I}{D} DI倍（ I , D I,D I,D 为整数），则先插值 I I I倍，再抽取 D D D倍。 也就是说，非整数倍采样率转换 = 插值 + 抽取。而值得注意的是，由于抽取会造成频率混叠，所以需要先进行插值，再进行抽取。 2.转换思路 例如：从44.1kHz 转为 48kHz，整体的采样率转换过程如下： 两个数不是整数倍关系，但我们可以找一个最小整数对： 48 44.1 = 160 147 \frac{48}{44.1}=\frac{160}{147} 44.148=147160插值 I = 160 I=160 I=160，即每个采样点之间插入 159个0低通滤波：使用抗混叠滤波器滤除插值导致的镜像分量抽取 D = 147 D=147 D=147：每间隔147点保留低通滤波：防止混叠 3.滤波器的设计 上文中的转换步骤给出了两次滤波，但实际上只需要进行滤波一次，回顾一下上文对于插值和抽取的频域分析： 插值 I I I倍：在频域复制原始频谱 I I I次（频谱压缩，周期变成原来的 1 / L 1/L 1/L）滤波：设计滤波器的带宽为 ω = min ( π I , π D ) \omega=\min(\frac{\pi}{I},\frac{\pi}{D}) ω=min(Iπ,Dπ)，这个设计既防止了插值产生的镜像，又准备好避免抽取造成混叠抽取 D D D倍：因为滤波器之前已经把带宽限制到 π D \frac{\pi}{D} Dπ，所以不会产生混叠 所以整体的非整数倍采样率转换过程应该如下： x[n] ──↑L──▶ 插值后信号 ──★滤波器★──▶ 抽取后信号 ──↓M──▶ y[n] 总结 本文中并没有给出实际的采样率转换的matlab或者C++的代码，是因为在之前的文章中已经给出了低通滤波器的介绍，并且做音视频开发中采样率转换实际上非常普遍，有兴趣的小伙伴如果设计不好，那么可以一起讨论。 反正收藏也不会看，不如点个赞吧！