%RLS算法 
%用噪聲訓(xùn)練濾波器權(quán)系數(shù)矢量，再進(jìn)行濾波
%dxt
%EI，HUST
%2008.12
clear all                      %清除內(nèi)存變量                      
clc                            %清除command框
%讀取混合有噪聲的語(yǔ)音信號(hào)
[xinit,Fs,nbits] = wavread('signalnoise.wav');
xinit = xinit';                 %轉(zhuǎn)置
xn  = xinit(1,:);               %取混有噪聲的語(yǔ)音信號(hào)的單聲道
n  = wavread('noise.wav');      %讀取噪聲信號(hào)
n  = n';                        %轉(zhuǎn)置
n  = n(1,:);                    %取噪聲信號(hào)的單聲道
  
N=length(n);					%噪聲信號(hào)的長(zhǎng)度
M=32;                           %濾波器的階數(shù)
L=0.98;                        %遺忘因子
Ns = length(xn);                %有噪聲的語(yǔ)音信號(hào)
d = 0.002*ones(1,N);             %期望輸出
x1 = 10*n+d;                       %RLS濾波器輸入信號(hào)

w=zeros(1,M);                   %權(quán)系數(shù)
x=zeros(1,M);                   %1*M，代表x(n)'
u=zeros(1,N);                   %臨時(shí)數(shù)組u=x(n)*T(n-1)*x(n)'
T=200*eye(M);                   %T(-1)
k=zeros(M,1);                   %增益k(n)
e=zeros(1,N);                   %e(n|n-1)
e1=e;                           %平方誤差
w1 = zeros(1,Ns-M+1);
for n=M:N
    x=x1(n:-1:n-M+1);           %x(n)
    u(n)=x*T*x';
    k=T*x'/(L+u(n));            %計(jì)算增益k(n)
    e(n)=d(n)-x*w';             %計(jì)算e(n|n-1)
    w=w+k'*e(n);                %計(jì)算權(quán)系數(shù)w(n)
    w1(N-n+1)=w(1);
    T=(T-x*k*T)/L;              %計(jì)算T(n)
    e1(n)=e(n)*e(n);
end
%利用訓(xùn)練的濾波器系數(shù)進(jìn)行語(yǔ)音信號(hào)的濾波    
s = zeros(1,Ns-M+1);
for n=M:Ns
   s(n-M+1) = 0;
   for j=1:M
       s(n-M+1)= s(n-M+1)+w(j)*xn(n-j+1);
   end
end
 
subplot(2,1,1);                   %繪制原始語(yǔ)音信號(hào)
plot(xn);
title('混有噪聲的原始語(yǔ)音信號(hào)');
subplot(2,1,2);                   %繪制濾波后的信號(hào)
plot(s,'k');
title('濾波后語(yǔ)音信號(hào)');
wavwrite(s,Fs,nbits,'filter.wav');
wavwrite(5*s,Fs,nbits,'filter5.wav');