作者：桂。
時間：2017-05-05 21:45:07
鏈接：http://www.cnblogs.com/xingshansi/p/6815217.html
前言
主要總結一下常用的音頻特征，並給出具體的理論分析及代碼。
一、過零率
過零率的表達式為：
其中N為一幀的長度，n為對應的幀數，按幀處理。
理論分析：過零率體現的是信號過零點的次數，體現的是頻率特性。因為需要過零點，所以信號處理之前需要中心化處理。
code(zcr1即為過零率):
for i=1:fn
z=X(:,i); % 取得一幀數據
for j=1: (wlen- 1) ; % 在一幀內尋找過零點
if z(j)* z(j+1)< 0 % 判斷是否為過零點
zcr1(i)=zcr1(i)+1; % 是過零點，記錄1次
end
end
end
二、短時能量
短時能量的表達式為：
理論分析：短時能量體現的是信號在不同時刻的強弱程度。
code:
for i=1 : fn
u=X(:,i); % 取出一幀
u2=u.*u; % 求出能量
En(i)=sum(u2); % 對一幀累加求和
end
三、短時自相關函數
短時自相關函數定義式為：
理論分析：學過信號處理的都應該知道，信號A與信號B翻轉的卷積，就是二者的相關函數。其中
是因為分幀的時候，加了窗函數截斷，w代表窗函數。
code:假設一幀截斷的信號
r = xcorr(signal);
這與直接利用卷積的方式等價：
給出卷積的實現：
function output_signal=my_direct_convolution(input_signal,impulse_response)
% Input:
% input_signal: the input signal
% impulse_response: the impulse response
% Output:
% output_signal:the convolution result
N=length(input_signal);%define length of signal
K=length(impulse_response);%define length of impulse response
output_signal=zeros(N+K-1,1);%initializing the output vector
xp=[zeros(K-1,1);input_signal;zeros(K-1,1)];
for i=1:N+K-1
output_signal(i)=xp(i+K-1:-1:i)'*impulse_response;
end
卷積也可以借助FFT快速實現。
調用卷積的函數：
r1 = my_direct_convolution(signal,signal(end:-1:1));
圖中可以看出r與r1完全等價：
四、短時平均幅度差
假設x是加窗截斷后的信號，短時平均幅度差定義：
理論分析：音頻具有周期特性，平穩噪聲情況下利用短時平均幅度差可以更好地觀察周期特性。
code:
取一幀信號，計算短時平均幅度差：
u = X(:,i) %取一幀信號
for k = 1:wlen
amdvec(k) = sum(abs(u(k:end)-u(1:end-k+1)));%求每個樣點的幅度差再累加
end
前面四個都是信號的時域分析，音頻信號更多是在時頻域分析(可借助tftb-0.2工具包分析)。
常用的有STFT（短時傅里葉變換）、小波變換、ST、W-V變換，以線性調頻信號為例：
左圖最下面為合成信號，右圖為四種變換對合成信號進行的時頻分析。
這里只分析利用短時傅里葉變換（Short time fourier transform, STFT）的情形。
又因為實數的傅里葉變換共軛對稱，有時也僅僅分析頻域的一半信息即可。
為什么要進行STFT呢，原因按我的理解可能有兩點：
- 傳統FFT只能看到信號頻率的特性，時域信號只能觀察時域特性，都是一維的情況，如果二維聯合觀察？這個時候STFT就可以實現;
- 語音是非平穩信號，比如求相關矩陣，理論上是E{.}求取均值的形式，通常無法得出概率密度，往往有數據近似：
這個式子能夠近似相關矩陣，有兩個前提條件：a)信號平穩，這樣才能保證統計特性一致；b)遍歷性，這個時候才能保證統計沒有以偏概全。
但語音信號是非平穩信號，直接求取相關矩陣理論上沒有意義，其他統計信息也有類似的特性。但語音變化緩慢，可以認為是短時平穩，即在短的時間內（如20~30ms）是平穩的，這個時候平穩+遍歷性的假設，就可以讓我們借助觀測數據估計統計信息。這個短時平穩的划分就是信號分幀。進一步：分幀信號分別FFT，就是STFT。
信號分幀的code:
Nframe = floor( (length(x) - wlen) / nstep) + 1; for k = 1:Nframe idx = (1:wlen) + (k-1) * nstep; x_frame = x(idx); end
加窗截斷分幀的示意圖：
為什么要有幀移量？即幀與幀之間有部分重疊？從上面中間圖可以看出，加窗階段后相鄰兩幀在端點處變化較大，對於變化較大的情況一般思路就是平滑，比如進行插值處理，其實幀移的操作就是插值呀。
五、語譜圖（基於FFT）
有時FFT也換成DCT實現，FFT延展與DCT是等價的，就不一一列出了。只分析FFT情況。
基於FFT語譜圖的定義：
就是分幀，對每一幀信號FFT
然后求絕對值/平方。
理論分析：
給出示意圖,語音分幀→每一幀分別FFT→求取FFT之后的幅度/能量，這些數值都是正值，類似圖像的像素點，顯示出來就是語譜圖。
語譜圖code:
function d=FrequencyCal(x,nw,ni)
n=nw; %幀長
h=ni; %幀移量
s0=length(x);
win=hamming(n)'; %加窗,hamming為例
c=1;
ncols=1+fix((s0-n)/h); %分幀，並計算幀數
d=zeros((1+n/2),ncols);
for b=0:h:(s0-n)
u=win.*x((b+1):(b+n));
t=fft(u);
d(:,c)=t(1:(1+n/2))';
c=c+1;
end
讀取語音調用語譜圖code:
[signal,fsc] = wavread('tone4.wav');
nw=512;ni=nw/4;
d=FrequencyCal(signal',nw,ni);
tt=[0:ni:(length(signal)-nw)]/fsc;
ff=[0:(nw/2)]*fsc/nw*2;
% imagesc(tt,ff,20*log10(abs(d)));
imagesc(tt,ff,abs(d).^2);
xlabel('Time(s)');
ylabel('Frequency(Hz)')
title('語譜圖')
axis xy
常用對數譜,修改上面的一句code即可：
imagesc(tt,ff,20*log10(abs(d))); % imagesc(tt,ff,abs(d).^2);
效果圖：
有時對數中添加常數項；
imagesc(tt,ff,20*log10(C+abs(d)));
效果圖：
六、短時功率譜密度
先來看看功率譜定義：
可見對於有限的信號，功率譜之所以可以估計，是基於兩點假設：1）信號平穩; 2)隨機信號具有遍歷性
A-功率譜密度
1）周期圖法
已知N個采樣點的信號
，對其進行傅里葉變換：
進一步得到功率譜密度：
按照上文分析相關函數的思路，給出一個分析：
就是信號u的相關函數就是u卷積上u的翻轉，而相關函數與功率譜密度是互為傅里葉變換。u對應傅里葉變換U，u的翻轉對應U的共軛，時域的卷積對應頻域的相乘，就得到了功率譜估計的表達式，同樣代碼實現依然可以借助上文分析相關函數的特性加以分析。
對應code，其中my_direct_convolution仍然調用上面的函數:
fs = 1000;%采樣率
f0 = 30;%信號頻率
t = 0:1/fs:1;
x = sin(2*pi*f0*t);
lenx = length(x)
subplot 121
plot(t,x,'k');
title('時域信號')
subplot 122
Prr = abs(fft(my_direct_convolution(x',x(end:-1:1)')))/lenx;
plot((0:fs)/2,Prr(1:length(Pxx)),'r--');
xlabel('Frequency(Hz)')
效果圖中可以看出信號30Hz可以明顯從功率譜圖觀察：
可以看出：功率譜密度與相關函數對應，而相關函數是統計信息，按前文提到的，它是建立在信號平穩的假設之上。如果信號不夠平穩呢？周期圖法的思路顯然是不適用的，Welch就是對這一問題的改進。
2）平均周期圖法（Welch）
與周期圖譜求功率譜密度不同，Welch不再從整個時間段考慮，而是做了三點改進：
- 截斷，將這個信號分成多個片段
- 加窗：因為截斷，截斷就要泄露，通常都選擇加窗處理
- 重疊：截斷之后，為了防止相鄰兩段差異過大，通常插值平滑處理，也就是取重疊
給出Welch定義，每一段的功率：
其中
，d2(n)是窗函數，M是每一段的長度，假設總長度為L，則Welch平均功率譜密度為：
簡單來理解就是：分段的每一段用周期圖法得到功率譜密度，然后加權平均，不再細說了。
B-短時功率譜密度
按前面分析，周期圖法針對的是平穩信號，而Welch雖然考慮了非平穩的特性，但分段的數量通常較小，每一段的長度較大，對音頻信號而言，這也是不夠的。音頻信號可以看作短時間的近似平穩（如一幀信號），對每一幀利用周期圖法分析，這個就是短時功率譜密度的思路。通常為了防止一幀的信號不夠平穩，每一幀也可以進一步利用Welch方法處理。
對應code，分幀、Welch都是上面的思路，直接調用了：
function [Pxx] = pwelch_2(x, nwind, noverlap, w_nwind, w_noverlap, nfft)
% 計算短時功率譜密度函數
% x是信號，nwind是每幀長度，noverlap是每幀重疊的樣點數
% w_nwind是每段的窗函數，或相應的段長，
% w_noverlap是每段之間的重疊的樣點數，nfft是FFT的長度
x=x(:);
inc=nwind-noverlap; % 計算幀移
X=enframe(x,nwind,inc)'; % 分幀
frameNum=size(X,2); % 計算幀數
%用pwelch函數對每幀計算功率譜密度函數
for k=1 : frameNum
Pxx(:,k)=pwelch(X(:,k),w_nwind,w_noverlap,nfft);
end
七、譜熵
譜熵的定義，首先對每一幀信號的頻譜絕對值歸一化：
這樣就得到了概率密度，進而求取熵：
理論分析：熵體現的是不確定性，例如拋骰子一無所知，每一面的概率都是1/6，信息量最大，也就是熵最大。如果知道商家做了手腳，拋出3的概率大，這個時候我們已經有一定的信息量，拋骰子本身的信息量就少了，熵也就變小。對於信號，如果是白噪聲，頻譜類似均勻分布，熵就大一些；如果是語音信號，分布不均勻，熵就小一些，利用這個性質也可以得到一個粗糙的VAD（有話幀檢測）。譜熵有許多的改進思路，濾波取特定頻段、設定概率密度上限、子帶平滑譜熵，自帶平滑通常利用拉格朗日平滑因子，這是因為防止某一段子帶沒有信號，這個時候的概率密度就沒有意義了，這個思路在利用統計信息估計概率密度時經常用到，比如朴素貝葉斯就用到這個思路。
譜熵code:
for i=1:fn
Sp = abs(fft(y(:,i))); % FFT變換取幅值
Sp = Sp(1:wlen/2+1); % 只取正頻率部分
Ep=Sp.*Sp; % 求出能量
prob = Ep/(sum(Ep)); % 計算每條譜線的概率密度
H(i) = -sum(prob.*log(prob+eps)); % 計算譜熵
end
八、基頻
基頻：也就是基頻周期。人在發音時，聲帶振動產生濁音(voiced)，沒有聲帶振動產生清音（Unvoiced）。濁音的發音過程是：來自肺部的氣流沖擊聲門，造成聲門的一張一合，形成一系列准周期的氣流脈沖，經過聲道（含口腔、鼻腔）的諧振及唇齒的輻射形成最終的語音信號。故濁音波形呈現一定的准周期性。所謂基音周期，就是對這種准周期而言的，它反映了聲門相鄰兩次開閉之間的時間間隔或開閉的頻率。常用的發聲模型：
基音提取常用的方法有：倒譜法、短時自相關法、短時平均幅度差法、LPC法，這里借用上面的短時自相關法，說一說基頻提取思路。
自相關函數：
通常進行歸一化處理，因為r(0)最大，
得到歸一化相關函數時候，歸一化的相關函數第一個峰值為k = 0，第二個峰值理論上應該對應基頻的位置，因為自相關函數對稱，通常取一半分析即可。
取一幀信號為例，整個說話段要配合有話幀檢測（VAD技術，這里不提了）自相關法估計基頻code：
[x, fs] = wavread('1.wav');
nw = 256;
signal = x(16001:(16000+nw),1);%取一幀信號
%相關函數計算
r = my_direct_convolution(signal,signal(end:-1:1)); %利用卷積計算相關函數
r = r/(signal'*signal);%相關函數歸一化
rhalf = r(nw:end); % 取延遲量為正值的部分
%提取基音，假設介於50~600Hz之間
lmin = round((50/fs)*nw);
lmax = round((600/fs)*nw);
[tmax,tloc] = max(rhalf(lmin:lmax)); % 在Pmin～Pmax范圍內尋找最大值
pos = lmin+tloc-1; % 給出對應最大值的延遲量
pitch = pos/nw*fs %估計得基頻
卷積還是調用上面的函數。波形的短時周期還是比較明顯的，給出一幀信號：
九、共振峰
首先給出共振峰定義：當聲門處准周期脈沖激勵進入聲道時會引起共振特性，產生一組共振頻率，這一組共振頻率稱為共振峰頻率或簡稱共振峰。
共振峰參數包括共振峰頻率和頻帶的寬度，它是區別不同韻母的重要參數，由於共振峰包含在語音的頻譜包絡中，因此共振峰參數的提取關鍵是估計自然語音的頻譜包絡，並認為譜包括的極大值就是共振峰，通常認為共振峰數量不超過4個。發聲模型：
對這個模型抽象，通常有聲管模型、聲道模型兩個思路，以聲道模型為例：認為信號經過與聲道的卷積，得到最終發出的聲音。聲道就是系統H。
共振峰的求解思路非常多，這里給出一個基於LPC內插的例子。
如果表達這個系統響應H呢？一個基本的思路就是利用全極點求取：
從按照時域里信號與聲道卷積的思路，頻域就是信號與聲道相乘，所以聲道就是包絡：
利用線性預測系數（LPC）並選擇適當階數，可以得到聲道模型H。LPC之前有分析過，可以點擊這里。
其實利用LPC得出的包絡，找到四個峰值，就已經完成了共振峰的估計。為了更精確地利用LPC求取共振峰有，兩種常用思路：拋物線內插法、求根法，這里以內插法為例。其實就是一個插值的思路，如圖
為了估計峰值，取p(k-1)、p(k)、p(k+1)哪一個點都是不合適的，插值構造拋物線，找出峰值就更精確了。
LPC內插法code:
[x,fs]=wavread(fle); % 讀入一幀語音信號
u=filter([1 -.99],1,x); % 預加重
wlen=length(u); % 幀長
p=12; % LPC階數
a=lpc(u,p); % 求出LPC系數
U=lpcar2pf(a,255); % 由LPC系數求出頻譜曲線
freq=(0:256)*fs/512; % 頻率刻度
df=fs/512; % 頻率分辨率
U_log=10*log10(U); % 功率譜分貝值
[Loc,Val]=findpeaks(U); % 在U中尋找峰值
ll=length(Loc); % 有幾個峰值
%拋物線內插修正
for k=1 : ll
m=Loc(k); % 設置m-1,m和m+1
m1=m-1; m2=m+1;
p=Val(k); % 設置P(m-1),P(m)和P(m+1)
p1=U(m1); p2=U(m2);
aa=(p1+p2)/2-p;
bb=(p2-p1)/2;
cc=p;
dm=-bb/2/aa;
pp=-bb*bb/4/aa+cc;
m_new=m+dm;
bf=-sqrt(bb*bb-4*aa*(cc-pp/2))/aa;
F(k)=(m_new-1)*df;
Bw(k)=bf*df;
end
為了更好地估計聲道，這里用了預加重，因為類似電磁波等信號，波信號傳播過程中高頻分量衰減更大，所以需要對高頻進行一定程度的提升，這個操作叫做：預加重。
對應的效果圖：
求解的共振峰頻率（Hz）:676.15 1372.53 2734.15 3513.69，基音周期與共振峰不是一回事啊！有時為了表征聲道特性，也可以直接利用LPC系數作為特征，而不必求取共振峰。
梅爾倒譜系數(MFCC)打算單拎出來寫了，涉及的概念有點多，其他特征用到再補充吧，上一篇文章也提到了很多音頻特征的概念。
參考：
- 宋知用《MATLAB在語音信號分析與合成中的應用》