www.pudn.com > NLMS.zip > NLMS.m

%比较NLMS对CCS信号的收敛性 
%   Rin 
%-------------->----------------------------->Rout 
%       X(N)    |                | 
%               |                | 
%               V                V 
%            -------          ------- 
%            |     |          |     | 
%            |     |ADAPTIVE  |     |ECHO 
%            |W(N) |FILTER    | H(N)|PATH 
%            |     |          |     | 
%            |     |          |     | 
%            -------          ------- 
%               |                | 
%               | - Dw(n)        |Dh(n) 
%      E(N)     V  +             V 
%<-------------- <---------------------------<-Sin 
%    Sout 
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ 
 
clear; 
tic; 
hold off; 
%******************** 
%Get the Xn Data 
%******************** 
fid=fopen('MySpeech.pcm','r'); 
data=fread(fid,inf,'int16'); 
Rin=data'; 
fclose(fid); 
%******************** 
%Get the ECHO data 
%******************** 
fid=fopen('MyEcho.Pcm','r'); 
data=fread(fid,inf,'int16'); 
Recho=data'; 
fclose(fid); 
%************************ 
%回声路径 
%************************ 
[pathmodel,scalor]=EC_GetModel(1,1); 
echofilter=pathmodel*scalor; 
DELAY=0*80; 
zerohead=zeros(1,DELAY); 
echofilter=[zerohead,echofilter]; 
 
%******************** 
%Processing NLMS  
%******************** 
 
N=DELAY+128;	%128是模型的最大长度 
 
NLMSW=zeros(1,N); 
X=zeros(1,N); 
NLMSStep=zeros(1,N); 
 
MIU=1;			% ALG STEP  MIU是1的效果最好 
SIGMA=0.5; 
Energ=0; 
RechoLength=length(Recho); 
DiskNLMS=zeros(1,length(Recho)); 
ERLENLMS=zeros(1,length(Recho)); 
%for(j=1:10) 
%    MIU=j*0.2; 
%    W=zeros(1,N); 
 %   X=zeros(1,N); 
 %   DiskNLMS=zeros(1,length(Recho)); 
    for(i=1:length(Recho)) 
             d=Recho(i); 
             Energ=Energ-round(X(N)^2/32768); 
             X(2:N)=X(1:N-1); 
             X(1)=Rin(i); 
             Energ=Energ+round(X(1)^2/32768);% 
             dn=round(NLMSW*X'/32768);	%32768 
             NLMSE(i)=(d-dn); 
             NLMSStep(i)=(MIU/(SIGMA+Energ)); 
             NLMSW=NLMSW+NLMSStep(i)*X*NLMSE(i);	%核心程序(2) 
             DiskNLMS(i)=10*log10(((echofilter-NLMSW/32768)*(echofilter-NLMSW/32768)')/(echofilter*echofilter')); 
    end 
%    figure(j+1); 
%    plot(DiskNLMS); 
%end 
ERLENLMS=10*log10(( NLMSE(1:RechoLength).^2+1e-20)./(Recho(1:RechoLength).^2+1e-20));%计算ERLE 
hold off; 
figure(1); 
subplot(2,1,1); 
plot(DiskNLMS,'r');xlabel('系统距离'); 
%subplot(3,1,1); 
%plot(echofilter,'r');title('回声路径 ');Xlabel('原始回声路径'); 
%subplot(3,1,2); 
%plot(W/8000);Xlabel('NLMS回声路径'); 
 
%figure(2) 
%subplot(2,1,1); 
%plot(e);title('误差信号 u=1/32768,Sigma=0.5/32768');Xlabel('NLMS'); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%processing LMS 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for(j=-5:5) 
		LMSW=zeros(1,N); 
		X=zeros(1,N); 
        %LMSStep=ones(1,length(Recho)); 
        %LMSStep=LMSStep*14/(32768*32768); 
		STEP=(j*0.02+2.6)*1e-5; 
		DiskLMS=zeros(1,length(Recho)); 
        ERLELMS=zeros(1,length(Recho)); 
		for(i=1:length(Recho)) 
            %if(i>18000) 
            %    STEP=STEP*2; 
           % end 
		    X(2:N)=X(1:N-1); 
		    X(1)=(Rin(i)); 
		    d=(Recho(i)); 
		    dn=(LMSW*X'/32768);%32768 
		    LMSE(i)=(d-dn); 
		    LMSW=LMSW+LMSE(i)*STEP*X; 
		    DiskLMS(i)=10*log10(((echofilter-LMSW/32768)*(echofilter-LMSW/32768)')/(echofilter*echofilter')); 
        end 
%        hold on; 
		figure(j+10) 
		plot(DiskLMS); 
end 
ERLELMS=10*log10((LMSE(1:RechoLength).^2+1e-20)./(Recho(1:RechoLength).^2+1e-20));%计算ERLE 
toc 
 
hold on; 
figure(1); 
%subplot(2,1,1); 
plot(DiskLMS,'g');xlabel('系统距离'); 
subplot(2,1,2); 
for i=20:20:(RechoLength) 
    ERLENLMS(i:-1:i-19)=mean(ERLENLMS(i:-1:i-19)); 
end 
plot(ERLENLMS,'r');xlabel('NLMS ERLE'); 
hold on; 
for i=20:20:(RechoLength) 
    ERLELMS(i:-1:i-19)=mean(ERLELMS(i:-1:i-19));	 
end 
plot(ERLELMS,'g');xlabel('LMS ERLE'); 
%plot(e);title('误差信号 Step=1/32768');Xlabel('LMS'); 
%figure(2); 
%hold on; 
%plot(NLMSStep,'r');axis([0,80000,0,14*3/(32768*32768)]); 
%plot(LMSStep,'g');axis([0,80000,0,14*3/(32768*32768)]); 
figure(2) 
subplot(3,1,1) 
plot(NLMSW);xlabel('NLMS'); 
subplot(3,1,2) 
plot(LMSW); xlabel('LMS'); 
subplot(3,1,3) 
plot(echofilter);xlabel('理想回声'); 
 
figure(3) 
subplot(2,1,1); 
plot(NLMSE);axis([0,15000,-2500,2500]); 
subplot(2,1,2); 
plot(LMSE);axis([0,15000,-2500,2500]); 
toc