语音增强--卡尔曼滤波介绍及MATLAB实现

语音增强--------------卡尔曼滤波

状态方程

$${\mathbf{x}}_{k+1}={\mathbf{\Phi }}_k{\mathbf{x}}_k+\mathbf{\Gamma }{\mathbf{u}}_k$$

观测方程

$${\mathbf{y}}_{k+1}={\mathbf{H}}_{k+1}{\mathbf{x}}_{k+1}+{\mathbf{n}}_{k+1}$$

原理介绍

假设在时刻$t_k$，基于$t_k$时刻以前所获得的全部知识，对状态变量${\mathbf{x}}_k$做出一个预测估计，记为${\hat{\mathbf{x}}}_k^{-}$，则预测估计的误差为
$${\mathbf{e}}_k^{-}={\mathbf{x}}_k-{\hat{\mathbf{x}}}_k^{-}$$
称为预测误差。预测误差是零均值的，其协方差矩阵为
$${\mathbf{C}}_k^{-}=E\left\{{\mathbf{e}}_k^{-}\mathbf{e}{{}_k^{-}}^T\right\}=E\{\left({\mathbf{x}}_k-{\hat{\mathbf{x}}}_k^{-}\right){\left({\mathbf{x}}_k-{\hat{\mathbf{x}}}_k^{-}\right)}^T\}$$
在预测估计${\hat{\mathbf{x}}}_k^{-}$的基础上，利用$t_k$时刻所获取的新观测数据${\mathbf{y}}_k$来进一步改善对${\mathbf{x}}_k$的估计，记为 ${\hat{\mathbf{x}}}_k$，称为更新估计，更新估计通过下式完成：
$${\hat{\mathbf{x}}}_k={\hat{\mathbf{x}}}_k^{-}+{\mathbf{K}}_k\left({\mathbf{y}}_k-{\mathbf{H}}_k{\hat{\mathbf{x}}}_k^{-}\right)$$
其中 ${\mathbf{K}}_k$为待定的增益矩阵，称为卡尔曼增益。
更新估计的误差记作${\mathbf{e}}_k$，则有
$${\mathbf{e}}_k={\mathbf{x}}_k-{\hat{\mathbf{x}}}_k={\mathbf{x}}_k-[{\hat{\mathbf{x}}}_k^{-}+{\mathbf{K}}_k\left({\mathbf{y}}_k-{\mathbf{H}}_k{\hat{\mathbf{x}}}_k^{-}\right)]$$
矢量卡尔曼滤波的实质就是寻找适当的增益矩阵${\mathbf{K}}_k$，使更新估计的均方误差达到最小。
具体推导过程可以参照参考文献[1]，这里只给出最终的递推过程
(1) 建立状态空间模型
$$\begin{array}{rl}& {\mathbf{x}}_{k+1}={\mathbf{\Phi }}_k{\mathbf{x}}_k+\mathbf{\Gamma }{\mathbf{u}}_k\\ & {\mathbf{y}}_{k+1}={\mathbf{H}}_{k+1}{\mathbf{x}}_{k+1}+{\mathbf{n}}_{k+1}\end{array}$$
(2) 设置初始化条件
x ^ 0 = E { x 0 } C 0 = V a r { x 0 } \begin{aligned} & {{{\mathbf{\hat{x}}}}_{0}}=E\left\{ {{\mathbf{x}}_{0}} \right\} \\ & {{\mathbf{C}}_{0}}=Var\left\{ {{\mathbf{x}}_{0}} \right\} \end{aligned} ​x^0​=E{x0​}C0​=Var{x0​}​
(3) 预测
$${\hat{\mathbf{x}}}_{k+1}^{-}={\mathbf{\Phi }}_k{\hat{\mathbf{x}}}_k$$
(4) 计算预测误差的协方差
$${\mathbf{C}}_{k+1}^{-}={\mathbf{\Phi }}_k{\mathbf{C}}_k{\mathbf{\Phi }}_k^T+\mathbf{\Gamma }{\mathbf{Q}}_k{\mathbf{\Gamma }}^T$$
(5) 计算卡尔曼增益
$${\mathbf{K}}_{k+1}={\mathbf{C}}_{k+1}^{-}{\mathbf{H}}_{k+1}^T{\left({\mathbf{H}}_{k+1}{\mathbf{C}}_{k+1}^{-}{\mathbf{H}}_{k+1}^T+{\mathbf{R}}_{k+1}\right)}^{-1}$$
(6) 更新
$${\hat{\mathbf{x}}}_{k+1}={\hat{\mathbf{x}}}_{k+1}^{-}+{\mathbf{K}}_{k+1}\left({\mathbf{y}}_{k+1}-{\mathbf{H}}_{k+1}{\hat{\mathbf{x}}}_{k+1}^{-}\right)$$
(7) 估计误差的协方差
$${\mathbf{C}}_{k+1}=\left(\mathbf{I}-{\mathbf{K}}_{k+1}{\mathbf{H}}_{k+1}\right){\mathbf{C}}_{k+1}^{-}$$
(8) 令$k=k+1$ ，重复步骤(3)-(8)直到当前时刻
其中
$$E\left\{{\mathbf{n}}_k{\mathbf{n}}_i^T\right\}=\{\begin{array}{rl}& {\mathbf{R}}_{k}i=k\\ & 0i\ne k\end{array}$$
$$E\left\{{\mathbf{u}}_k{\mathbf{u}}_i^T\right\}=\{\begin{array}{rl}& {\mathbf{Q}}_{k}i=k\\ & 0i\ne k\end{array}$$
仿真参数为

参数名称	参数值
信噪比	5dB
采样率	16KHz

仿真结果如下

随着卡尔曼滤波收敛的过程中，误差的协方差矩阵会变得越来越小，从图中也可以看出，经过卡尔曼滤波之后的语音信号，噪声很大程度地被滤除，进一步地提高了语音的质量。

主函数

clear;
close all;
clc;
%% 读入数据
[signal,~]=audioread('clean.wav');      %读入干净语音
[noise,fs]=audioread('noise.wav');       %读入噪声
N=3*fs;                                              %选取3秒的语音
signal=signal(1:N);
noise=noise(1:N);
t=(0:N-1)/fs;
SNR=5;                 %信噪比大小
noise=noise/norm(noise,2).*10^(-SNR/20)*norm(signal);
x=signal+noise;          %产生固定信噪比的带噪语音
Time = (0:1/fs:(length(signal)-1)/fs)';      %时间轴
Noise=x(1:fs,1);                     %选取前1秒语音作为噪声方差估计
len_win = 0.0025;       % 窗长2.5ms
shift_percent = 1;       % 窗移占比
AR_order = 20;           % 滤波器阶数
iter = 7;                      %迭代次数设置
%% 分帧加窗处理
len_winframe = fix(len_win * fs);
window = ones(len_winframe,1);
[y, num_frame] = KFrame(x, len_winframe, window, shift_percent);

%% 初始化
H = [zeros(1,AR_order-1),1];   % 观测矩阵
R = var(Noise);                       % 噪声方差
[filt_coeff, Q] = lpc(y, AR_order);              % LPC预测，得到滤波器的系数
C = R * eye(AR_order,AR_order);              % 误差协方差矩阵
enhanced_speech = zeros(1,length(x));    % 增强后的语音信号
enhanced_speech(1:AR_order) = x(1:AR_order,1)';   %初始化
updata_x = x(1:AR_order,1);

% 迭代器的次数.
i = AR_order+1;
j = AR_order+1;

%% 卡尔曼滤波
for k = 1:num_frame   %一次处理一帧信号
    jStart = j;     %跟踪每次迭代AR_Order+1的值.
    OutputOld = updata_x;    %为每次迭代保留第一批AROrder预估量
    
    for l = 1:iter               %迭代次数
        fai = [zeros(AR_order-1,1) eye(AR_order-1); fliplr(-filt_coeff(k,2:end))];
        
        for ii = i:len_winframe
            %% 卡尔曼滤波
            predict_x = fai * updata_x;
            predict_C = (fai * C * fai') + (H' * Q(k) * H);
            K = (predict_C * H')/((H * predict_C * H') + R);
            updata_x = predict_x + (K * (y(ii,k) - (H*predict_x)));
            enhanced_speech(j-AR_order+1:j) = updata_x';
            C = (eye(AR_order) - K * H) * predict_C;
            j = j+1;
        end
        i = 1;
        if l < iter
            j = jStart;
            updata_x = OutputOld;
        end
        % 更新滤波后信号的lpc
        [filt_coeff(k,:), Q(k)] = lpc(enhanced_speech((k-1)*len_winframe+1:k*len_winframe),AR_order);
    end
end
enhanced_speech = enhanced_speech(1:N)';
figure(1)
subplot(321);
plot(t,signal);ylim([-1.5,1.5]);title('干净语音');xlabel('时间/s');ylabel('幅度');
subplot(323);
plot(t,x);ylim([-1.5,1.5]);title('带噪语音');xlabel('时间/s');ylabel('幅度');
subplot(325);
plot(t,real(enhanced_speech));ylim([-1.5,1.5]);title('卡尔曼滤波增强后的语音');xlabel('时间/s');ylabel('幅度');
subplot(322);
spectrogram(signal,256,128,256,16000,'yaxis');
subplot(324);
spectrogram(x,256,128,256,16000,'yaxis');
subplot(326);
spectrogram(enhanced_speech,256,128,256,16000,'yaxis');

KFrame.m

function [Output, NumSegments] = KFrame(Input, WindowLength, Window, HoppingSize)
% Chopper windows the signal based on window length, shift percantage and
% uses Hamming windowing technique.

% Number of samples to hop.
HoppingSamples = fix(WindowLength.*HoppingSize);

% Number of segments.
NumSegments = fix(((length(Input)-WindowLength)/HoppingSamples) + 1);

% Index matrix which guides the signal through chopping process.
Index = (repmat(1:WindowLength,NumSegments,1) + repmat((0:(NumSegments-1))'*HoppingSamples,1,WindowLength))';

% Final window which multiplies with original signal to give pieces of it.
FinalWindow = repmat(Window,1,NumSegments);

% Ta-da... 
Output = Input(Index).*FinalWindow;
end

关于语音及噪声文件，具体请参考：语音信号处理常用语料库下载地址

参考文献：
[1] 叶中付. 统计信号处理[M]. 中国科学技术大学出版社, 2013.