1.算法仿真效果

matlab2022a仿真結(jié)果如下：

2.算法涉及理論知識(shí)概要

基音是語(yǔ)音信號(hào)的基本頻率成分，它決定了語(yǔ)音的音調(diào)和聲音的音高。在語(yǔ)音信號(hào)處理中，基音估計(jì)是一個(gè)重要的任務(wù)，它可以用于語(yǔ)音合成、語(yǔ)音識(shí)別、語(yǔ)音增強(qiáng)等應(yīng)用。擴(kuò)展卡爾曼濾波（Extended Kalman Filter, EKF）是一種用于非線性系統(tǒng)的濾波方法，它可以用于基音的估計(jì)。

在語(yǔ)音信號(hào)中，周期性的振動(dòng)成分被稱為基音?；舻闹芷谑侵赶噜弮蓚€(gè)周期波形的時(shí)間間隔，也稱為基音周期。頻率是指每秒鐘振動(dòng)的周期數(shù)，它的倒數(shù)稱為周期。對(duì)于一個(gè)周期為T的基音，其頻率f = 1/T。基音的頻率范圍通常在50Hz-500Hz之間。

卡爾曼濾波（Kalman Filter, KF）是一種用于線性系統(tǒng)的濾波方法，它可以在有噪聲的觀測(cè)數(shù)據(jù)中，根據(jù)已知的系統(tǒng)模型和初始狀態(tài)，推斷出系統(tǒng)的狀態(tài)。擴(kuò)展卡爾曼濾波是一種用于非線性系統(tǒng)的濾波方法，它通過(guò)在每個(gè)時(shí)間步驟使用局部線性化來(lái)近似非線性系統(tǒng)，并使用卡爾曼濾波來(lái)進(jìn)行狀態(tài)估計(jì)。

擴(kuò)展卡爾曼濾波需要一個(gè)系統(tǒng)模型，它描述了基音的演化規(guī)律。在基音估計(jì)中，系統(tǒng)模型可以表示為：

x(k) = A(k-1)x(k-1) + w(k-1)

其中，x(k)表示在時(shí)間k時(shí)的狀態(tài)向量，A(k-1)表示狀態(tài)轉(zhuǎn)移矩陣，w(k-1)表示系統(tǒng)噪聲。在基音估計(jì)中，狀態(tài)向量可以表示為：

x(k) = [p(k), T(k)]

其中，p(k)表示基音周期，T(k)表示基音的相位。狀態(tài)轉(zhuǎn)移矩陣A(k-1)可以表示為：

A(k-1) = [1 0; 0 1]

這個(gè)矩陣表示基音周期和相位在每個(gè)時(shí)間步驟中保持不變。系統(tǒng)噪聲w(k-1)可以表示為：

w(k-1) = [w1(k-1), w2(k-1)]

其中，w1(k-1)和w2(k-1)分別表示基音周期和相位的噪聲。

擴(kuò)展卡爾曼濾波還需要一個(gè)觀測(cè)模型，它描述了觀測(cè)數(shù)據(jù)和狀態(tài)向量之間的關(guān)系。在基音估計(jì)中，觀測(cè)模型可以表示為：

y(k) = H(k)x(k) + v(k)

其中，y(k)表示在時(shí)間k時(shí)的觀測(cè)向量，H(k)表示觀測(cè)矩陣，v(k)表示觀測(cè)噪聲。在基音估計(jì)中，觀測(cè)向量可以表示為：

y(k) = [y1(k), y2(k)]

其中，y1(k)和y2(k)分別表示基音周期和相位的觀測(cè)值。觀測(cè)矩陣H(k)可以表示為：

H(k) = [1 0; 0 1]

這個(gè)矩陣表示我們可以直接觀測(cè)到基音周期和相位。觀測(cè)噪聲v(k)可以表示為：

v(k) = [v1(k), v2(k)]

其中，v1(k)和v2(k)分別表示基音周期和相位的噪聲。

擴(kuò)展卡爾曼濾波算法可以分為兩個(gè)步驟：預(yù)測(cè)和更新。在預(yù)測(cè)步驟中，我們使用系統(tǒng)模型來(lái)預(yù)測(cè)下一個(gè)時(shí)間步驟的狀態(tài)向量和協(xié)方差矩陣。在更新步驟中，我們使用觀測(cè)模型來(lái)根據(jù)觀測(cè)數(shù)據(jù)來(lái)更新預(yù)測(cè)值。下面是擴(kuò)展卡爾曼濾波算法的詳細(xì)步驟：

初始化狀態(tài)向量和協(xié)方差矩陣：

x(0) = [p(0), T(0)]

P(0) = diag([p_var(0), T_var(0)])

對(duì)于每個(gè)時(shí)間步驟k：

a. 預(yù)測(cè)步驟：

根據(jù)系統(tǒng)模型，預(yù)測(cè)下一個(gè)時(shí)間步驟的狀態(tài)向量：

x(k|k-1) = A(k-1)x(k-1|k-1)

根據(jù)系統(tǒng)模型，預(yù)測(cè)下一個(gè)時(shí)間步驟的協(xié)方差矩陣：

P(k|k-1) = A(k-1)P(k-1|k-1)A(k-1)^T + Q(k-1)

b. 更新步驟：

計(jì)算卡爾曼增益K(k)：

K(k) = P(k|k-1)H(k)^T(H(k)P(k|k-1)H(k)^T + R(k))^(-1)

根據(jù)觀測(cè)數(shù)據(jù)，計(jì)算當(dāng)前時(shí)間步驟的狀態(tài)向量：

x(k|k) = x(k|k-1) + K(k)(y(k) - H(k)x(k|k-1))

根據(jù)觀測(cè)數(shù)據(jù)，計(jì)算當(dāng)前時(shí)間步驟的協(xié)方差矩陣：

P(k|k) = (I - K(k)H(k))P(k|k-1)

其中，Q(k-1)表示系統(tǒng)噪聲的協(xié)方差矩陣，R(k)表示觀測(cè)噪聲的協(xié)方差矩陣。對(duì)于基音估計(jì)，我們可以將Q(k-1)和R(k)設(shè)置為常數(shù)，如下所示：

Q(k-1) = diag([q1, q2])

R(k) = diag([r1, r2])

其中，q1和q2分別表示基音周期和相位的噪聲方差，r1和r2分別表示基音周期和相位的觀測(cè)噪聲方差。

3.MATLAB核心程序

%pitch tracking

for ii=2:size(datass,2)

%基于先前估計(jì)的均值一步預(yù)測(cè)

One_step_state=F*(state(:,ii-1));

P_OneStep(:,:,ii)=F*P(:,:,ii-1)*F'+C*Q*C';

H=cos((B*One_step_state)'+pha')*G-(G*One_step_state)'*diag(sin(B*One_step_state+pha))*(B);

O_covariance=(H*P_OneStep(:,:,ii)*H'+R);

% Kalman gain

K=P_OneStep(:,:,ii)*H'*O_covariance^(-1);

% 計(jì)算一步預(yù)測(cè)殘差

h=(G*One_step_state)'*cos(B*One_step_state+pha);

correction_factor=K*(datass(:,ii)-h);

state(:,ii)= One_step_state+correction_factor;

P(:,:,ii) ?= P_OneStep(:,:,ii)-K*H*P_OneStep(:,:,ii); ?

end

%卡爾曼平滑器;

N=size(datass,2);

pitch(:,N) ????= state(:,N);

P_upS(:,:,N) ??= P(:,:,N);

for k = (N-1):-1:1

%計(jì)算除最后一個(gè)步驟外的所有步驟的預(yù)測(cè)步驟

sgain = (P(:,:,k)*F')/(F*P(:,:,k)*F' + C*Q*C');

pitch(:,k) ????= state(:,k) ?+ sgain*(pitch(:,k+1) ?- F*(state(:,k)));

P_upS(:,:,k) ??= P(:,:,k)+ sgain*(P_upS(:,:,k+1) - P_OneStep(:,:,k+1))*sgain';

end

end