卡尔曼滤波算法及MATLAB实现

1、基于matlab的卡尔曼信号滤波设计卡尔曼滤波的基本思想是：以最小均方误差为最佳估计准则，采用信号与噪声的状态空间 模型，利用前一时刻的估计值和当前时刻的观测值来更新对状态变量的估计，求出当前时刻的估计值，算法根据建立的系统方程和观测方程对需要处理的信号做出满足最小均方误差的 估计。语音信号在较长时间内是非平稳的，但在较短的时间内的一阶统计量和二阶统计量近似为常量，因此语音信号在相对较短的时间内可以看成白噪声激励以线性时不变系统得到的稳态输 出。假定语音信号可看成由一AR模型产生：W 户=工 akS(t-k)+u(t时间更新方程：却出力超maP|t|M)=FP(t-l|t-lFTiggT测量更

2、新方程：出t|上 都t-l)K为卡尔曼增益，其计算公式为:出止其中Xfn)=x(n-p+lj xn-p+2)-上GOO i0 W、 分别为过程模型噪声协方差和测量模型噪声协方差，测量协方差可以通过观测得到,则较难确定，在本实验中则通过与两者比较得到。由于语音信号短时平稳，因此在进行卡尔曼滤波之前对信号进行分帧加窗操作，在滤波之后对处理得到的信号进行合帧，这里选取帧长为256,而帧重叠个数为128;卜图为原声音信号与加噪声后的信号以及声音信号与经卡尔曼滤波处理后的信号原声音信号与加噪声后的信号MATLAB程序实现如下：%基于LPC全极点模型的最大后验概率估计法,采用卡尔曼滤 波%clear;cl

3、c;% 力口 载 声 音 9据%loadvoice.maty=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=leng

4、th(x); %输入信号的长度%AR模型的阶数%加汉明窗函数%转移矩阵count=floor(lenth/m)-1; %处理整个信号需要移动的帧数%先不考虑补零的问题 p=11;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,

5、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=1a,VS=lpc(x_frame(:),p);elsea,VS=lpc(T(r-2)*m+1:(r-2)*m+256),p); end% 帧长内过程噪声协方差 %if (VS-y_temp0)VS=VS-y_temp;elseVS=0

6、.0005; endF(p,:)=-1*a(p+1:-1:2);for j=1:256if(j=1)S=F*S0;Pn=F*P*F+G*VS*G;elseS=F*S; %时间更新方程Pn=F*P*F+G*VS*G;endK=Pn*H*(y_temp+H*P*H).A(-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:countif r=1s_out(1:128)=sss(1:128,r);else if r=counts_out(r*m+1:r*m+m)=sss(129:256,r);elses_out(r-1)*m+1):(r-1)*m+m)=sss(129:256,r-1)+sss(1:128,r);endend endfigure(2)subplot(211);plot(m1(1,:),m1(2,:);xlabel(时间);yl