kalmansignaldenoiser.m

学习卡尔曼滤波的最好事例

function output=KalmanSignalDenoiser(Noisy,Clean,fs)

% OUTPUT=KALMANSIGNALDENOISER(NOISY,CLEAN,FS)
% this purpose of this function is to demonstrate the capability of kalman
% filter for denoising noisy speech (corrupted by white noise). Kalman
% filtering of noisy speech usually have two steps: 
% 1 . Estimating the AR parameters of speech segment
% 2 . Filtering the segment
% There are different approaches for extracting AR parameters of noisy
% speech in the literature, however, in this function none is implemented.
% Clean speech signal should be provided for this purpose.
%
% ARGUMENTS
% NOISY : noise contaminated speech
% CLEAN : clean speech signal
% FS    : Sampling frequency which should be the same for both signals
%
% Output is the denoised speech signal
%
% Required functions:
% SEGMENT
%Sep-04
%Esfandiar Zavarehei

W=fix(.025*fs); %Window length is 25 ms
SP=1; %Shift percentage is 40% (10ms) %Overlap-Add method works good with this value(.4)
SpecP=13;
Window=ones(W,1);

x=segment(Clean,W,SP,Window);
y=segment(Noisy,W,SP,Window);
n=segment(Noisy-Clean,W,SP,Window);
R=var(n);

H=[zeros(1,SpecP-1) 1];
G=H'; GGT=G*H;
FUpper=[zeros(SpecP-1,1) eye(SpecP-1)];
I=eye(SpecP);

[A Q]=lpc(x,SpecP);
P=diag(repmat(R(1),1,SpecP));

o=zeros(1,W*size(x,2));% allocating memory to the output in advance save a lot of computation time
o(1:SpecP)=y(1:SpecP,1)';
hwb = waitbar(0,'Please wait...','Name','Processing');
start=SpecP+1;
Sp=o(1:SpecP)';
t=SpecP+1;
for n=1:size(x,2)
    waitbar(n/size(x,2),hwb,['Please wait... ' num2str(fix(100*n/size(x,2))) ' %'])
    F=[FUpper; fliplr(-A(n,2:end))];
    for i=start:W               
        S_=F*Sp;
        e=y(i,n)-S_(end);%innovation
        P_=F*P*F'+GGT*Q(n);    
        K=(P_*H')/(H*P_*H' + R(n));
        SOut=S_+K*e;        
        o(t-SpecP+1:t)=SOut'; %Notice that the previous SpecP-1 output samples are updated again
        P=(I-K*H)*P_;
        Sp=SOut;
        t=t+1;
    end
    start=1;
end
close(hwb)
output=o;