MATLAB 卡尔曼滤波

PIZI

   这一段时间对现代滤波进行了学习，对自适应滤波器和卡尔曼滤波器有了一定认识，并对它们用MATLAB对语音信号进行了滤波，发现卡尔曼滤波器还是比较有用，能够在较大的噪声中还原原来的信号。新的学期马上就开始了，由于TI的开发板一直在维修，所以学习TI开发板的计划搁置，但是对声音信号的处理及滤波器的认识有了进一步提高。新的学期继续努力！  
    卡尔曼滤波的基本思想是：以最小均方误差为最佳估计准则，采用信号与噪声的状态空间模型，利用前一时刻的估计值和当前时刻的观测值来更新对状态变量的估计，求出当前时刻的估计值，算法根据建立的系统方程和观测方程对需要处理的信号做出满足最小均方误差的估计。
语音信号在较长时间内是非平稳的，但在较短的时间内的一阶统计量和二阶统计量近似为常量，因此语音信号在相对较短的时间内可以看成白噪声激励以线性时不变系统得到的稳态输出。假定语音信号可看成由一AR模型产生：
                     
     时间更新方程：
                        
    测量更新方程：
                       
    K(t)为卡尔曼增益，其计算公式为：
                        
其中
                                    
、 分别为过程模型噪声协方差和测量模型噪声协方差，测量协方差可以通过观测得到，则较难确定，在本实验中则通过与两者比较得到。
     由于语音信号短时平稳，因此在进行卡尔曼滤波之前对信号进行分帧加窗操作，在滤波之后对处理得到的信号进行合帧，这里选取帧长为256，而帧重叠个数为128；
     下图为原声音信号与加噪声后的信号以及声音信号与经卡尔曼滤波处理后的信号：


                                原声音信号与加噪声后的信号

                               原声音信号与经卡尔曼滤波处理后的信号
MATLAB程序实现如下：
%%%%%%%%%%%%%%%%%基于LPC全极点模型的最大后验概率估计法,采用卡尔曼滤波%%%%%%%%%%%%%%
clear;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%加载声音数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load voice.mat
y=m1(2,:);
x=y+0.08*randn(1,length(y));
%%%%%%%%%%%%%%%原声音信号和加噪声后的信号%%%%%%%%%%%%%%%
figure(1);
subplot(211);plot(m1(1,:),m1(2,:));xlabel('时间');ylabel('幅度');title('原声音信号');
subplot(212);plot(m1(1,:),x);xlabel('时间');ylabel('幅度');title('加噪声后的信号');
%%%%%%%%%%%%%%%%%%%%%%%%%输入参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fs=44100;                      %信号采样的频率
bits=16;                 %信号采样的位数
N=256;                      %帧长
m=N/2;                      %每帧移动的距离
lenth=length(x);            %输入信号的长度
count=floor(lenth/m)-1;     %处理整个信号需要移动的帧数%%%先不考虑补零的问题
p=11;                             %AR模型的阶数
a=zeros(1,p);
w=hamming(N);                   %加汉明窗函数
y_temp=0;
F=zeros(11,11);           %转移矩阵
F(1,2)=1;
F(2,3)=1;
F(3,4)=1;
F(4,5)=1;
F(5,6)=1;
F(6,7)=1;
F(7,8)=1;
F(8,9)=1;
F(9,10)=1;
F(10,11)=1;
H=zeros(1,p);                        %
S0=zeros(p,1);
P0=zeros(p);
S=zeros(p);
H(11)=1;
s=zeros(N,1);
G=H';
P=zeros(p);
%%%%%%%%%%%%%%%%测试噪声协方差%%%%%%%%%%%%%%%%%%%%%%
y_temp=cov(x(1:7680));
x_frame=zeros(256,1);
x_frame1=zeros(256,1);
T=zeros(lenth,1);
for r=1:count
%%%%%%%%%%%%%%%%%%%5%%%%%分帧处理%%%%%%%%%%%%%%%%%%%%%
        x_frame=x((r-1)*m+1:(r+1)*m);
%%%%%%%%%%%%%%%%采用LPC模型求转移矩阵参数%%%%%%%%%%%%%%   
          if r==1
           [a,VS]=lpc(x_frame(:),p);  
         else
           [a,VS]=lpc(T((r-2)*m+1:(r-2)*m+256),p);
        end
%%%%%%%%%%%%%%%%帧长内过程噪声协方差%%%%%%%%%%%%%%%%%%
        if (VS-y_temp>0)   
            VS=VS-y_temp;
        else
            VS=0.0005;
        end
  
        F(p,:)=-1*a(p+1:-1:2);

        for j=1:256
            if(j==1)
            S=F*S0;
            Pn=F*P*F'+G*VS*G';
            else
            S=F*S;      %时间更新方程
            Pn=F*P*F'+G*VS*G';
            end
            K=Pn*H'*(y_temp+H*P*H').^(-1); %卡尔曼增益
            P=(eye(p)-K*H)*Pn;                     %测量更新方程
            S=S+K*[x_frame(j)-H*S];
            T((r-1)*m+j)=H*S;
        end
%%%%%%%%%%%%%%%%对得到的每帧数据进行加窗操作%%%%%%%%%%%%%%%%%%%%%%%%
        ss(1:256,r)=T((r-1)*m+1:(r-1)*m+256);
         sss(1:256,r)=ss(1:256,r).*w;
     
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%合帧操作%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         for r=1:count
             if r==1
            s_out(1:128)=sss(1:128,r);
            else if r==count
            s_out(r*m+1:r*m+m)=sss(129:256,r);
            else
             s_out(((r-1)*m+1):((r-1)*m+m))=sss(129:256,r-1)+sss(1:128,r);
             end
            end
       end

figure(2)
subplot(211);plot(m1(1,:),m1(2,:));xlabel('时间');ylabel('幅度');title('原声音信号');
subplot(212);plot((1:1109760)/Fs,s_out);xlabel('时间');ylabel('幅度');title('经卡尔曼滤波后的声音信号');

tanyueye

楼主能把voice.mat文件贴上来吗。想练练手~~~~

造字仓颉

熟悉一下程序的编写。怎么同化呢？

EMD

求共享经验啊

panshi

学习一下{:5_213:}{:5_213:}

Howie721

谢谢楼主分享，正在学习卡尔曼滤波~

Howie721

谢谢楼主分享，正在学习卡尔曼滤波~

一路前行

请问一下，如何把ar模型与卡尔曼，总觉得没训练就直接给了系统的状态方程

无锡气象台

要是有Fortran版的就好了

沐..

为什么这个滤波程序跟其他帖子的差异有点大，其他的好多个代码文件我都看不懂