rls.m

数字信号处理中常用的两种算法在语音信号处理中的应用程序

%This is the RLS algorithm of two_channel acoustic echo cancellation
%clear all
tic
%-load the  impulse response-------------%
load g_1 Imp                             % one input signal in the receiving room
g1=Imp;
load g_2 Imp                             % another input signal in the receiving room
g2=Imp;
load h_1  Imp                            % one impulse response in the receiving room
h1=Imp;
load h_2  Imp                            % another impulse response in the receiving room
h2=Imp;
load h_21  Imp                           % one impulse response in the receiving room
h21=Imp;
load h_22  Imp                           % another impulse response in the receiving room 
h22=Imp;

%------read the acoustic signal----------%
% y=wavread('chinese_2.wav',[1 4000]);     % read the two channel scoustic signal
% s=y(:,1);                                % speech source signal
%%%%%%%%%----2.white gauss noise signal
s=randn(3000,1);
L1=length(s);                            % length of the input signal
N=length(g1);                            % order of the modeling filters
for i=1:L1-N+1
    s0=s(i+N-1:-1:i);
    s1(i)=s0'*g1';                       % get the signal filtered by impuse response of the transmition room
    s2(i)=s0'*g2';
end

%----preprocess the speech signal--------%
% beta=0.1;                              % factor of preprocessing
% s1=s1+beta*(s1+abs(s1));               % preprocessing of nonlinear transform
% s2=s2+beta*(s2-abs(s2));

%--generate a background noise signal----%
L=length(s1);                            % length of the useful signal
p1=sum(abs(s1).*abs(s1))/L;
p2=sum(abs(s2).*abs(s2))/L;
a=10*log10((p1+p2)/(10^4));              % background noise(40db)
noise=wgn(1,L,a);

%-------set another parameters needed----%
M=length(h1);
hh1=[h1,h2]';                            % concatenate the two room inpulse response as one vector
hh2=[h21,h22]';
w=zeros(2*M,1);                          % weight vector of the adaptive filter
d=zeros(1,L-M);                          % desired signal
y=zeros(1,L-M);                          % adaptive filter output signal
e=zeros(1,L-M);                          % defference of the y and d signal
mse=zeros(1,L-M);                        % mean square error 
mis=zeros(1,L-M);                        % misalignment 
stepsize=0.9;                            % adaptive stepsize                                                 
alpha=0.99;                              % forgetting factor           
R=alpha*eye(2*M,2*M);                    % correlation matrix
r=zeros(2*M,1);                          % cross-correlation vector

%------RLS algorithm---------------------%
for i=1:L-M
    xx1=s1(i+M-1:-1:i);                  % one receiving singal of the receiving room                       
    xx2=s2(i+M-1:-1:i);                  % another receiving singal of the receiving room
    xx=[xx1 xx2]';                       % concatenate the two room inpulse response as one vector          
    h=hh1;
%     if i>=(L-M)/2                        % alter the impulse response in the receiving room to another one
%         h=hh2;
%     end
    d(i)=h'*xx+noise(i);                 % desired signal               
    y=w'*xx;                             % filter signal
    e(i)=d(i)-y;                         % error signal
    R=(R-(R*xx)*(xx'*R)/(alpha+xx'*R*xx))/alpha;   % iterate the correlation matrix
    r=alpha*r+d(i)*xx;                   % iterate cross-correlation vector
    w=R*r;                               % iterate the coefficient of modeling filters
    mis(i)=norm(h-w)/norm(h);            % misalignment
    i
end

%--------compute the mse-----------------%
NN=length(e);
block=500;
E=[zeros(1,fix(block/2)),e,zeros(1,fix(block/2)+1)];
D=[zeros(1,fix(block/2)),d,zeros(1,fix(block/2)+1)];
for i=1:NN                               % mean square error and smoothed with 300 data
    mse(i)=E(i:i+block-1)*E(i:i+block-1)'/(D(i:i+block-1)*D(i:i+block-1)');   
end                                    

%------plot------------------------------%
figure(1);    
plot(10*log10(mis(1:NN)));                     
ylabel('misalignment');xlabel('samples');
title('Misalignment of RLS algorithm');
figure(2);
plot(10*log10(mse(1:NN)));                      
ylabel('mse');xlabel('samples');
title('Mean Square Error of RLS algorithm');

disp('norm w=');disp(norm(w));
disp('norm h=');disp(norm(h));
toc