Mathematica 歌姬計劃（2）—— 聽力訓練

02-13

歡迎來到 Mathematica 歌姬計劃，這一計劃的目標是用 MMA 像 Vocaloid 那樣調教自己想要的人聲電子音樂。

前情提要

在上篇文章中，我們研究了如何在 Mathematica 中調用 Cortana 音源：

yxlllc：Mathematica 歌姬計劃（1）—— 與 Cortana 的聯合zhuanlan.zhihu.com

本文為本系列的第二篇文章，在前一篇的基礎上，介紹如何運用 Mathematica 測量並獲取 Cortana 音源的基頻與共振峰信息。 2.1 節介紹將 Speak 函數產生的音頻流記錄為波形的方法， 2.2 節介紹基頻曲線（PIT）的測量方法， 2.3 節介紹正弦模型以及共振曲線的測量方法，以及最後總結。

由於本系列的核心重點是MMA編程，本文假設讀者對 DSP (即 Digital signal processing，數字信號處理) 以及物理聲學相關的基礎知識有所了解，將不會涉及理論部分的一系列推導。本文會盡量避免數學公式的出現，將重點放在編程實現的過程上面。

2.1 錄下 Cortana 的聲音

在上篇文章最後，我們新定義了 MySpeak 函數，使得 Cortana 可以成功說中文，解決了MMA 自帶的 Speak 函數不支持中文輸入的問題。然而，不論是原來的 Speak 函數還是新定義的 MySpeak 函數，它們的輸出結果均是音頻流的形式。也就是說，這種情況下Cortana 姐姐的聲音是轉瞬即逝的，就像英語聽力考試一樣，說完了，就沒了，然後麥醬一臉懵逼。

考慮到麥醬的計算力有限（相比於C語言來說），要做歌聲合成的話，基於流的音頻實時處理是一件幾乎不可能完成的事情。所以一個直接的想法就是把音頻流產生的聲音「存起來」，變成波形數據的形式，可以被程序實時調用並做下一步計算。

一個最樸素的解決方案是，把音頻流直接錄製下來。"一遍沒聽懂就錄下來回家慢慢聽～"，於是麥醬從肚兜里掏出了 AudioCapture 函數（11.2 版本最新函數）：

展開上面畫圈的箭頭：

點擊這個「話筒」圖標可以設置音頻輸入設備：

我的電腦目前只能支持外錄，也就是利用麥克風硬體錄製環境中的聲音。有些電腦是支持內錄的，即定位音效卡上臨時存放音頻流波形信息的緩衝區，在播放時不斷將它們提取出來，相當於做一個「音頻流 ->文件流"的操作。就算電腦默認不支持內錄，我們也可以學習各種屏幕錄製軟體所採用的方法，配置虛擬音效卡欺騙操作系統，假裝外錄，實際上是內錄。這裡就不介紹具體怎麼配置虛擬音效卡了，總之配置好了之後，剛才的音頻輸入設備中會多出專門用來內錄的"虛擬設備"。選擇這個設備，然後在Cortana 發聲時點擊紅色按鈕進行錄製就行了。

正當麥醬調好錄音機準備開始給 Cortana 姐姐錄音時，Cortana 姐姐表示："不用麻煩你翻錄了，我這有現成的錄音帶"。

原來，微軟校長早已考慮到像麥醬這樣的學生對聽力素材有著巨大的需求，於是在設計Cortana 這樣的 TTS 平台之初就提供了文件流輸出的功能。讓我們回顧一下上篇文章提到的微軟官方說明文檔：

SpeechSynthesizer 類msdn.microsoft.com

注意看這一行：

也就是說我們直接調用這個函數就可以生成 Cortana 聲音的波形文件，直接拿到"錄音帶"！

於是我們稍微修改一下上次的 MySpeak 函數成為 VoiceData 函數

Needs["NETLink`"];LoadNETAssembly["System.Speech"];Options[VoiceData] = {"Rate" -> 0, "Volume" -> 100, "SampleRate" -> 44100, "SampleDepth" -> 16, "Channels" -> 1};VoiceData[string_, OptionsPattern[]] := Module[{synth, format, tmpfile, data}, synth = NETNew["System.Speech.Synthesis.SpeechSynthesizer"]; synth@Rate = OptionValue["Rate"]; synth@Volume = OptionValue["Volume"]; format = NETNew["System.Speech.AudioFormat.SpeechAudioFormatInfo", OptionValue["SampleRate"], OptionValue["SampleDepth"], OptionValue["Channels"]]; tmpfile = Close@OpenWrite[] <> ".wav"; synth@SetOutputToWaveFile[tmpfile, format]; synth@Speak[ToString[string]]; synth@Dispose[]; data = First@AudioData@AudioTrim[Import[tmpfile]]; DeleteFile[tmpfile]; SampledSoundList[data, OptionValue["SampleRate"]] ]

大概思路就是利用 NetLink 調用上述 SetOutputToWaveFile 方法指定 Cortana 輸出音頻文件的地址 (Close@OpenWrite[] <> ".wav"），作為一個臨時文件 (tmpfile) ，然後用 Import 函數導入這個 WAV 波形至 MMA 中，成功導入後將它刪掉。同樣像上次的 MySpeak 函數那樣，我們可以分別用 "Rate" 和 "Volume" 選項設置發音速度和音量，另外增加的選項是 "SampleRate" ，可用來指定採樣率，比如輸入以下代碼：

ssl = VoiceData["你好", "Rate" -> 0, "Volume" -> 100, "SampleRate" -> 44100]

運行後變數 ssl 被賦值為一個 SampledSoundList 對象：

在 MMA 中，SampledSoundList 對象可以被 Sound 函數作用，生成音頻播放界面：

Sound[ssl]

點擊左下角的三角形播放鍵就可以聽到 Cortana 姐姐說「你好」了。

2.2 抑揚頓挫

聲音的三要素是音調、音色和響度，人類語音自然也不例外。雖然麥醬現在獲得了錄音帶，直接聽還是聽不懂 Cortana 姐姐在說什麼，哪怕是漢語中最簡單的「你好」，於是開始思考如何把三要素信息從這些數據中提取出來。

我們首先考慮音調。無論什麼語言，「抑揚頓挫」都是體現說話者語氣和情緒的關鍵。在漢語中，語調的重要性相對更高，因為不同的「平仄」幾乎直接對應不同的漢字，表示完全不同的信息。在物理中，音調對應的是頻率，更準確的說，是基本頻率，簡稱基頻 (fundamental frequency)。在 Vocaloid 中，這便是 PIT 參數。

在 MMA 中，AudioLocalMeasurements 函數為這種測量提供了可能，首先打開幫助頁面：

展開"更多信息和選項" ，找到"頻率方面屬性"：

果然提供了基頻的測量選項，再往下還有個說明：

其中參數 t 的含義比較令人費解，經過研究，它體現的是靈敏度，值越小對快速變化的信號等越不敏感，一般取 0.5 左右比較合適。對於 Cortana 這樣的女聲音源，minfreq 一般取 100 左右，maxfreq 一般取 400 左右。

於是我們把這些經驗性的參數都封裝一下，得到如下的 FBase 函數：

Options[FBase] = {"Range" -> {0.4, 110, 440}, "Partition" -> {1024/44100, 512/44100, HannWindow}};FBase[ssl_, OptionsPattern[]] := AudioLocalMeasurements[ssl, Flatten[{"FundamentalFrequency", OptionValue["Range"]}], PartitionGranularity -> OptionValue["Partition"], MissingDataMethod -> {"Interpolation", InterpolationOrder -> 1}]

其中 "Range" 選項的三個參數分別依次對應上述的 t , minfreq 和 maxfreq. "Partition" 參數表示對音頻的劃分和加窗，經過大量實驗，它們的默認值基本可以達到最好的效果，因此調用時可不做任何修改。比如對之前那段"你好"（變數ssl）, 我們直接調用函數：

fbase=FBase[ssl]

運行後變數 fbase 被賦值一個 TimeSeries 對象：

在 MMA 中，TimeSeries 對象可以被 ListPlot 或 ListLinePlot 等函數直接作用，直接生成圖像，這便是基頻 (Hz) - 時間 (s) 曲線（Vocaloid里的PIT曲線）

ListLinePlot[fbase]

「你好"（Cortana）的 PIT 曲線

從上圖可以明顯看出一個雙「3聲」結構，跟漢語拼音長得是一模一樣。注意到「你」字「3聲」開頭的下降段非常短促乃至非常接近「2聲」，這是漢語口語的連讀習慣。

2.3 聽音識人識字

音調考慮完了，我們來考慮音色。對於人類語言來說，音色有兩個最大的作用：

區分不同的講話者
區分不同的音素

其中第2條直接決定我們為什麼能聽到」你好「。不同的母音或輔音，在頻域上各個頻段的能量分配差別是很大的。在激勵源（比如聲帶）不變的情況下，這種能量分配主要由發聲腔體的物理結構決定（比如是張大嘴還是抿著嘴），可以近似認為與激勵源本身的頻率（語調）無關。

現在我們不得不用數學來精確描述這個模型了。設 t 時刻激勵源的頻譜係數為關於頻率 f 的兩個函數 A(t,f) 和 $phi(t,f)$ ，分別表示振幅和相位。再設 t 時刻時腔體結構對頻率 f 的放大作用（共振）為 g(t ,f), 最終的波形是 a(t)，則可以認為：

$a(t)=int_{0}^{infty } g(t,f)A(t,f)cos[2pi ft+phi(t,f)]df$

相當於認為跟腔體共振有關的 g 與跟激勵源有關的 A 和 $phi$ 是獨立的，可以分別寫開。這個近似對於母音是非常好的，因為根據語音學的定義，母音在發音過程中產生的氣流不會摩擦腔體，意味著激勵源與腔體之間不會產生無規則的相互作用（即輔音中的雜訊成分，比如齒音）。

由於任一小段時間內激勵源產生的是有周期的信號，傅里葉變換的結果 A 應該是離散的，也就是說：

$A(t,f)=sum_{k=1}^{infty}{a_{k}(t)delta(f-k f_{base}(t))}$

這裡的 $f_{base}$ 就是上一節測出來的 PIT 曲線， $a_{1},a_{2}...a_{n}$ 描述的是離散化後的傅里葉級數中不同倍頻的振幅大小。

將它代入原式，在狄拉克 $delta$ 函數的作用下連續的積分會被轉化為離散求和：

$a(t)=sum_{k=0}^{infty}{g(t,k f_{base}(t))a_k(t)cos[2pi f_{base}(t)t+phi(t,kf_{base}(t))]}$

激勵信號一般具有很強的隨機性。為了進一步簡化模型，我們假設激勵源在不同倍頻處的平均能量相等（類似於白雜訊），即給定 t 時刻不同的 $a_k$ 都當成一個定值，這樣我們可以直接把所有的 $a_k$ 都可以吸收至前面的 g 函數中，記為 G 函數：

$a(t)=sum_{k=0}^{infty}{G(t,k f_{base}(t))cos[2pi f_{base}(t)t+phi(t,kf_{base}(t))]}$

這便是語音的正弦模型。給定 t 時刻，關於 f 的函數 G(t, f) 稱為該時刻的共振曲線。

在 MMA 中，並沒有測量共振曲線相關的函數，我們只能自己構造了。

WindowList[windowfunction_, n_] := WindowList[windowfunction, n] = Array[windowfunction, n, {-0.5, 0.5}];LSFKernel[scanf_, n_, samplerate_] := LSFKernel[scanf, n, samplerate] = LeastSquaresFilterKernel[{"Bandpass", (2 [Pi] scanf)/samplerate}, n];Formant[ssl_, fbase_, nwidth_: 2049, nshift_: 512, windowfunction_: HannWindow, fwidth_: 220, fshift_: 20, flimit_: 8000] := Module[{datalist, wdatalist, amplist, tlist, LineFit, func}, LineFit = Interpolation[#, InterpolationOrder -> 1] &; datalist = Partition[ssl[[1]], nwidth, nshift, {-nwidth, 1}, 0]; wdatalist = WindowList[windowfunction, nwidth]*# & /@ datalist; amplist = MapThread[ Sqrt[2/nwidth*Total[#1^2]/Total[#2^2]] &, {datalist, wdatalist}]; tlist = (nshift (Range@Length[datalist] - 1) + nwidth/2)/ssl[[2]]; func = LineFit@MapThread[{#1, LineFit@ Table[{f, #2 Total[ ListConvolve[ LSFKernel[ If[f == 0, {0, fwidth/2}, fshift*Floor[(f - fwidth/2)/fshift + 1/2] + {0, fwidth}], nwidth, ssl[[2]]], #3, Floor[(nwidth + 1)/2], 0]^2]^0.5}, {f, 0, flimit, Quiet@fbase[#1]}]} &, {tlist, amplist, wdatalist}]; Quiet[Expand[func[#1]][#2] /. (a_*func1_ + b_*func2_)[#2] :> a*func1[#2] + b*func2[#2]] &];

構造出的 Formant 函數將直接計算得到的正弦模型中的 G，大概思路是首先對音頻輸入信號劃分（Partition函數）並加窗（乘上 windowfunction）, 然後根據前面計算過的基頻利用最小方均 FIR 濾波器（LeastSquaresFilterKernel 函數）分離各個頻段的信號（帶通濾波），計算能量值並擬合出共振曲線。

你可能注意到 Formant 函數有非常多的參數，不過跟上一節一樣，經過大量實驗，它們的默認值基本可以達到最好的效果，因此調用時可不做任何修改。因此只需要前兩個變數作為輸入，分別是對應語音信號 SampledSoundList 對象（ssl）和一個 PIT 曲線 (fbase)，比如我們將「你好」的 ssl 和 fbase 作為輸入：

G = Formant[ssl, fbase];

然後我們可以創建個動態來觀察共振曲線隨時間的變化情況：

Manipulate[ Plot[20 Log10@G[t, f], {f, 0, 8000}, PlotRange -> {{0, 8000}, {-100, 0}}], {t, 0, 0.57}]

為了更好的觀察，建議採用的分貝坐標來觀察曲線的演化，即做一個 20*Log10 操作，這裡就不上傳 GIF 結果了，只給一個截圖作為舉例：

母音「ao」(Cortana) 的共振曲線

上圖取自「你好」音頻的第 0.423 秒，此時 Cortana 正在發「ao」這個母音，可以看出有3個共振峰，分別在 1000 Hz, 4000 Hz 和 7000 Hz 左右。

總結

經過大量的聽力訓練，麥醬現在已經大體上熟悉了 Cortana 姐姐的聲音（建模完成），終於能聽懂講話的內容了！這為之後的首次歌唱做出了關鍵鋪墊。