1.算法仿真效果
matlab2022a仿真结果如下：
2.算法涉及理论知识概要
基音是语音信号的基本频率成分，它决定了语音的音调和声音的音高。在语音信号处理中，基音估计是一个重要的任务，它可以用于语音合成、语音识别、语音增强等应用。扩展卡尔曼滤波（Extended Kalman Filter, EKF）是一种用于非线性系统的滤波方法，它可以用于基音的估计。
在语音信号中，周期性的振动成分被称为基音。基音的周期是指相邻两个周期波形的时间间隔，也称为基音周期。频率是指每秒钟振动的周期数，它的倒数称为周期。对于一个周期为T的基音，其频率f = 1/T。基音的频率范围通常在50Hz-500Hz之间。
卡尔曼滤波（Kalman Filter, KF）是一种用于线性系统的滤波方法，它可以在有噪声的观测数据中，根据已知的系统模型和初始状态，推断出系统的状态。扩展卡尔曼滤波是一种用于非线性系统的滤波方法，它通过在每个时间步骤使用局部线性化来近似非线性系统，并使用卡尔曼滤波来进行状态估计。
扩展卡尔曼滤波需要一个系统模型，它描述了基音的演化规律。在基音估计中，系统模型可以表示为：
x(k) = A(k-1)x(k-1) + w(k-1)
其中，x(k)表示在时间k时的状态向量，A(k-1)表示状态转移矩阵，w(k-1)表示系统噪声。在基音估计中，状态向量可以表示为：
x(k) = [p(k), T(k)]
其中，p(k)表示基音周期，T(k)表示基音的相位。状态转移矩阵A(k-1)可以表示为：
A(k-1) = [1 0; 0 1]
这个矩阵表示基音周期和相位在每个时间步骤中保持不变。系统噪声w(k-1)可以表示为：
w(k-1) = [w1(k-1), w2(k-1)]
其中，w1(k-1)和w2(k-1)分别表示基音周期和相位的噪声。
扩展卡尔曼滤波还需要一个观测模型，它描述了观测数据和状态向量之间的关系。在基音估计中，观测模型可以表示为：
y(k) = H(k)x(k) + v(k)
其中，y(k)表示在时间k时的观测向量，H(k)表示观测矩阵，v(k)表示观测噪声。在基音估计中，观测向量可以表示为：
y(k) = [y1(k), y2(k)]
其中，y1(k)和y2(k)分别表示基音周期和相位的观测值。观测矩阵H(k)可以表示为：
H(k) = [1 0; 0 1]
这个矩阵表示我们可以直接观测到基音周期和相位。观测噪声v(k)可以表示为：
v(k) = [v1(k), v2(k)]
其中，v1(k)和v2(k)分别表示基音周期和相位的噪声。
扩展卡尔曼滤波算法可以分为两个步骤：预测和更新。在预测步骤中，我们使用系统模型来预测下一个时间步骤的状态向量和协方差矩阵。在更新步骤中，我们使用观测模型来根据观测数据来更新预测值。下面是扩展卡尔曼滤波算法的详细步骤：
初始化状态向量和协方差矩阵：
x(0) = [p(0), T(0)]
P(0) = diag([p_var(0), T_var(0)])
对于每个时间步骤k：
a. 预测步骤：
根据系统模型，预测下一个时间步骤的状态向量：
x(k|k-1) = A(k-1)x(k-1|k-1)
根据系统模型，预测下一个时间步骤的协方差矩阵：
P(k|k-1) = A(k-1)P(k-1|k-1)A(k-1)^T + Q(k-1)
b. 更新步骤：
计算卡尔曼增益K(k)：
K(k) = P(k|k-1)H(k)^T(H(k)P(k|k-1)H(k)^T + R(k))^(-1)
根据观测数据，计算当前时间步骤的状态向量：
x(k|k) = x(k|k-1) + K(k)(y(k) - H(k)x(k|k-1))
根据观测数据，计算当前时间步骤的协方差矩阵：
P(k|k) = (I - K(k)H(k))P(k|k-1)
其中，Q(k-1)表示系统噪声的协方差矩阵，R(k)表示观测噪声的协方差矩阵。对于基音估计，我们可以将Q(k-1)和R(k)设置为常数，如下所示：
Q(k-1) = diag([q1, q2])
R(k) = diag([r1, r2])
其中，q1和q2分别表示基音周期和相位的噪声方差，r1和r2分别表示基音周期和相位的观测噪声方差。
3.MATLAB核心程序
%pitch tracking
for ii=2:size(datass,2)
%基于先前估计的均值一步预测
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);
% 计算一步预测残差
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
%计算除最后一个步骤外的所有步骤的预测步骤
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
