重采樣常用於音頻處理。在用麥克風對音頻進行采集的時候,常見的采樣率有8k(電話)、44.1k(CD)、48k(視頻音軌)、96k/192k(Hi-Res),而某些系統會有默認固定的輸出采樣率(如Android的默認輸出采樣率為44.1k),此時就需要對輸入音頻數據進行重采樣。
重采樣的源樣本序列為$x[n]$
$x[n] = x_c(nT)$
重采樣的目標序列為$x'[n]$
$x'[n] = x_c(nT')$
如何通過$x[n]$得到$x'[n]$就是本文的討論內容。
本文假設以采樣周期為$T$對$x_c(t)$進行采樣滿足奈奎斯特采樣定律。
減采樣(downsampling)
減小采樣率的過程被稱為減采樣,這一小節討論的是按整數倍減小采樣率。
按照我們一般的思維來說,按整數倍(倍數為$M$)減少采樣率應該是直接對源樣本序列每隔$M$個樣本提取一個值
$x_d[n] = x[nM] = x_c(nMT)$
這種提取方法被稱為采樣率壓縮器,簡稱壓縮器(compressor)。可以看到所得的新序列是原始連續信號的一部分,並且新序列的采樣周期為$T_d = MT$。對於該新序列,我們可以分為兩種情況進行討論:
- $T_d$符合奈奎斯特采樣定理,即新序列能通過一個低通濾波器還原為原始的連續信號
- $T_d$不符合奈奎斯特采樣定理,即新序列發生混疊,無法還原為原始的連續信號
如下圖假設信號在$M=2$時恰好滿足奈奎斯特采樣定理,那么在$M=3$時則會發生混疊
如果在采用了壓縮率為$M$的壓縮器后,序列仍然符合奈奎斯特采樣定理,我們可以直接進行使用$x_d[n] = x[Mn]$來得到減采樣序列。而發生混疊的情況則稍微復雜一點。觀察混疊的頻譜,可以發現只有低頻部分保持了與原始信號頻譜的一致性,而相當多的高頻由於混疊而失去了原始頻譜。
頻譜丟失得越多說明信號的失真越大,因此為了減少失真,需要盡可能保留更多的原始信號頻譜。我們可以先對元素信號進行低通濾波,然后再對濾波后的信號進行周期為$MT$的采樣即可得到失真更少的序列。按照這種思想,采樣周期固定為$MT$,如果一個被采樣信號的采樣周期為$MT$,那么采樣后不會混疊的條件就是該信號的截至頻率為$\frac{\pi}{MT}$,因此低通濾波的截至頻率為$\frac{\pi}{MT}$。
很明顯,先低通濾波后采樣的這種方法能最大限度地減低頻譜的丟失,從而降低信號失真。
$\tilde{X}_d(j\Omega) = X_c(j\Omega)H_{M}(j\Omega)$
其中$\tilde{X}_d(j\Omega)$就是對原始信號進行低通濾波后的信號$\tilde{x}_d(t)$的傅里葉變換,低通濾波器的傅里葉變換為$H_{M}(j\Omega)$。不過我們手中的並不是原始信號$x_c(t)$,而是序列$x[n]$,為了進入上述流程,我們需要先用$x[n]$重構出$x_c(t)$,才能進行低通濾波以及后續采集
$\begin{align*}
\tilde{x}_d(t)
&= \mathcal{F}^{-1}\tilde{X}_d(j\Omega)\\
&= \mathcal{F}^{-1}\{X_c(j\Omega)H_{M}(j\Omega)\}\\
&= \mathcal{F}^{-1}\{X(e^{j\Omega T})H_r(j\Omega)H_{M}(j\Omega)\}\\
&= \mathcal{F}^{-1}\{X(e^{j\Omega T})TH_{M}(j\Omega)\}
\qquad
\left\{\begin{matrix}
H_r(j\Omega)&=\left\{\begin{matrix}T, & |\Omega|<\pi/T\\ 0, & else\end{matrix}\right.\quad\\
H_M(j\Omega)&=\left\{\begin{matrix}1, &|\Omega|<\pi/MT \\ 0, & else\end{matrix}\right.\\
\end{matrix}\right .\\
&= \mathcal{F}^{-1}\{X_s(j\Omega)TH_M(j\Omega)\}\\
&= x_s(t)*[MT\cdot h_m(t)]/M\qquad fourier\ convolution\ theorem\\
&= \left\{\sum_{n=-\infty}^{\infty}x[n]\delta(t-nT)\right\}*\left\{\frac{sin(\pi t/MT)}{\pi t/MT} \right\}/M\\
&= \sum_{n=-\infty}^{\infty}x[n]\frac{sin[\pi(t-nT)/MT]}{\pi(t-nT)/MT}/M
\end{align*}$
然后從$\tilde{x}_d(t)$中以$MT$為周期采集得到$\tilde{x}_d[n]$
$\color{red}{\begin{align*}
\tilde{x}_d[n] &= \tilde{x}_d(nMT)\\
&= \left .\sum_{k=-\infty}^{\infty}x[k]\frac{sin[\pi(t-kT)/MT]}{\pi(t-kT)/MT}/M\right|_{t=nMT}\\
&= \sum_{k=-\infty}^{\infty}x[k]\frac{sin[\pi(nMT-kT)/MT]}{\pi(nMT-kT)/MT}/M\\
&= \sum_{k=-\infty}^{\infty}x[k]\frac{sin[\pi(nM-k)/M]}{\pi(nM-k)/M}/M
\end{align*}}$
如果要從連續時間系統來理解的話,我們可以發現重建的連續信號$\tilde{x}_d(t)$是由無數個$x[k]$分別對相應位置的截至頻率為$\pi/MT$的sinc函數進行加權后疊加,然后對疊加得到的信號的幅度除以$M$得到的。
如果要從離散時間系統來理解$\tilde{x}_d[n]$的構造公式的話,可以發現上面公式中的sinc函數有如下規律(下二)
$\begin{matrix}
H(e^{j\omega}) =\left\{\begin{matrix}M, &|\omega|<\pi/M \\0, &else \end{matrix}\right .& \Leftrightarrow &
h[n] = \displaystyle{\frac{sin[\pi n/M]}{\pi n/M}}\\
H(e^{j\omega}) =\left\{\begin{matrix}1, &|\omega|<\pi/M \\0, &else \end{matrix}\right .& \Leftrightarrow &
h[n]=\displaystyle{\frac{sin[\pi n/M]}{\pi n/M}/M}
\end{matrix}$
可以發現該sinc函數是一個增益為1、截至頻率為$\pi/M$的低通濾波器,該濾波器后接因子為M的壓縮器。那么,前面的公式可以理解為:用該濾波器對$x[k]$進行濾波后再用因子為M的壓縮器即可得到$\tilde{x}_d[n]$。這一過程又被稱為抽取(decimation)。
增采樣(upsampling)
增加采樣率的過程被稱為增采樣,這一小節討論的是按整數倍增加采樣率。
假設增加采樣率的倍數為$L$,那么增加采樣率后的采樣周期為$T_i = LT$,有
$x_i[n] = x_c(nT_i) = x_c(nLT)$
不過實際上我們只有采樣周期為$T$的序列$x[n]$,因此需要先進行$x_c(t)$的重建,然后對$x_c(t)$進行周期為$T_i$的采樣
$\color{red}{\begin{align*}
x_i[n] &= x_c(nT/L)\\
&= \left .\sum_{k=-\infty}^{\infty}x[k]\frac{sin[\pi(t-kT)/T]}{\pi(t-kT)/T}\right|_{t=nT/L}\\
&= \sum_{k=-\infty}^{\infty}x[k]\frac{sin[\pi(nT/L-kT)/T]}{\pi(nT/L-kT)/T}\\
&= \sum_{k=-\infty}^{\infty}x[k]\frac{sin[\pi(n-kL)/L]}{\pi(n-kL)/L}
\end{align*}}$
下面我們從頻域來展開討論
增采樣即采樣頻率提高了,因此不會出現混疊的情況,即不會改變原始的連續時間信號。現假設用采樣周期$T$對$x_c(t)$進行采樣恰好滿足奈奎斯特采樣定理,即有
增采樣后采樣周期為$T_i = T/L$,那么頻譜將有如下變化
觀察兩張圖右下角的$X(e^{j\omega})$以及$X_i(e^{j\omega})$,它們的時域表示分別就是我們的源序列$x[n]$與目標序列$x_i[n]$。從$X(e^{j\omega})$變為$X_i(e^{j\omega})$只需進行執行兩個步驟:
- 對$X(e^{j\omega})$的變量$\omega$進行倍數為$L$的擴展,得到$X_e(e^{j\omega}) = X(e^{j\omega L})$
- 對所得的新頻譜進行低通濾波,濾波器的增益為$L$,截至頻率為$\frac{\pi}{L}$
對第一步有如下分析:
$\begin{align*}
X(e^{j\omega}) &= \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}\\
\Rightarrow X(e^{j\omega L}) &= \sum_{k=-\infty}^{\infty}x[k]e^{-j\omega Lk}\\
&= \sum_{n/L=-\infty}^{\infty}x[n/L]e^{-j\omega n} \qquad n=Lk\\
\Rightarrow x_e[n] &= \left\{\begin{matrix} x[n/L], &n = 0,\pm L, \pm 2L\cdot\cdot\cdot\\0, &else\end{matrix}\right.
\end{align*}$
如下圖所示,對序列$x[n]$進行步長為$L$的擴展,即可得到$x_e[n]$,這種轉換稱為擴展器(expander)
第二步的低通濾波器的增益為$L$、截至頻率為$\pi/L$,即其脈沖響應為(從前面的增采樣公式同樣也能得出該結論)
$h_i[n] = \frac{sin(\pi n/L)}{\pi n/L}$
因此增采樣系統分解如下圖,這一過程又被稱為內插(interpolation)。
簡單的內插濾波采樣
在對序列進行內插時,某些情況下,我們並不用追求很准確的信號還原,此時用一些簡單的濾波器即可達到不錯的效果,如常見的線性內插、三次樣條內插等。
$h_{lin}[n] =\left\{ \begin{matrix} 1-|n|/L , & |n|\leqslant L\\ 0, &else\end{matrix} \right.
\qquad\qquad h_{cu}[n] = \left\{ \begin{matrix}
(a+2)|n/L|^3 -(a+3)|n/L|^2+1, &0\leqslant |n|\leqslant L\\
a|n/L|^3-5a|n/L|^2+8a|n/L|-4a, &L\leqslant |n|\leqslant 2L\\
0, &else
\end{matrix} \right.$
- 線性內插只用到內插的位置兩旁的兩個樣本,在圖上表示的話,就是把兩個樣本用直線連接后,該直線在對應內插位置上的的值就是所求的內插值。
- 三次樣條內插會用到內插位置兩旁的四個樣本,即左右各兩個,曲線類似於對低通濾波器進行截取后的曲線(下面取$L=5, a=-0.5$)。
觀察上面兩個脈沖響應,我們能發現這些簡單內插器脈沖響應的一些規律:
- 零點的左右兩邊對稱,$\tilde{h}[n] = \tilde{h}[-n]$
- 長度為$2KL-1,K=1,2,\cdot\cdot\cdot$,也就是說左右兩邊距離零點超過KL時的值為0,即$\tilde{h}[n] = 0,|n|\geqslant KL$·
- 零點處的值為1,$h[0] = 1$
- KL處的值為0,$h[KL] = 0$
既然有這些簡單但不夠准確的內插器,那么就應該有相應的方法來判斷這些內插器的內插效果。內插效果可以通過觀察這些內插器的頻譜來進行分析。
對比增采樣所需要的低通濾波器跟簡單的線性內插器頻譜,可以發現在頻譜的$\omega<\frac{\pi}{L}$處本應做大小固定為$L$的增益,但是線性內插在這部分有衰減,另外在$\omega>\frac{\pi}{L}$處本應該是對原始序列頻譜進行截斷,但是由於線性內插器的頻譜在$\omega>\frac{\pi}{L}$仍然有較大的能量(幅度較高),因此可見線性內插可能不會獲得很令人滿意的內插效果。不過如果原本的采樣頻率就遠大於原始信號的截至頻率,則表明采樣所得的序列的頻譜較為集中,此時采用線性內插則會得到較好的效果。
對比線性內插以及三次樣條內插,可以發現在$\omega<\frac{\pi}{L}$處,三次樣條的旁瓣較寬,而在$\omega>\frac{\pi}{L}$處,三次樣條的能量更小(幅度更低),即三次樣條更接近於低通濾波器,因此三次樣條內插會比線性內插的效果更好。
圖中藍色實線為線性內插頻譜,紅色虛線為三次樣條內插頻譜。為了更突出兩者的頻譜差異,右圖對頻譜的幅度進行了對數運算。
非整數因子重采樣
我們前面討論了減采樣以及增采樣,它們的系統分別被稱為抽取器與內插器。對於非整數因子,我們可以通過先接一個內插器,后接一個抽取器這種級聯的方式來得到該重采樣系統。例如,若$L=100,M=101$,則重采樣的周期為$1.01T$
由於中間兩個都是低通濾波器,因此可以合並成一個低通濾波器,取其中的最小值$min(\pi/L, \pi/M)$作為截至頻率,增益為L。