rls.m

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

%This is the RLS algorithm of two_channel acoustic echo cancellation
% clear all
% clc
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('english.wav',[1 6000]);      % read the two channel scoustic signal
s=y(:,1);                                % speech source signal
%%%%%%%%%----2.white gauss noise signal
% load white_noise wn;
% s=wn(1:30000);

%-----filter the input signal------------%
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^5));              % background noise(40db)
noise=wgn(1,L,a);
power_n=(p1+p2)/(10^5);
%-------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                                                 
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>=20000                        % 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;
    mis(i)=10*log10(norm(h-w)/norm(h));            % misalignment
    i
end

%--------compute the mse-----------------%
NN=length(e);
block=500;
for i=1:NN-block+1                               % 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)');   
    mse(i)=10*log10(MSE(i));
    MSE2(i)=[(d(i:i+block-1)*d(i:i+block-1)')-block*power_n]/[e(i:i+block-1)*e(i:i+block-1)'-block*power_n]; 
    mse2(i)=10*log10(MSE2(i));
end                                    

%------plot------------------------------%
figure(1);    
plot(mis);                     
ylabel('misalignment');xlabel('samples');
title('Misalignment of RLS algorithm');
figure(2);
plot(mse);                      
ylabel('mse');xlabel('samples');
title('Mean Square Error of RLS algorithm');
figure(3);
plot(mse2);title('mse 2');
disp('norm w=');disp(norm(w));
disp('norm h=');disp(norm(h));

save RLS_mse mse;
save RLS_mis mis;
toc