基于扩展卡尔曼滤波EKF的语音信号基音估计算法matlab仿真

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