基于RLS算法的语音噪声回声消除处理MATLAB代码

clear; clc; Ts=0.1; num=100;%?采样点数? k=1:num; x=10*cos(2*pi*k*Ts+0.25*pi);%正弦干扰 %s=0.5*randn(1,num); %x=x+s; xr=cos(2*pi*k*Ts);%参考信号 p=2;%滤波器长度 xita=zeros(p,num-p+1); %保存权值 err=zeros(1,num-p+1); %保存误差 kn=zeros(2,num-p+1); lamda=0.99; %遗忘因子 diag_I=[1,0;0,1]; for i=1:num-p+1 if i==1 hn=[xr(1,i);0]; corr=10^5*diag_I; kn(:,1)=(corr*hn)./(lamda^i+(hn')*corr*hn);%是列向量 err(1,i)=x(1,i)-xita(1,1)*xr(1,i); xita(:,i)=err(1,i).*kn(:,i); corr=(diag_I-kn(:,i)*hn')*corr; aaa=1; else hn=[xr(1,i); xr(1,i-1)]; kn(:,i)=(corr*hn)./(lamda^i+(hn')*corr*hn);%是列向量 % fprintf( '%d',(lamda^i+(hn')*corr*hn)); err(1,i)=x(1,i)-xita(:,i-1)'*hn; fprintf('xi=%d tmp=%d\n',x(1,i),xita(:,i-1)'*hn); disp(err(1,i)); xita(:,i)=xita(:,i-1)+err(1,i).*kn(:,i); corr=(diag_I-kn(:,i)*hn')*corr; disp(corr); end end subplot(2,1,2); plot(k(1,1:num-p+1),xita,'r-'); subplot(2,1,1); plot(k(1,1:num-p+1),err,'b-'); sysorder=50;%抽头数?? [inp,sr]=audioread('huiyin2.wav');%输入含回声的信号?? inp=inp(:); [d1,sr] = audioread('2.wav');%输入原始信号? d=[d1;zeros(4000,1)]; %需要使得d和inp行数完全一致才能运算 d=d(:); totallength=size(d,1);%步长，返回d的行数?? N=length(inp); %k=0:1382493;?? k=0:691246; %算法开始?? w=zeros(sysorder,1);%初始化?? w1=zeros(sysorder,1);%初始化?? error=zeros(N,1); %初始化?? EE=zeros(N,1); Loop=15000; RMSE=zeros(1,476000); for n=sysorder:N %u=inp(n:-1:n-sysorder+1);? ?? ?? ? %u的矩阵?? u=inp(n-sysorder+1:1:n);%u矩阵为输入的信号中的一小段?? %NLMS算法?? y(n)=w'*u; %系统输出?? r(n)=u'*u; %自相关矩阵?? e(n)=d(n)-y(n);%误差?? srr(n)=10*log10(d(n)/(inp(n)-d(n))); srr1(n)=10*log10(d(n)/(y(n)-d(n))); %e(n)=inp(n)-y(n);? ?? ?? ?? ?? ???%误差? ?? ?? ?? fai=1e-10;%修正参数，防止r(n)过小导致步长值太大?? if n<2000 mu=0.32; else mu=0.15; end w=0.8*w+mu*u*e(n)/(r(n)+fai);%NLMS迭代方程?? www(:,:,n)=w;%显示每一步的W,但由于太多而无法显示 Ai=mean(w); % RMSE(1,n)=sqrt(sum((w-Ai).^2)/50) ; % error(n)=error(n)+e(n)^2;?? %LMS算法?? y1(n)=w1'*u; e1(n)=d(n)-y1(n); w1=0.8*w1+0.1*u*e1(n);%LMS迭代方程?? end EE=EE+error; error=EE/Loop; error=error.^2; error(n)=error'*error; y=y(:); y1=y1(:); e=e(:); e1=e1(:); srr=srr(:); srr1=srr1(:); error=10*log10(error); figure(1); subplot(5,1,1); plot(inp); title('系统输出波形对比'); xlabel('样本n'); ylabel('混响9秒信号波形'); subplot(5,1,2); plot(y,'r'); xlabel('样本n'); ylabel('NLMS输出信号波形'); subplot(5,1,3); plot(y1,'g'); axis([0 476000 -1 1]); xlabel('样本n'); ylabel('LMS输出信号波形'); subplot(5,1,4); plot(d); xlabel('样本n'); ylabel('原始信号波形'); figure; plot(error); figure; plot(e) %sound(inp,sr); %sound(d,sr) %sound(y,sr);?? %audiowrite(nlms.wav,y,sr);? ? sound(y1,sr); %audiowrite(lms.wav,y1,sr);?