能用三角函數表示聲音嗎——正弦模型綜述

摘要：正弦模型(Sinusoidal Modeling)指的是用一系列振幅、頻率和相位不斷變化的正弦波來拟合音頻。其有非常豐富的應用場合。相比于非常成熟的線性預測模型(Linear Prediction)，中文技術社群對于正弦模型的介紹并不足夠。本文參考了十餘篇英文文獻，闡述這一模型的思想和實作思路，解釋其中的技術細節。

文章目錄

• 1.前言
• • 1.1.聲碼器(The Digital Phase Vocoder)
• 2. 分析(Analysis)
• • 2.1. 短時傅裡葉變換(STFT)
  • • 2.1.1.零延拓(Zero-padding)
  • 2.2.峰值檢測與估計(Peak Detection and Estimation)
  • • 2.2.1. 二次曲線/抛物線插值(Quadratic/Parabolic Interpolation)
  • 2.3. 音高檢測(Pitch Detection)
  • 2.4. 峰值連接配接(Peak Continuation/Peak Matching)
  • 2.5. “正弦+噪聲+暫态”模型(Sinusoidal+Noise+Transient)
  • • 2.5.1. “正弦+噪聲”模型(Sinusoidal+Noise)
    • 2.5.2. 暫态模型(Transient Model)
• 3. 合成(Synthesis)
• • 3.1. “僅振幅合成法”(Magnitude-Only Synthesis)
  • 3.2. 相位内插合成法(Phase Interpolation Synthesis)
  • • 3.2.1. 三次多項式内插法(Cubic Polynomial Interpolation)[^McAulay]
    • 3.2.2 相位的意義
  • 3.3. 幀重疊相加法(Overlap-Add/OLA Synthesis)
• 4. 結語
• 注腳

本文由@EthanLifeGreat/@EthanUnbeaten原創發表于CSDN

1.前言

模型的目的在于用相對熟悉的屬性來描繪未知的事物。而将音頻看作是不同的正弦信号的疊加是一個看起來比較自然的想法（其中存在的不合理之處将會在後面提及）。

如果能夠實作，那麼模組化者對于音頻的整體屬性就能形成直覺的了解，同時有利于對音頻進行想要的調整。

在介紹模型之前，請允許筆者先對相關概念進行簡單介紹。

1.1.聲碼器(The Digital Phase Vocoder)

用正弦波來拟合聲波的辦法最早在70年代由Moor提出1。後來經過完善，發展出了聲碼器2。

聲碼器将聲音信号與多個濾波器相卷積（其本質是離散傅裡葉變換）得到多個頻率的子帶(subband / DFT bin)，以及子帶上的振幅和相位。

後來發明的MPEG音頻編碼方式也是基于這個思想——對不同的子帶配置設定不同比特（資訊計量機關）。清晰記錄音頻中人耳敏感的部分；模糊記錄不敏感的部分。

盡管聲碼器的實踐和正弦模型的實踐方式接近，其本質的思維方式是不同的。DPV注重于固定頻率通道(channel)的内容，沒法溝通相鄰的頻率通道，故對刻畫跨通道變化的頻率并不完美。

對此，Smith等人3和McAulay等人4獨立地提出峰值跟蹤(Peak Tracking)的解決方案，正式地将正弦波的概念引進模型。這将是本文讨論的重點。

2. 分析(Analysis)

2.1. 短時傅裡葉變換(STFT)

如前面所言，需要用正弦波來表示音頻信号，必須找到正弦波們的頻率、對應的振幅和對應的相位。對此，我們需要的頻率分析工具是傅裡葉變換(Fourier Transform, FT)，具體而言，是離散傅裡葉變換(DFT)。後者是前者在計算機上的版本。

然而僅僅進行一次DFT是不足夠的，我們首先需要對音頻在時間上進行分幀(frame)，對每一幀獨立地進行DFT，這被稱為短時傅裡葉變換(STFT)/離散短時傅裡葉變換(DSTFT). 這麼做的原因在于：我們希望在每一個幀裡，具有平穩性——例如（在語音分析裡）最多存在一個音素(phone). 誇張地講，如果含有兩個字的音頻片段被同時送入DFT進行分析，我們便難以分析字的内容，也難以辨識兩個字的先後順序。

由于穩定性的影響，我們需要讓每一幀盡可能短。然而對于DFT而言，過短的信号意味着更少的資訊，将會使頻譜的分辨率下降，使得我們難以區分頻率相近的成分。其中，FT的時間分辨率和頻率分辨率滿足反比例關系：分析幀時間越短(時間分辨率高)，則頻率分辨度越低；反之亦然。這被稱為Gabor極限(Gabor–Heisenberg Limit)5.

是以，對幀率（每秒分割的幀數）的選擇也是STFT的一門藝術。新竹清華大學的劉奕汶教授6對此有一句總結：“（窗長）盡量長，不能太長。”

其中的窗長等于幀長，由于截取分析幀後通常要加分析窗，而分析窗的長度與幀長相等。對于分析窗的介紹，可以參考Harris的文章7. 對此我之前的 文章 也有簡要論述。

更詳盡地了解STFT，或可參見Matlab官方封裝的函數 stft.m . 受篇幅限制，本文不再贅述。

2.1.1.零延拓(Zero-padding)

進行離散傅裡葉變換前，為了提高變換後的頻率分辨率——即更接近的DTFT的連續的結果——我們需要對加窗後的信号進行零延拓，即在信号後面補零。

由于FFT的運算性質，一般将補零後的長度為2的某次方（視情況而定，如1024，4096等）。

另外，為了使得相位一緻，應使加窗信号居于DFT輸入的正中間（x取0處），其餘部分為補足的零。如下圖（圖源8）：

從上至下以次為：(a)原信号 (b)加窗信号 (c)延拓信号 (d)信号頻譜。

2.2.峰值檢測與估計(Peak Detection and Estimation)

對幀内容進行N點DFT/FFT後，可以得到一個自變量為頻率的離散的複值函數。函數值的絕對值為頻率處對應振幅，虛部與實部之比的反正切值為相角。

峰值指的是對于振幅譜上的極大值點。

一般認為，隻有較大的峰值才是頻譜中有價值的内容。而那些較最大峰值低80dB以上的峰值，将難以被聽到且其分辨率很低9。對存在這些峰值的音頻分析結果進行調整（如變調）可能帶來音頻品質下降。

同時，某些峰值可能并非真正存在于原本的音頻内容中，而是源于離散傅裡葉變換過程。例如一個被正弦窗截斷的正弦波，經過DFT後會出現許多的峰值：

是以，峰值檢測的流程并不僅僅包括找到極大值點，還在于對找到的極大值點進行篩選，即：

• 删掉不重要的極大值點
• 删掉本身不存在于音頻中的峰值等

一個峰值檢測的結果例子如下圖，叉為峰值（圖源9）：

然而如第一段所言，盡管此函數在圖上看起來連續，但其仍是離散函數，要高效地找到隐藏的極大值點（即峰值估計），可以使用插值法(interpolation).

2.2.1. 二次曲線/抛物線插值(Quadratic/Parabolic Interpolation)

在對數譜上取振幅最值點附近三點，進行二次曲線/抛物線插值(Quadratic/Parabolic Interpolation)，得到估計出的最值點10，再據此對相位進行線性拟合。二次插值示意圖如下（圖源新竹清華大學課件6）：

注意插值進行的縱坐标是分貝dB（對數譜），經驗證明，在對數譜上進行二次插值的精确度比線上性譜上的精确率高一倍8。

另外，除了使用DFT分析出頻率成分，還有學者使用最優化的思想，疊代計算出音頻中主要成分的振幅、頻率和相位，這被稱為“使用合成來分析”11(Analysis-by-synthesis, ABS)12。

2.3. 音高檢測(Pitch Detection)

找到了峰值之後，我們可以選擇進行音高檢測。這可以對之後的分析帶來某些友善。

在此之前我們需要先介紹一些音樂方面的概念。

音高(Pitch)的是音樂領域裡比較模糊的一個概念，筆者比較認同的一個解釋是——聽感上最相近的純音的頻率（a percept that can be compared against that of a pure tone）6。

音高通常被直覺地定義為基頻(Fundamental Frequency)，基頻可以被定義為泛音/諧波(partial/harmonic)序列裡的公因數。例如，鋼琴鍵C4的基頻是261.6Hz，但其頻譜成分卻包含許多（近似于13）261.6的倍數的成分，如下（圖源14）：

基于以上知識，我們不難想象檢測音高（基頻）的算法設計過程。但是，音高檢測又有什麼好處呢？

• 判斷噪音。如果某個峰值離基頻的任何倍數都較遠，那麼這個峰值很可能不屬于這個音。
• 改善分析窗的長度。知道了某個幀中的音的音高，便可以設定合适的分析窗長度，以獲得更好的時間-頻率效益取舍9。這種分析法也被稱為“音高-同步分析”11(Pitch-synchronous Analysis).
• 便于調整音高。如果對音頻音高進行調整時直接将所有的頻率成分乘以某一倍數，則會放大諧波與基頻倍數之間的差距13，使得改變後音頻變得更加不和諧12。

2.4. 峰值連接配接(Peak Continuation/Peak Matching)

峰值連接配接的本質是把相鄰幀内的峰值相連，更确切地說是将相鄰幀内對應頻率的峰值相連，是以也被稱為峰值比對(matching). 而對應頻率一般意味着彼此在各自幀内的頻率最接近。

這麼做的理論基礎在于，我們确信一個穩定的聲音由多個頻率近似穩定，而振幅不斷變化的成分構成。那麼，為了得到對聲音的更連續的描述，我們可以基于這個假設，對相鄰幀内的峰值内插，也相當于去掉了“幀”這個離散的存在。以下是連接配接的示意圖（圖源15）：

其操作大緻流程是，為幀内的每一個峰值在下一個幀内找一個最相近的峰值。在滿足相連接配接的峰值頻率差小于一個給定值的前提下，對産生的沖突（多連一、一連多）按照一定規則進行解決。如果一個峰值在下一個幀找不到對應的連接配接，則被視為是一個軌道(track/trajectory)的死亡(death)；而如果一個峰值沒有被前一個幀中的峰值對應，則被視為是一個軌道的誕生(birth)。

關于峰值連接配接中類似的概念表述有許多，其最終功能近似，不一一列舉。

2.5. “正弦+噪聲+暫态”模型(Sinusoidal+Noise+Transient)

本節簡要解決一下音頻中正弦波不能很好拟合的部分。

前面提到，正弦波隻适合于拟合一個穩定的聲音（如樂音）。而對于噪聲（如雨聲）或者暫态聲音（如打擊音/Attack）都沒辦法很好地拟合。

例如，以下是一段人聲的峰值連接配接結果（圖源15）：

其中的清音部分(Unvoiced Segment)峰值多、軌道短，這樣模型拟合效率低、效果差。而相比下濁音部分(Voiced Segment)顯得合适很多。

2.5.1. “正弦+噪聲”模型(Sinusoidal+Noise)

了解過語音識别的同學可能知道，清音其實就近似于噪音。于是，我們有了“正弦+噪聲”模型，或叫“決定+随機”模型(Deterministic+Stochastic)16。

它首先拟合出正弦波成分，又叫确定性(Determinisitic)成分，剩下的部分稱為剩餘(Residual)。

對于正弦波拟合不好的剩餘，該模型視其為在不同頻率部分能量不同的噪聲17。我們隻記錄其大緻位置的振幅，也就是描繪出它頻譜的包絡，如下圖b子圖中的折線（圖源9）。

這種做法會模糊頻率的精度，同時丢失相位資訊，但好處是可以用較少的資料量表示正弦波無法拟合的噪聲。

2.5.2. 暫态模型(Transient Model)

暫态主要指敲擊瞬間發出的聲音，即ADSR包絡中的（時間很短的）Attack。

這一部分無法很好地用正弦波或噪聲的辦法拟合18。有一種解決辦法是先把暫态時間段剝離出來，不對之進行分析，等合成的時候直接給拼回來19。

Verma等人18則認為這種做法不符合正弦模型的“變通精神”(the flexible spirit)，同時這種在時域上剝離暫态也會帶走某些噪聲。是以，他們提出了對暫态模組化的辦法。

Verma等人的做法是，對暫态部分先進行離散餘弦變換(Discrete Cosine Transform, DCT)20，定義如下：

C [ k ] = β [ k ] ⋅ x [ n ] ⋅ c o s [ ( 2 n + 1 ) k π 2 N ] C[k]=\beta[k] \cdot x[n]\cdot cos[\frac{(2n+1)k\pi}{2N}] C[k]=β[k]⋅x[n]⋅cos[2N(2n+1)kπ​]

其中 k = 1 k=1 k=1 時 β [ k ] = 1 N \beta[k]=\sqrt{\frac{1}{N}} β[k]=N1​

​，否則 β [ k ] = 2 N \beta[k]=\sqrt{\frac{2}{N}} β[k]=N2​

​.

簡單而言，DCT可以将沖激信号轉換為餘弦函數。進而便于進行正弦拟合。如，一個指數速度衰減的信号和其經DCT變換後形式如：

暫态模型分析/合成流程如下：

與“正弦+噪聲”模型相結合時，可以按暫态、正弦、噪聲的順序，也可以按正弦、暫态、噪聲的順序進行分析。

其實寫完這一節才發現有許多的暫态模型，可以參考斯坦福大學CCRMA21的網頁了解更多

3. 合成(Synthesis)

合成的意思是用分析得到的參數化表示，合成出近似于原來的 或 任何想要的音頻。

為了印證“聲音能用三角函數來表示”，我們當然需要用正餘弦波來合成信号。所使用的正餘弦波參數（頻率、振幅和相位）就是我們在分析中得到的資料。

在此，我們主要探讨正弦部分的合成（即不讨論2.5.中的“噪聲”和“暫态”部分）。介紹三種合成辦法，它們分别是：“僅振幅合成法”11(Magnitude-Only Synthesis)、相位内插合成法(Phase Interpolation Synthesis)和幀重疊相加法(Overlap-Add Synthesis)。

3.1. “僅振幅合成法”(Magnitude-Only Synthesis)

顧名思義，這個辦法僅使用分析得到的頻率和對應的振幅（不用相位）進行合成。

這個辦法（相較于相位内插）可以大大簡化合成的流程，了解起來也比較直覺——

比如做實驗的時候，我們想生成一個頻率為f的正弦波，我們一般會忽略相位：

t = linspace(0,1,44100);		% 時間線(s)
A = 1;							% 振幅
f = 1000;						% 頻率(Hz)
signal = A*sin(2*pi*f.*t);		% 頻率1000Hz，相位為0的正弦波
           

但其實相位(0)隐藏在了裡面。

在許多場合下，人耳對于相位是不敏感的，例如你基本沒法用耳朵區分這兩個信号：

signal1 = A*sin(2*pi*f.*t)
signal2 = A*cos(2*pi*f.*t)
           

這就賦予了隻用振幅合成法的現實意義。

“僅振幅合成”的具體的操作辦法是：對于2.4.中連接配接好的每一條軌迹，設其初始相位為0（或者其它某值，效果相同），逐采樣點疊代相位并内插振幅後，取各點正弦值，最後按時間順序将所有軌迹的正弦值對應相加。 其表達式為：

∑ 所 有 軌 迹 j A j [ n ] ⋅ c o s ( ϕ j , n ) \sum_{所有軌迹j}^{}A_j[n]\cdot cos(\phi_{j,n}) 所有軌迹j∑​Aj​[n]⋅cos(ϕj,n​)

其中，（省略 j j j）

ϕ n = ϕ n − 1 + 2 π f s ⋅ f [ n ] \phi_{n}=\phi_{n-1}+\frac{2\pi}{f_s}\cdot f[n] ϕn​=ϕn−1​+fs​2π​⋅f[n]

其中，n表示時域順序，即n=0表示該軌迹的起點， ϕ 0 \phi_0 ϕ0​=0， f s f_s fs​表示采樣頻率（下同）.

A[n]與f[n]是由分析結果（A[m]與f[m]）線性内插而來的振幅和頻率序列，滿足：

A [ n ] = A [ m ] + A [ m + 1 ] − A [ m ] Δ t ⋅ f s ⋅ n A[n]=A[m]+\frac{A[m+1]-A[m]}{\Delta t\cdot f_s}\cdot n A[n]=A[m]+Δt⋅fs​A[m+1]−A[m]​⋅n

f [ n ] = f [ m ] + f [ m + 1 ] − f [ m ] Δ t ⋅ f s ⋅ n f[n]=f[m]+\frac{f[m+1]-f[m]}{\Delta t\cdot f_s}\cdot n f[n]=f[m]+Δt⋅fs​f[m+1]−f[m]​⋅n

這兩條等式非常直覺，就是直線的兩點式，其中 Δ t \Delta t Δt為相鄰分析幀的中心位置時間差。

關于不考慮相位的弊端和适用場合的讨論将放在下一節進行。

3.2. 相位内插合成法(Phase Interpolation Synthesis)

這節我們讨論考慮相位的合成辦法。

考慮相位的合成法需要将相位也進行内插，即确定幀與幀間各點的相位資訊。

3.2.1. 三次多項式内插法(Cubic Polynomial Interpolation)4

不同于振幅内插和3.1.中的頻率内插的線性辦法，相位内插所使用的内插辦法是三次方(cubic)的。

為什麼是三次？首先我們講講為什麼不能是線性（一次）或二次的。這是因為：

• 頻率是相位的導數，要保證相位曲線在兩端的導數值是對應的頻率。
• 分析時得到的相位隻是對 2 π 2\pi 2π 取模的結果，在 [ − π , π ] [-\pi,\pi] [−π,π] 間。

這時我們需要對相位進行加 M ⋅ 2 π M\cdot 2\pi M⋅2π 展開(unwrap)，M 為某待定整數。再用其它先驗知識來求出最合适的M。

例如4文中的圖6，展示了五種可能 M 帶來的内插結果：

如何從中選擇最合适的M呢？作者McAulay等人提出的辦法是，找一個“最平滑的”(maximally smooth)曲線，也就是曲線的二次導數的平方的線積分最小。即求整數M，使得

f ( M ) = ∫ 0 T [ θ ′ ′ ( t ; M ) ] 2 d t f(M)=\int_0^T{[\theta''(t;M)]^2dt} f(M)=∫0T​[θ′′(t;M)]2dt

最小。

得到相位後取餘弦，點乘線性内插後的振幅即可，再把所有軌迹相加即可：4

3.2.2 相位的意義

盡管我們聽不出正弦波和餘弦波的差别，在許多的情況下，人對于相位是敏感的，這些情況是：

• 對象分析過程不當9
• 非常低沉的樂器聲音(低于30Hz9，低于100Hz4)
• 某些人聲(含有較高的泛音)9
• 含有噪音的語音片段(會導緻合成的噪音部分有一種不真實和令人厭煩的音效)4

一個大提琴波形、帶相位和不帶相位合成的結果對比圖如下：16

可以看出帶相位的合成法基本保持了原來的波形（而且幾乎完美複刻），反觀不帶相位的就不知道是哪跟哪了（雖然聽起來可能差不多）。

3.3. 幀重疊相加法(Overlap-Add/OLA Synthesis)

這種辦法事實上并沒有進行峰值連接配接，它隻是單純地把每一幀的内容合成出來，再按一定規則拼加在一起。其中疊加的意義在于使幀和幀之間過度自然。

每個幀内的每個峰值對應一段正弦波，其振幅、頻率恒定，相位由中心相位和頻率共同決定，也是用頻率逐點疊代。

友善起見，這裡直接寫代碼了，畢竟夠短。

leftHalf = floor(N/2);	  % N是幀内的采樣點數
startPha = pha - (leftHalf+1)*freq;	% pha是分析得到的相位，freq是分析得到的頻率，均視作中心相位、頻率。
phas = startPha + freq.*(1:N);
y = sum(amp.*cos(phas), 1);       % 所有正弦波相加
           

合成好了幀内的内容便可以進行疊加了。疊加前要進行加窗，這裡要保證窗疊加的結果是常數，這被稱為常數疊加(Constant Overlap-Add, COLA).

關于COLA，中文社群裡已經有很多介紹了，不再贅述。

4. 結語

本文主要介紹了Sinusoidal Modeling的分析、合成過程。讨論了實踐中的一些細節。

關于本文的程式設計實作，可以參考筆者的Gitee項目基于Matlab的正弦模型項目。到底能不能用三角函數表示聲音？運作代碼聽聽差別就知道啦。

下一篇文章我們将讨論正弦模型的一些實際應用場合。

歡迎讀者就文章及相關内容留言探讨。

注腳

1. Moorer, J. A. 1973. “The Hetrodyne Filter as a Tool for Analysis of Transient Waveforms.” MemoAIM-208, Stanford Artificial Intelligence Laboratory, Computer Science Dept., Stanford University. ↩︎

2. Portnoff, M.R. 1976. “Implementation of the Digital Phase Vocoder Using the Fast Fourier

   Transform.” IEEE Transactions on Acoustics, Speech and Signal Processing 24(3):243–248. ↩︎

3. Smith, J.O. and X. Serra. 1987. “PARSHL: An Analysis/Synthesis Program for Non-HarmonicSounds based on a Sinusoidal Representation.” Proceedings of the 1987 International Computer Music Conference. San Francisco: Computer Music Association. ↩︎
4. McAulay, R.J. and T.F. Quatieri. 1986. “Speech Analysis/Synthesis based on a Sinusoidal Representation.” IEEE Transactions on Acoustics, Speech and Signal Processing 34(4):744–754. ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ ↩︎
5. GABOR, D. Acoustical Quanta and the Theory of Hearing. Nature 159, 591–594 (1947). https://doi.org/10.1038/159591a0 ↩︎
6. http://ocw.nthu.edu.tw/ocw/index.php?page=course&cid=130& ↩︎ ↩︎ ↩︎
7. F.J. Harris, On the use of windows for harmonic analysis with the discrete Fourier transform, Proc. IEEE 66 (1) (1978) 51–83. ↩︎
8. Smith, Julius & Serra, Xavier. (2005). PARSHL: An Analysis/Synthesis Program for Non-Harmonic Sounds Based on a Sinusoidal Representation. Proceedings of the International Computer Music Conference. ↩︎ ↩︎
9. Serra, X. . Musical Sound Modeling with Sinusoids plus Noise. Musical Signal Processing. 1997. ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ ↩︎
10. M. Abe and J. Smith, “Design criteria for simple sinusoidal parameter estimation based on quadratic interpolation of FFT magnitude peaks”, (AES 2004) ↩︎
11. 未找到對應中文譯名，筆者譯。 ↩︎ ↩︎ ↩︎
12. George, E. Bryan, Smith, Mark J. T. An Analysis-by-Synthesis Approach to Sinusoidal Modeling Applied to the Analysis and Synthesis of Musical Tones[J]. journal of the audio engineering society, 1991. ↩︎ ↩︎
13. 聲音的泛音/諧波頻率通常并非基頻的完美倍數。例如本文圖6中，鋼琴C4音的頻譜的前五個諧波的頻率比是1.0000 : 2.0000 : 3.0033 : 4.0075 : 5.016314。導緻這種非整數的原因是琴弦僵硬的彎折22。 ↩︎ ↩︎
14. Szeto, Wai Man , and K. H. Wong . “Sinusoidal modeling for piano tones.” Signal Processing, Communication and Computing (ICSPCC), 2013 IEEE International Conference on IEEE, 2013. ↩︎ ↩︎
15. McAulay, R. J., & Quatieri, T. F. (1988). Speech Processing Based on a Sinusoidal Model. The Lincoln Laboratory Journal, 1(2), 153–168. ↩︎ ↩︎
16. Serra, X. 1989. “A System for Sound Analysis/Transformation/Synthesis Based on a Deterministic Plus Stochastic Decomposition.” PhD thesis, Stanford University. ↩︎ ↩︎
17. 注：如果一個噪聲在所有頻率成分上能量都相等，則是白噪聲(white noise)；否則為有色噪聲。相關内容參見百科詞條“有色噪聲”。 ↩︎
18. Verma, T. S. , & Meng, T. H. Y. . (2000). Extending spectral modeling synthesis with transient modeling synthesis. Computer Music Journal, 24(2), 47-59. ↩︎ ↩︎
19. Scott Nathan. Levine, Smith, & Julius Orion. (1998). Audio representations for data compression and compressed domain processing /. ↩︎
20. Rao, Kamisetty & Yip, P… (1990). Discrete cosine transform. Algorithms, advantages, applications. ↩︎
21. https://ccrma.stanford.edu/ ↩︎
22. A. Askenfelt, Ed., Five Lectures on the Accoutics of the Piano. Royal Swedish Academy of Music, 1990, available online at http://www.speech.kth.se/music/5_lectures/. ↩︎