基音頻率檢測

一、概念

何為基音周期？人在發音時，根據聲帶是否振動可以將語音信號分為清音和濁音兩種。濁音攜帶大量的能量，因此又被稱為有聲語音，其在時域上有明顯的周期性。而清音類似于白噪聲，沒有明顯的周期性。發濁音時，氣流通過聲門使聲帶產生張弛震蕩式振動，產生準周期的激勵脈沖串。這種聲帶振動的頻率稱為基音頻率；相應的周期就稱為基音周期。

基音頻率與個人聲帶的情況有關，包括聲帶長短、薄厚、韌性、勁度和發音習慣，總的來說基音頻率就是說話人的特征之一。而且基音頻率還隨著人的性別、年齡不同而有所不同。男性大概在70-200Hz，女性大概在200-450Hz之間。

二、檢測方法

盡管基音周期在目前有非常多的方法，但這些方法都具有局限性，沒有一種檢測方法能夠適用于不同的說話人、不同的要求環境，主要原因歸納為如下方面：

• 語音信號變化復雜，聲門激勵的波形并不是完全的周期脈沖串，語音的尾部也不具有聲帶振動的周期性，對有些清濁音的過渡幀很難判定其的周期性。
• 聲道共振峰有時會影響激勵信號的諧波結構。
• 在濁音語音段很難對每個基音周期的開始和結束位置進行精確的判斷
• 語音信號常常混有噪聲
• 基音頻率變化范圍大，從低音男聲的70Hz到兒童女性的450Hz，接近3個倍頻程給基音檢測帶來了一定的困難。

目前基音檢測算法大致可以分為兩大類：非基于事件檢測方法和基于事件檢測方法，事件指的是聲門閉合。

非基于事件的檢測方法主要有：自相關函數法、平均幅度差函數法、倒譜法等。非基于事件的檢測方法是利用語音信號短時平穩性這一特點，先將語音信號分為長度一定的語音幀，然后對每一幀語音求基音周期。它的優點是：算法簡單，運算量小，但缺點在于：無法檢測幀內基音周期的非平穩變化，檢測精度不高。

基于事件的檢測方法主要有：小波變換、Hilbert-Huang變換。基于事件的檢測方法是通過聲門閉合時刻來對基音周期進行估計，而不需要對語音信號進行短時平穩假設。優點是：在時域和頻域上有良好的局部特性，能跟蹤基音周期的變化，并能將微小的基音周期變化檢測出來，檢測精度高。缺點是：計算量較大

三、估計一幀信號的基音頻率

倒譜法

由于語音$x(i)$是由聲門脈沖激勵$u(i)$經聲道響應$v(i)$濾波而得，即
$$x(i)=u(i)?u(i)$$
設這三個量的倒譜分別為$x(i)、u(i)、v(i)$,則有：
$$x(i)=u(i)?v(i)$$
由于在倒譜域中$u(i)和v(i)$是相對分離的，說明包含有基音信息的聲脈沖倒譜可以與聲道響應倒譜分離，因此從倒譜域分離$u(i)和u(i)$。在計算出倒譜后，就在倒譜頻率為$P_{min}～P_{max}$之間尋找倒譜函數的最大值，倒譜函數最大值對應的樣本的點數就是當前幀語音信號的基音周期$T_0(n)$，基音頻率為$F_0(N)=f_s/T_0(n)$

從matlab中取出一幀信號，進行256點的fft變換，再對其求模并取對數后得到幅度譜，取兩個相鄰的峰值點，如下圖所示

clc; close all; clear all; [x1,Fs] = audioread('voice/a9.wav'); wlen=256; inc=128; % 給出幀長和幀移 N=length(x1); % 信號長度 time=(0:N-1)/Fs; % 計算出信號的時間刻度 signal=enframe(x1,wlen,inc)'; % 分幀 framedata = signal(:,15); x2 = fft(framedata); amp = 20*log10(abs(x2)+0.0000001); plot(amp)

從圖中可以看出兩個峰值點的周期相差9，根據FFT的性質，相鄰點頻率間隔是：
$$\delta =\frac{f_s}{N}=\frac{8000}{256}=31.25Hz$$
因此相差9點周期對應的基頻就是：
$$f_b=9\delta =281.25Hz$$
由于FFT變換后兩邊是對稱的，因此用matlab中fftshift（）將對稱一邊的頻譜搬移過來，再取點進行計算，不難發現對應的基頻也在281.25Hz。

signal=enframe(x1,wlen,inc)'; % 分幀 framedata = signal(:,15); x2 = fft(framedata); amp = 20*log10(abs(x2)+0.0000001); plot(fftshift(amp));

如果用多個間隔來估計相鄰點頻率檢測，那就是
$$f_b=\frac{65?20}{5}\delta =281.25Hz$$
如果將這個幅度譜看做是時域信號，則其中周期性對應時間周期是：
$$t=9\frac{1}{f_s}=0.001125s$$
所對應的頻率就是：
$$f=\frac{f_s}{9}\approx 888.9Hz$$
對幅度譜進行fft變換，獲得的新的頻域幅度上對應888.9Hz的點數就是
$$\frac{888.9}{\delta }=\frac{888.9}{31.25}=28.4448\approx 28點$$

framedata = signal(:,15); x2 = fft(framedata); amp = 20*log10(abs(x2)+0.0000001); % plot(fftshift(amp)); x3 = fft(fftshift(amp)); amp2 = abs(x3); plot(amp2);

從上圖中不難看出，他們之間在29~116之間均有峰值，并且間隔大約在28。

過程推導：
$$\begin{gathered} \frac{f_s}{m}?\frac{f_s}{N}=\frac{N}{m}=\frac{888.9}{31.25}=28點 \\ {f}_b=m\frac{f_s}{N}=281.25Hz \end{gathered}$$

所以
$$f_b=\frac{f_s}{28.4448}=281.25Hz$$
想法：

要求一幀信號中的基頻，假設如下幾個方法：

• 求一幀信號對數幅度譜的FFT譜，取一半之后，找幾個峰值點的位置，求出峰值點的間隔

其他軟件驗證

將這個語音放進Adobe audition中，從語譜圖來看，第一根能量柱大約是在300Hz往下一點點。

從頻譜圖來看，第一個峰值所對應的頻率是281Hz，與上面求出來的數值誤差不大

算法實現

根據上面的想法，自己編寫了一個程序來計算基頻，大致思路就是多個間隔來估計相鄰點頻率檢測，取兩次FFT變換的幅度譜區間的一半，找出這個區間內的極大值點，根據點的個數算出基音頻率。

clc; close all; clear all; [x1,Fs] = audioread('voice/a9.wav'); wlen=256; inc=128; % 給出幀長和幀移 N=length(x1); % 信號長度 time=(0:N-1)/Fs; % 計算出信號的時間刻度 signal=enframe(x1,wlen,inc)'; % 分幀 framedata = signal(:,15); x2 = fft(framedata); amp = 20*log10(abs(x2)); x3 = abs(fft(amp)); x4 = x3(1:128); %找出極大值點 datacursormode on [y,x] = findpeaks(x4,'MinPeakHeight',194); plot(x4); % axis([0,128,0,600]); text(x+.02,y,num2str((1:numel(y))')) hold on; plot(x,y,'g.');%計算fs x_len = length(x); x_min = min(x); x_max = max(x); fs = (x_max - x_min)/(x_len - 1)*31.25

運行結果如下圖所示，算出的基頻是437.5Hz，與Au中找出的181Hz有較大差距，從找出極大值的點來看，顯然還是點沒找對，比如1點，3點應該是排除的，但由于這些點并不是均勻從大到小的，如果將閾值設置的比1還高，那么4,5,6三個點就會被過濾，因此采取設置閾值來尋找極大值的方法不是十分合適的。

最后根據老師的思路，用話音的基音頻率的范圍來確定x軸范圍中的最大值

clc; close all; clear all; [x1,Fs] = audioread('voice/a9.wav'); wlen=256; inc=128; % 給出幀長和幀移 N=length(x1); % 信號長度 time=(0:N-1)/Fs; % 計算出信號的時間刻度 signal=enframe(x1,wlen,inc)'; % 分幀 framedata = signal(:,15); x2 = fft(framedata); amp = 20*log10(abs(x2)); x3 = abs(fft(amp)); x3(1:27) = 0;%在話音基頻范圍外的都取零 x3(115:256) = 0; [M,idx] = max(x3); f_b=Fs/(idx-1); plot([0:255],x3); hold on; plot(idx,M,'g.');

所求出的基音頻率是285.7Hz，與我們剛開始人工求出的281.25Hz還是十分接近的。

短時自相關法

短時自相關法基音檢測主要是利用短時自相關函數的性質，通過比較原始信號及其延遲后信號間的類似性來確定基音周期。歸一化自相關函數的最大幅值是其他延遲量時，幅值都小于1.如果延遲量等于基音周期，那兩個信號具有最大類似性，直接找出短時自相關函數的兩個最大值間的距離，就可以作為基音周期的初估值。再求出短時自相關后和倒譜法尋找最大值一樣，用相關函數法時也在$P_{min}～P_{max}$見尋找歸一化相關函數的最大值，最大值對應的延遲量就是基音周期。

還是用原來的元音a的音頻，取出一幀來進行短時自相關計算

clc; close all; clear all; [x1,Fs] = audioread('voice/a9.wav'); %plot(x1); lag=255; %觀察一段范圍的元音 wlen=256; inc=128; % 給出幀長和幀移 N=length(x1); % 信號長度 time=(0:N-1)/Fs; % 計算出信號的時間刻度 signal=enframe(x1,wlen,inc)'; % 分幀 %取出一幀數據 % frametime = (frame-1)*inc+1:(frame-1)*inc+wlen; framedata = signal(:,15); % plot(frametime,framedata); %對數據進行短時自相關計算 [c2,lags2,bound]=autocorr(framedata,lag); plot(lags2,c2);

人工計算

取x軸為12,56,113,141的四個點進行計算
$$\begin{gathered} f_{b1}=\frac{8000}{56?28}=285.71Hz \\ {f}_{b2}=\frac{8000}{141?113}=285.71Hz \end{gathered}$$
機器計算

與倒譜法中用到的利用話音的基音頻率的范圍來確定x軸范圍中的最大值來計算基音頻率

clc; close all; clear all; [x1,Fs] = audioread('voice/a9.wav'); %plot(x1); lag=255; %觀察一段范圍的元音 wlen=256; inc=128; % 給出幀長和幀移 N=length(x1); % 信號長度 time=(0:N-1)/Fs; % 計算出信號的時間刻度 signal=enframe(x1,wlen,inc)'; % 分幀 %取出一幀數據 % frametime = (frame-1)*inc+1:(frame-1)*inc+wlen; framedata = signal(:,15); % plot(frametime,framedata); %對數據進行短時自相關計算 [c2,lags2,bound]=autocorr(framedata,lag); % plot(lags2,c2); %在話音基頻范圍外的都取零 c2(1:27) = 0; c2(115:256) = 0; [M,idx] = max(c2); f_b=Fs/(idx-1); plot([0:255],c2); hold on; plot(idx,M,'g.');

算出的基頻是285.7Hz，與手工算出的基頻是一致的

短時幅度差

短時平均幅度差函數：
$$F_n(k)=\sum _{m=0}^{N?1?k}\left|x_n(m)?x_n(m+k)\right|,0<k\le K$$

clc; close all; clear all; [b,Fs] = audioread('voice/a9.wav'); b1=b(3000:3256); % 選取一段語音 N=128; % 窗口大小 A=[]; for k=1:128; % 延遲長度sum=0;for m=1:N;sum=sum+abs(b1(m)-b1(m+k-1));endA(k)=sum; end % A(1:27) = 100;%在話音基頻范圍外的都取零 % A(115:256) = 100; % [M,idx] = min(A); % f_b=Fs/(idx-1); % s=b(1:320) figure(1) subplot(2,1,1) plot(b1); xlabel('樣點') ylabel('幅度') % axis([0,320,-0.2,0.2]) subplot(2,1,2) plot(A); xlabel('延時k') ylabel('AMDF') % axis([0,320,0,60]);

計算x=57和x=29兩點間的基音頻率
$$f_b=\frac{8000}{57?29}=285.71Hz$$

四、估計一段語音的基音頻率

將分幀之后的每一幀數據都進行基頻估計，得到該語言每一幀的基音頻率圖。

倒譜法

%求出一段語音的基音頻率 clc; close all; clear all; [x1,Fs] = audioread('voice/a9.wav'); subplot(2,1,1); plot(x1); wlen=256; inc=128; % 給出幀長和幀移 N=length(x1); % 信號長度 time=(0:N-1)/Fs; % 計算出信號的時間刻度 signal=enframe(x1,wlen,inc)'; % 分幀 [n,m]=size(signal); for i=1:mframedata = signal(:,i); % f_b=pitch_cep(framedata,Fs); %倒譜法 % f_b=pitch_cor(framedata,Fs); %自相關法f_b=pitch_admf(framedata,Fs); %平均幅度差x(i)=f_b; end subplot(2,1,2); plot(x);

自相關法

短時幅度差法

畫出的出現了基頻為8000的時候，究其原因是找兩個最小值的點重合了

查看第一幀的圖，由于除法公式可以知道，當分母趨于0的時候，結果就會接近最大值

五、利用voicebox驗證

還沒有辦法做到將基頻疊加到語譜圖中，初步猜想可能需要插值

總結

以上是生活随笔為你收集整理的基于倒谱法、自相关法、短时幅度差法的基音频率估计算法（MATLAB及验证）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。