卡爾曼濾波算法及MATLAB實現(xiàn)

1、基于matlab的卡爾曼信號濾波設計卡爾曼濾波的基本思想是：以最小均方誤差為最佳估計準則，采用信號與噪聲的狀態(tài)空間 模型，利用前一時刻的估計值和當前時刻的觀測值來更新對狀態(tài)變量的估計，求出當前時刻的估計值，算法根據(jù)建立的系統(tǒng)方程和觀測方程對需要處理的信號做出滿足最小均方誤差的 估計。語音信號在較長時間內(nèi)是非平穩(wěn)的，但在較短的時間內(nèi)的一階統(tǒng)計量和二階統(tǒng)計量近似為常量，因此語音信號在相對較短的時間內(nèi)可以看成白噪聲激勵以線性時不變系統(tǒng)得到的穩(wěn)態(tài)輸 出。假定語音信號可看成由一AR模型產(chǎn)生：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、 分別為過程模型噪聲協(xié)方差和測量模型噪聲協(xié)方差，測量協(xié)方差可以通過觀測得到,則較難確定，在本實驗中則通過與兩者比較得到。由于語音信號短時平穩(wěn)，因此在進行卡爾曼濾波之前對信號進行分幀加窗操作，在濾波之后對處理得到的信號進行合幀，這里選取幀長為256,而幀重疊個數(shù)為128;卜圖為原聲音信號與加噪聲后的信號以及聲音信號與經(jīng)卡爾曼濾波處理后的信號原聲音信號與加噪聲后的信號MATLAB程序?qū)崿F(xiàn)如下：%基于LPC全極點模型的最大后驗概率估計法,采用卡爾曼濾 波%clear;cl

3、c;% 力口 載 聲 音 9據(jù)%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(加噪聲后的信號);%輸入?yún)?shù) %Fs=44100;%信號采樣的頻率bits=16;%信號采樣的位數(shù)N=256;% 幀長m=N/2;%每幀移動的距離lenth=leng

4、th(x); %輸入信號的長度%AR模型的階數(shù)%加漢明窗函數(shù)%轉(zhuǎn)移矩陣count=floor(lenth/m)-1; %處理整個信號需要移動的幀數(shù)%先不考慮補零的問題 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);% 測試噪聲協(xié)方差 % 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模型求轉(zhuǎn)移矩陣參數(shù) %if r=1a,VS=lpc(x_frame(:),p);elsea,VS=lpc(T(r-2)*m+1:(r-2)*m+256),p); end% 幀長內(nèi)過程噪聲協(xié)方差 %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% 對得到的每幀數(shù)據(jù)進行力口窗操 作%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