上篇(基於sinc的音頻重采樣(一):原理)講了基於sinc方法的重采樣原理,並給出了數學表達式,如下:
(1)
本文講如何基於這個數學表達式來做軟件實現。軟件實現的細節很多,這里主要講核心部分。函數srcUD()和filterUD()就是實現的主要函數(這兩個函數是在源碼基礎上作了一定的改動,核心思想沒變)。srcUD()是實現一幀中點的重采樣,一個點一個點的做。filterUD()被srcUD()調用。數學表達式就體現在函數filterUD()里。粗看肯定會懵,怎么也跟上面的數學表達式聯系不起來。下面就講講實現細節,讓代碼和數學表達式聯系起來。
先看函數srcUD( ),下圖是其實現,是將一幀中原采樣點的值轉換成新的采樣率下的值,圖中對輸入參數的意思已做了解釋。
主要看3個細節。標1的紅框處是算新的采樣率下兩個點之間的采樣間隔。設定Tx為原采樣率下的采樣間隔,Ty為新的采樣率下的采樣間隔,factor = Tx / Ty,所以 Ty = Tx / factor。軟件實現時對采樣間隔做歸一化處理,設定Tx = 1,即原采樣點x(n)的采樣間隔歸一化到1,則新采樣點y(m)的采樣間隔為Ty = 1/factor。把采樣間隔用Q0.15的定標表示,所以原采樣率下的采樣間隔為32768,新的采樣率下的采樣間隔是dtb = dt *(1 << Np) +0.5 (加上0.5是做四舍五入)。標2的紅框處是算數學表達式1中的D。D = A = min(1, factor) = 代碼中的dh,同dtb,再用Q0.15定標就是dhb。標3的紅框處是具體算生成的新采樣率下的點的值。這一幀開始時間是*Time,原采樣率下采樣點間隔已歸一化為1,定標后變成了32768,幀內有Nx個點,所以這幀的結束時間是endTime = *Time + (1<<Np)*(WORD)Nx。每生成一個新的采樣率下的點后時間就會向后移動一個新的采樣率下的采樣間隔(即dtb)開始算新的采樣點的值,直到時間超過了endTime。算每個新的采樣點時先要找到它對應的數學表達式中的x(Km),即對時間取整,代碼中就是Xp = &X[*Time>>Np],這里Xp就是對應的x(Km)。(*Time&Pmask) 就是算Δm。算重采樣后的一個新產生的點的值時,先算左翼(left wing)的值,再算右翼(right wing)的值,然后兩個值相加就得到了這個新的采樣點的值了。關於左右半翼,下文算sinc函數值時會講。算左右半翼值的函數是一樣的,即filterUD()。先看看怎么把算左右半翼的值用一個函數表示的。回到上面的數學表達式1,假設sinc函數左右都取hend個過零點,上式就變成了
可以看出把計算變成了兩部分相加,前半部分是左翼,后半部分是右翼,即:
所以y(m) = y_left + y_right。着重看右翼。令i=k+1,則
再令i = -i,則
由於sinc函數是偶對稱函數,即sinc(-x) = sinc(x)。所以上式變成
為了公式統一,令 Δ′m= 1 - Δm,i=k,則
再看y_left的表達式:
可以看出計算左右兩翼的表達式統一了,只是傳給函數的參數有差異。
再來看函數filterUD( ),下圖是其實現,是生成新的采樣率下的一個點的半翼的值,圖中對輸入參數的意思已做了解釋。
再看上面的數學表達式1。在式子中,有變量Km和Δm,當新舊采樣率以及m確定時它們就是已知變量,在原理篇中講過它們的求法。y(m)是A、sinc((k+Δm)D)和x(km - k)的乘累加,A和x(km - k)均很容易得到,唯獨sinc((k+Δm)D)需要計算求出。下面看怎么求sinc((k+Δm)D)。
為了減少計算量,求sinc函數值通常采用查表的方法,即先做好表,使用時根據索引查表得到函數值。sinc函數是連續的,計算機處理時要先變成離散的。sinc函數是偶函數,關於Y軸對稱,右半軸(或叫右翼,right wing)的值得到了,左半軸(或叫左翼,left wing)的值也就得到了。這樣也就減小了表的大小(表的大小減為一半)。在右翼中,對每個過零點之間做256次采樣,就得到了sinc函數值的離散表示。以從0到第一個過零點為例,256均分X軸從0到1的函數值,如下圖:
可以根據原理篇中sinc函數的表達式算出每個采用點上的函數值,具體如下:
類似的可以算出右軸的每兩個過零點之間的樣本的函數值。如果截斷后右翼共有6個過零點,則表中共有1536(1536 = 256 * 6)個值,分別為(1,0.99997458,0999989969, ........),把這些值叫做sinc函數濾波器系數。開源實現中分small filter和large filter兩個表,兩個表中過零點之間均是256采樣,但small filter是6個過零點,表中共1536個值,而large filter是32個過零點,表中共8192個值。由於用了更多的sinc函數值和原采用率下的樣本值,large filter的效果更好,同時運算量也大好多。具體使用時需要評估用哪個表。
sinc函數的表做好了,但是sinc((k+Δm)D)中的(k+Δm)D有可能不落在那些離散的樣本上,即(k+Δm)D不等於1/256、2/256等,而是落在兩個樣點之間,例如落在1/256和2/256之間的藍點,如下圖:
這時的sinc函數值怎么求呢?軟件實現中用了線性插值法。線性插值是一種針對一維數據的插值方法,它根據一維數據序列中需要插值的點的左右鄰近兩個數據點來進行數值的估計。它是根據到這兩個點的距離來分配它們的比重的。如下圖:
已知點(x0,y0)、(x1,y1),如果在x處插值,則y的值可用下式求出:
(2)
這樣sinc函數值就求出來了。
算sinc函數值的表是浮點的,為了減少運算時的load,需要將浮點數定點化,即用Q格式表示。用Q0.15給表做定點化,表就變成了(32767,32766,……..),即代碼中的數組Imp(實現中也把兩個相鄰的樣本的Y值的差(即2式中的y1-y0)做成了表,即數組ImpD。這樣處理減少了運算,屬於典型的用空間換時間)。相應的其他浮點數的值也要做定點化,如Δm。sinc的橫軸值x也要做定點化,也是用的Q0.15,這樣第一個過零點的值定點化后就變為了32768。sinc函數兩個過零點之間共256個樣本,所以樣本的間隔是128(32768 / 256 = 128),即2式中的x1-x0的值。x除以128(即右移7位)就可算出落在哪兩個樣本之間。x & 128表示與相鄰左樣本的距離,即2式中的x-x0的值。y0也是定點化后的值,根據2式,落在兩個樣本之間的sinc的定點化后的值也就得到了。
現在看代碼中的實現細節。標1的紅框處是算Δm*D(Ho = (Ph*(UWORD)dhb)>>Np)以及每次加上D后的值(Ho += dhb)。不管是算左翼還是右翼的值,從表達式看出k都是從0開始,所以sinc(k*D + Δm*D)的第一個值是sinc(Δm*D),即k=0時的值Ho。后面經過一個原采樣點,即k加1,Ho就會加上D(即dhb)得到一個新的k*D + Δm*D。標2的紅框處,Ho>>Na(即Ho右移7位)可知落在哪兩個sinc樣本之間,*Hp = Imp[Ho>>Na]表示左邊的樣本的值,即表達式2中的y0的值。同樣*Hdp = ImpD[Ho>>Na]表示左右兩邊樣本的差值,即表達式2中的y1-y0的值。a = Ho & Amask 表示與相鄰左樣本的距離,即x-x0。上面說過定標后兩個樣本之間的間隔是128,即x1-x0的值是128,代碼中表現出來就是右移7位(>>Na)。所以(*Hdp)*a)>>Na就表示(y1-y0)*(x-x0)/(x1-x0)。t = *Hp表示想要求的值左邊的樣本,即y0,因而t += (((WORD)*Hdp)*a)>>Na就是實現了表達式2,sinc函數值就求出來了。再乘以相對應的原采樣點(即t *= *Xp)就得到了sinc(k*D + Δm*D)*x(km -k)的值。算好一個值后Ho會加上D,原采樣點也會左移或右移一個,從而算下一個值,直到sinc樣本的結束處。把這些算好的值加起來就是半翼的值了,再把左右半翼的值加起來就得到了新采樣率下的一個點的值了。
把函數srcUD()和filterUD()搞明白了基於sinc重采樣的軟件實現就好理解了。