這一段時(shí)間對(duì)現(xiàn)代濾波進(jìn)行了學(xué)習(xí),對(duì)自適應(yīng)濾波器和卡爾曼濾波器有了一定認(rèn)識(shí),并對(duì)它們用MATLAB對(duì)語(yǔ)音信號(hào)進(jìn)行了濾波,發(fā)現(xiàn)卡爾曼濾波器還是比較有用,能夠在較大的噪聲中還原原來(lái)的信號(hào)。新的學(xué)期馬上就開(kāi)始了,由于TI的開(kāi)發(fā)板一直在維修,所以學(xué)習(xí)TI開(kāi)發(fā)板的計(jì)劃擱置,但是對(duì)聲音信號(hào)的處理及濾波器的認(rèn)識(shí)有了進(jìn)一步提高。新的學(xué)期繼續(xù)努力!
卡爾曼濾波的基本思想是:以最小均方誤差為最佳估計(jì)準(zhǔn)則,采用信號(hào)與噪聲的狀態(tài)空間模型,利用前一時(shí)刻的估計(jì)值和當(dāng)前時(shí)刻的觀測(cè)值來(lái)更新對(duì)狀態(tài)變量的估計(jì),求出當(dāng)前時(shí)刻的估計(jì)值,算法根據(jù)建立的系統(tǒng)方程和觀測(cè)方程對(duì)需要處理的信號(hào)做出滿足最小均方誤差的估計(jì)。
語(yǔ)音信號(hào)在較長(zhǎng)時(shí)間內(nèi)是非平穩(wěn)的,但在較短的時(shí)間內(nèi)的一階統(tǒng)計(jì)量和二階統(tǒng)計(jì)量近似為常量,因此語(yǔ)音信號(hào)在相對(duì)較短的時(shí)間內(nèi)可以看成白噪聲激勵(lì)以線性時(shí)不變系統(tǒng)得到的穩(wěn)態(tài)輸出。假定語(yǔ)音信號(hào)可看成由一AR模型產(chǎn)生:
時(shí)間更新方程:
測(cè)量更新方程:
K(t)為卡爾曼增益,其計(jì)算公式為:
其中
、分別為過(guò)程模型噪聲協(xié)方差和測(cè)量模型噪聲協(xié)方差,測(cè)量協(xié)方差可以通過(guò)觀測(cè)得到,則較難確定,在本實(shí)驗(yàn)中則通過(guò)與兩者比較得到。
由于語(yǔ)音信號(hào)短時(shí)平穩(wěn),因此在進(jìn)行卡爾曼濾波之前對(duì)信號(hào)進(jìn)行分幀加窗操作,在濾波之后對(duì)處理得到的信號(hào)進(jìn)行合幀,這里選取幀長(zhǎng)為256,而幀重疊個(gè)數(shù)為128;
下圖為原聲音信號(hào)與加噪聲后的信號(hào)以及聲音信號(hào)與經(jīng)卡爾曼濾波處理后的信號(hào):

原聲音信號(hào)與加噪聲后的信號(hào)
原聲音信號(hào)與經(jīng)卡爾曼濾波處理后的信號(hào)
MATLAB程序?qū)崿F(xiàn)如下:
%%%%%%%%%%%%%%%%%基于LPC全極點(diǎn)模型的最大后驗(yàn)概率估計(jì)法,采用卡爾曼濾波%%%%%%%%%%%%%%
clear;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%加載聲音數(shù)據(jù)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load voice.mat
y=m1(2,:);
x=y+0.08*randn(1,length(y));
%%%%%%%%%%%%%%%原聲音信號(hào)和加噪聲后的信號(hào)%%%%%%%%%%%%%%%
figure(1);
subplot(211);plot(m1(1,:),m1(2,:));xlabel('時(shí)間');ylabel('幅度');title('原聲音信號(hào)');
subplot(212);plot(m1(1,:),x);xlabel('時(shí)間');ylabel('幅度');title('加噪聲后的信號(hào)');
%%%%%%%%%%%%%%%%%%%%%%%%%輸入?yún)?shù)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%
Fs=44100; %信號(hào)采樣的頻率
bits=16; %信號(hào)采樣的位數(shù)
N=256;%幀長(zhǎng)
m=N/2;%每幀移動(dòng)的距離
lenth=length(x);%輸入信號(hào)的長(zhǎng)度
count=floor(lenth/m)-1;%處理整個(gè)信號(hào)需要移動(dòng)的幀數(shù)%%%先不考慮補(bǔ)零的問(wèn)題
p=11;%AR模型的階數(shù)
a=zeros(1,p);
w=hamming(N);%加漢明窗函數(shù)
y_temp=0;
F=zeros(11,11);%轉(zhuǎn)移矩陣
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);
%%%%%%%%%%%%%%%%測(cè)試噪聲協(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==1
[a,VS]=lpc(x_frame(:),p);
else
[a,VS]=lpc(T((r-2)*m+1:(r-2)*m+256),p);
end
%%%%%%%%%%%%%%%%幀長(zhǎng)內(nèi)過(guò)程噪聲協(xié)方差%%%%%%%%%%%%%%%%%%
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;%時(shí)間更新方程
Pn=F*P*F'+G*VS*G';
end
K=Pn*H'*(y_temp+H*P*H').^(-1); %卡爾曼增益
P=(eye(p)-K*H)*Pn;%測(cè)量更新方程
S=S+K*[x_frame(j)-H*S];
T((r-1)*m+j)=H*S;
end
%%%%%%%%%%%%%%%%對(duì)得到的每幀數(shù)據(jù)進(jìn)行加窗操作%%%%%%%%%%%%%%%%%%%%%%%%
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('時(shí)間');ylabel('幅度');title('原聲音信號(hào)');
subplot(212);plot((1:1109760)/Fs,s_out);xlabel('時(shí)間');ylabel('幅度');title('經(jīng)卡爾曼濾波后的聲音信號(hào)');
愛(ài)華網(wǎng)本文地址 » http://www.klfzs.com/a/25101016/299459.html
愛(ài)華網(wǎng)


