pseudo.m

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

% This is the minimum norm solution algotithm of two_channel acoustic echo cancellation.
clear all
tic
%-load the useful signal and 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;

load g1_128 Imp                             % one input signal in the receiving room
g1=Imp;
load g2_128 Imp                             % another input signal in the receiving room
g2=Imp;
load h1_128  Imp                            % one impulse response in the receiving room
h1=Imp;
load h2_128  Imp                            % another impulse response in the receiving room 
h2=Imp;

%------read the acoustic signal----------%
y=wavread('chinese_2.wav',[1 1000]);
s=y(:,1);                                % speech source 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

%----ilitiate another prameters----------%
L=length(s1);
M=length(h1);                            % order of the adaptive filters
R=zeros(2*M,2*M);                        % correlation matrix
w=zeros(2*M,1);                          % initializing  the modeling filters               
h=[h1 h2]';                              % concatenate the two room inpulse response as one vector 
F=zeros(M,M);
B=zeros(2*M,2*M);
I=eye(2);
alpha=0.9;                               % forgetting factor
beta=0.001;                              % stepsize

for i=1:M
    for r=1:M
        F(i,r)=exp(-j*2*pi*(i-1)*(r-1)/M);
    end
end
B=kron(I,F);
BH=conj(B');

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
    x=[xx1 xx2]';                        % concatenate the two room inpulse response as one vector     
    d(i)=h'*x;                           % desired signal               
    y=w'*x;                              % filter signal        
    e(i)=d(i)-y;                         % error signal
    R=alpha*R+(1-alpha)*x*x';            % iterate the correlation matrix
    P=BH*R*B/M;                          % 利用DFT矩阵将相关矩阵对角化
    r=rank(P);                           % rank
    [U,V]=eig(P,'nobalance');            % 特征值分解
    
    V0=zeros(r,r);
    U0=zeros(2*M,r);
    t=1;
    for n=1:N
        if abs(V(n,n))~=0
            V0(t,t)=V(n,n);
            U0(:,t)=U(:,n);
            t=t+1;
        end
    end
    
    for k=1:r
        V0(k,k)=1/V0(k,k);       
    end
    
    P0=U0*V0*conj(U0');     
    R0=B*P0*BH/M;                        % 极小范数伪逆
    w=w+beta*R0*x*e(i);                  % iterate the coefficient of modeling filters
    mis(i)=norm(h-w)/norm(h);            % misalignment
    i
end

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)'); 
    MSE(i)=10*log10(mse(i));
end  
                          
figure(1);    
plot(10*log10(mis));                     % plot misalignment curve                   
ylabel('mis');xlabel('samples');
title('Misalignment of pseudo inverse matrix');
figure(2);
plot(MSE);                               % plot misalignment curve                
ylabel('mse');xlabel('samples');
title('Mean Square Error of pseudo inverse matrix');

toc