min_norm_direct.m

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

% This is the minimum norm direct solution of two_channel acoustic echo cancellation.
%clear all
%clc
   function y=min_norm_direct();
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;

%------generate the input signal---------%
%%%%%%%%%----1.two channel acoustic signal
y=wavread('chinese_2.wav',[1 20000]);
s=y(:,1);                                % select one channel of the two
%%%%%%%%%----2.white gauss noise signal
% s=randn(40000,1);
%%%%%%%%%----3.aotoregresive(AR) signal
% n=0.001*randn(40000,1);
% den=[1 0.6 -0.7 -0.5];
% num=1;
% s=filter(num,den,n);

%-----filter the source signal with impulse response---%
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 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    minimum norm direct solution BEGIN 
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 impulse 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);                 % compute the desired signal
    amp=xx'*xx;                          % get the square of the norm (or energy) of the input samples
    damp=d(i)/amp;                 
    w=xx*damp;                           % compute the weight value of the adaptive filter
    y(i)=w'*xx;                          % calculate the adaptive filter output signal
    e(i)=d(i)-y(i);                      % error of the two output signal
    mis(i)=10*log10(norm(h-w)/norm(h));  % misalignment
end 
%    minimum norm direct solution END    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-------compute the mean square error----%
NN=length(e);
block=500;                               % blocksize used to smooth mse
E=[zeros(1,fix(block/2)),e,zeros(1,fix(block/2))];
D=[zeros(1,fix(block/2)),d,zeros(1,fix(block/2))];
for i=1:NN                               % mean square error and smoothed with 500 data
    mse(i)=E(i:i+block-1)*E(i:i+block-1)'/(D(i:i+block-1)*D(i:i+block-1)'); 
    %mse2(i)=E(i:i+block-1)*E(i:i+block-1)'/block;
    MSE(i)=10*log10(mse(i));
    %MSE2(i)=10*log10(mse2(i));
end  

%----------plot--------------------------%
figure(1);
plot(mis,'r');grid on;
ylabel('misalignment (dB)');xlabel('samples');
title('Misalignment of MNDS algorithm');
figure(2);
plot(1:NN,MSE,'r'); axis([0 L1 -400 0]);
grid on;  
ylabel('mse (dB)');xlabel('samples');
title('Mean Square Error of MNDS algorithm');

disp('norm w=');disp(norm(w));
disp('norm h=');disp(norm(hh1));
disp('norm h=');disp(norm(hh2));

toc