aec.rar_aec_matlab 实现webrtc_webrtc_回声消除_回声消除aec

clear;clc;close all; adap=1; if(adap==1) [ssin,fs1]=audioread('near1.wav','native'); % [ssin,fs1]=audioread('speak_female_echo.wav','native'); [rrin,fs2]=audioread('speak_female.wav','native'); % [rrin,fs2]=audioread('far2.wav','native'); [echo,fs_echo]=audioread('speak_female_echo.wav','native'); elseif(adap==0) [ssin,fs1]=audioread('near.wav','native'); [rrin,fs2]=audioread('speak_female.wav','native'); end ssin=double(ssin); rrin=double(rrin); fs=fs1; l1=length(ssin); l2=length(rrin); if(l1<=l2) ssin=[ssin;zeros(l2-l1,1)]; len=l2; else rrin=[rrin;zeros(l1-l2,1)]; len=l1; end len_echo=length(echo); if(len_echo<=len) echo=[echo;zeros(len-len_echo,1)]; end echo1=zeros(len,1); if fs == 8000 threshold=2e-6; else %threshold=0.7e-6; threshold=1.5e-6; end M=16; N=64; L=M*N; Nb=floor(len/N); mufb=0.5; xo=zeros(N,1); do=zeros(N,1); zo=zeros(N,1); eo=zo; zm=zeros(N,1); Se = zeros(N+1,1); Sd = zeros(N+1,1); Sx = zeros(N+1,1); Sed = zeros(N+1,1); Sxd = zeros(N+1,1); XFm=zeros(N+1,M); xfwm=zeros(N+1,M); WFb=zeros(N+1,M); Ek2En=zeros(Nb,1); ekfbEn=zeros(Nb,1); Pn=zeros(N+1,1); am=zeros(N+1,1); mbuf=zeros(2*N,1); W=zeros(Nb,M); alp=0.15; adap_on=1; gamma=0.9; hnlLocalMin = 1; cohxdLocalMin = 1; mult=fs/8000; ovrd = 2; ovrdSm = 2; divergeState = 0; wins=[0;sqrt(hanning(2*N-1))]; cohedAvg=zeros(Nb,1); cohxdAvg=zeros(Nb,1); echoBandRange = ceil(500/fs*N):floor(3500/fs*N); suppState=0; uek=zeros(Nb,1); j=1; aa=0.98; eta= 0.15; mu=0.98; c=sqrt(pi)/2; qk=0.3; qkr=(1-qk)/qk; ksi_min=10^(-25/10); img=sqrt(-1); noise_mean=zeros(N+1,1); win=hamming(N+1); for kk=1:Nb pos=N*(kk-1)+1; xk = rrin(pos:pos+N-1); dk = ssin(pos:pos+N-1); xx=[xo;xk]; dd=[do;dk]; xo=xk; do=dk; tmp=fft(xx); XX=tmp(1:N+1); tmp=fft(dd); DD=tmp(1:N+1); %flitering XFm(:,1)=XX; YFb=XFm.*WFb; yfk = sum(YFb,2); tmp = [yfk ; flipud(conj(yfk(2:N)))]; ykt = real(ifft(tmp)); ykfb = ykt(end-N+1:end); ekfb = dk - ykfb; ekfbEn(kk)=sqrt(sum(ekfb.*conj(ekfb))); echo1(pos:pos+N-1)=ykfb; tmp = fft(xx.*wins); XX = tmp(1:N+1); %Adapting Pn=(1-alp)*Pn+alp*XX.*conj(XX); tmp = fft([zm;ekfb]); Ek = tmp(1:N+1); Ek2 = Ek ./(Pn + 0.001); absEf = max(abs(Ek2), threshold); absEf = ones(N+1,1)*threshold./absEf; Ek2 = Ek2.*absEf; Ek2En(kk)=sum(sqrt(Ek2.*conj(Ek2))); mEk = mufb.*Ek2; PP = conj(XFm).*(ones(M,1) * mEk')'; WFb = WFb + adap_on*PP; WFbEn = sum(real(WFb.*conj(WFb))); W(kk,:)=WFbEn; % [~, dIdx] = max(WFbEn); %NLP ee = [eo;ekfb]; eo = ekfb; window = wins; tmp = fft(xx.*window); xf = tmp(1:N+1); tmp = fft(dd.*window); df = tmp(1:N+1); tmp = fft(ee.*window); ef = tmp(1:N+1); % xfwm(:,1) = xf; % xf = xfwm(:,dIdx); Se = gamma*Se + (1-gamma)*real(ef.*conj(ef)); Sd = gamma*Sd + (1-gamma)*real(df.*conj(df)); Sx = gamma*Sx + (1-gamma)*real(xf.*conj(xf)); Sxd = gamma*Sxd + (1-gamma)*xf.*conj(df); Sed = gamma*Sed + (1-gamma)*ef.*conj(df); cohed = real(Sed.*conj(Sed))./(Se.*Sd + 1e-10); cohedAvg(kk) = mean(cohed(echoBandRange)); cohxd = real(Sxd.*conj(Sxd))./(Sx.*Sd + 1e-10); cohxdAvg(kk) = mean(1 - cohxd(echoBandRange)); cohedMean = mean(cohed(echoBandRange)); NcodxdMean = mean(1 - cohxd(echoBandRange)); hnled = min(1 - cohxd, cohed); [hnlSort2, hnlSortIdx2] = sort(hnled(echoBandRange)); hnlQuant = 0.75; hnlQuantLow = 0.5; qIdx = floor(hnlQuant*length(hnlSort2)); qIdxLow = floor(hnlQuantLow*length(hnlSort2)); hnlPrefAvg = hnlSort2(qIdx); hnlPrefAvgLow = hnlSort2(qIdxLow); if cohedMean > 0.98 && NcodxdMean > 0.9 suppState = 1; elseif cohedMean < 0.95 && NcodxdMean < 0.8 suppState = 0; end if NcodxdMean < cohxdLocalMin && NcodxdMean < 0.75 cohxdLocalMin = NcodxdMean; else cohxdLocalMin=min(cohxdAvg(kk-1)+0.0006*mult,1); end if cohxdLocalMin == 1 || suppState==1 uek(kk)=0; else uek(kk)=1; end if uek(kk)==0 if suppState == 1 && cohxdLocalMin==1 hnled = cohed; hnlPrefAvg = cohedMean; hnlPrefAvgLow = cohedMean; end else hnled = 1-cohxd; hnlPrefAvg = NcodxdMean; hnlPrefAvgLow = NcodxdMean; end if hnlPrefAvgLow < hnlLocalMin && hnlPrefAvgLow < 0.6 hnlLocalMin = hnlPrefAvgLow; hnlMin = hnlPrefAvgLow; else hnlMin=min(hnlLocalMin + 0.0008/mult, 1); end ovrd = max((-11.5)/(log(hnlMin + 1e-10) + 1e-10), 2); if ovrd < ovrdSm ovrdSm = 0.99*ovrdSm + 0.01*ovrd; else ovrdSm = 0.9*ovrdSm + 0.1*ovrd; end if cohxdLocalMin == 1 ovrdSm = 2; hnled = 1-cohxd; hnlPrefAvg = NcodxdMean; hnlPrefAvgLow = NcodxdMean; end ekEn = sum(Se); dkEn = sum(Sd); if divergeState == 0 if ekEn > 5*dkEn ef = df; divergeState = 1; end else if ekEn < 3*dkEn divergeState = 0; else ef = df; end end if ekEn > dkEn*19.95 WFb=zeros(N+1,M); end aggrFact = 0.3; wCurve = [0; aggrFact*sqrt(linspace(0,1,N))' + 0.1]; weight = wCurve; for i=1:length(hnled) if hnlPrefAvg<hnled(i) hnled(i) = weight(i)*hnlPrefAvg + (1 - weight(i))*hnled(i); end end od = ovrdSm*(sqrt(linspace(0,1,N+1))' + 1); hnled = hnled.^(od); ef = ef.*(hnled); tmp = [ef ; flipud(conj(ef(2:N)))]; mixw = real(ifft(tmp)); mola = mbuf(end-N+1:end) + mixw(1:N); mbuf = mixw; ercn(pos:pos+N-1) = mola; XFm(:,2:end) = XFm(:,1:end-1); % YFm(:,2:end) = YFm(:,1:end-1); % xfwm(:,2:end) = xfwm(:,1:end-1); % dfm(:,2:end) = dfm(:,1:end-1); end audiowrite('source.wav',int16(ercn),fs2); figure; subplot(411); plot(ssin); title('near') subplot(412); plot(rrin); title('far') subplot(413); plot(ercn); title('e') [source,fs4]=audioread('speak_male.wav','native'); subplot(414); plot(source); title('speak source') tt1=1:N:len_echo; echo=double(echo); tt2=1:N:Nb*N; figure(2) plot(tt2,W); title('WFbEn') figure(4) plot(tt2,Ek2En); title('Ek2En') figure(3) subplot(311); plot(tt1,cohedAvg(1:length(tt1)),'r'); title('cohedAvg'); hold on; plot(echo/max(echo)); subplot(312); plot(tt1,cohxdAvg(1:length(tt1)),'r'); title('cohxdAvg'); hold on; plot(echo/max(echo)); subplot(313); plot(tt1,uek(1:length(tt1)),'r'); hold on; plot(echo/max(echo)); title('uek'); erle=((echo-echo1).*(echo-echo1))./(echo.*echo+0.00001); erle=10*log(erle); tt3=1:len; figure(4) plot(tt3,erle);