webrtc aec matlab.zip_WebRTC_AEC_webrtc aec matlab_webrtc aec_

% Partitioned block frequency domain adaptive filtering NLMS and % standard time-domain sample-based NLMS [ssin,fs1]=audioread('near.wav','native'); [rrin,fs2]=audioread('far.wav','native'); ssin=double(ssin); rrin=double(rrin); fs=8000; mult=fs/8000;%%使用理论未知 % Flags NLPon=1; % NLP on M = 16; % Number of partitions N = 64; % Partition length L = M*N; % Filter length if fs == 8000 mufb = 0.6; else mufb = 0.8;%%步长，位置其值如何设置 end alp = 0.15; % Power estimation factor? alc = 0.1; % Coherence estimation factor %% Changed a little %% %% if fs == 8000 threshold=2e-6; % DTrob threshold else %threshold=0.7e-6; threshold=1.5e-6; %%门限 end if fs == 8000 echoBandRange = ceil(300*2/fs*N):floor(1800*2/fs*N); else echoBandRange = ceil(60*2/fs*N):floor(1500*2/fs*N);%%回声频谱范围 end suppState = 1;%%回声抑制标志 len=length(ssin); WFb=zeros(N+1,M); % Block-based FD(frequency domain) NLMS YFb=zeros(N+1,M); erfb=zeros(len,1); ercn=zeros(len,1);%%回声消除后的结果 zm=zeros(N,1); XFm=zeros(N+1,M); YFm=zeros(N+1,M); pn0=10*ones(N+1,1); pn=zeros(N+1,1); NN=len; Nb=floor(NN/N)-M; start=1; xo=zeros(N,1); do=xo; eo=xo; cohxdAvg=zeros(Nb+1,1);%%程序中有用到，不知什么意思 cohxdSlow=zeros(Nb+1,N+1);%%算法中未使用 cohedSlow=zeros(Nb+1,N+1);%%同上 cohedAvg=zeros(Nb+1,1);%%算法未使用，貌似用来记录ed相关性 hnledAvg=zeros(Nb+1,1); hnlxdAvg=zeros(Nb+1,1); ovrdV=zeros(Nb+1,1); dIdxV=zeros(Nb+1,1); hnlSortQV=zeros(Nb+1,1); hnlPrefAvgV=zeros(Nb+1,1); hnled = zeros(N+1, 1); weight=zeros(N+1,1); hnl = zeros(N+1, 1); xfwm=zeros(N+1,M); dfm=zeros(N+1,M); WFbD=ones(N+1,1); hnlLocalMin = 1; cohxdLocalMin = 1; hnlLocalMinV=zeros(Nb+1,1); cohxdLocalMinV=zeros(Nb+1,1); hnlMinV=zeros(Nb+1,1); dkEnV=zeros(Nb+1,1); ekEnV=zeros(Nb+1,1); ovrd = 2; ovrdSm = 2; hnlMin = 1; dIdx = 1; divergeState = 0; Sym=1e7*ones(N+1,1); wins=[0;sqrt(hanning(2*N-1))]; mbuf=zeros(2*N,1); cohxd = zeros(N+1,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); hh=waitbar(0,'Please wait...'); %%变量初始化 for kk=1:Nb %%循环每一帧，每帧64点 pos = N * (kk-1) + start;%%位置记录 % FD block method % ---------------------- Organize data %far is speaker played music xk = rrin(pos:pos+N-1);%%远处声源64点数据 %near is micphone captured signal dk = ssin(pos:pos+N-1);%%近处麦克风捕捉64点数据 %----------------------- far end signal process xx = [xo;xk];%%xo为上帧远端数据，初始为64*1的零向量 xo = xk; tmp = fft(xx); XX = tmp(1:N+1); dd = [do;dk]; % Overlap %%do为上帧远端数据，初始为64*1的零向量， do = dk; tmp = fft(dd); % Frequency domain DD = tmp(1:N+1);%%对xx及dd进行FFT变换，xx的频谱用于计算xx的功率，用于NLMS公式 % ------------------------far end Power estimation pn0 = (1 - alp) * pn0 + alp * real(XX.* conj(XX)); pn = pn0;%%XX功率计算 % ---------------------- Filtering XFm(:,1) = XX;%%XFm为65*16矩阵，存储XX的频谱，用于计算回声估计，2-16列为前15帧数据 YFb = XFm .* WFb; %%YFb为回声信号频谱 yfk = sum(YFb,2);%%将YFb按行求和，得到65*1的矩阵 tmp = [yfk ; flipud(conj(yfk(2:N)))]; ykt = real(ifft(tmp)); ykfb = ykt(end-N+1:end);%%得到回声信号时域，为64*1矩阵 % ---------------------- Error estimation ekfb = dk - ykfb;%%时域误差信号 erfb(pos:pos+N-1) = ekfb;%%误差信号保存，作图使用 tmp = fft([zm;ekfb]); % FD version for cancelling part (overlap-save) Ek = tmp(1:N+1);%%得到误差信号频谱，65*1矩阵 % ------------------------ Adaptation Ek2 = Ek ./(pn + 0.001); % Normalized error%%归一化误差信号 absEf = max(abs(Ek2), threshold);%%为65*1矩阵，每一点最小为门限 absEf = ones(N+1,1)*threshold./absEf; Ek2 = Ek2.*absEf;%%对误差频谱进行门限限制，使之最大为门限值 mEk = mufb.*Ek2; PP = conj(XFm).*(ones(M,1) * mEk')';%%系数调节大小，误差*步长*参考信号，为系数频谱 WFb = WFb + PP;%%更新权重 WFbEn = sum(real(WFb.*conj(WFb))); [tmp, dIdx] = max(WFbEn);%%存储16块中权重最大的位置及大小 WFbD = sum(abs(WFb(:, dIdx)),2);%%未使用 WFbD = min(max(WFbD, 0.5), 4);%%未使用 dIdxV(kk) = dIdx;%%记录每次最大块位置，算法未使用 % NLP %%非线性过程，重点部分 if (NLPon) ee = [eo;ekfb];%%重叠分段FFT处理，为了与XX,DD相匹配以计算相关性等 eo = ekfb; window = wins; if fs == 8000 gamma = 0.9; else gamma = 0.93; end tmp = fft(xx.*window);%%xx为远端信号 xf = tmp(1:N+1); tmp = fft(dd.*window);%%dd为近端mic捕捉信号 df = tmp(1:N+1); tmp = fft(ee.*window);%%ee为估计的近端声源信号 ef = tmp(1:N+1);%%对EE,XX,DD加窗进行FFT变换，加窗为了防止频谱泄露 xfwm(:,1) = xf;%%将本次XF的数据传入XFWM，其2-16为前15帧数据，为65*16矩阵 xf = xfwm(:,dIdx);%%得到最大能量列，认为其为本次处理的远端参考信号 dfm(:,1) = df;%%近端mic捕捉信号 SxOld = Sx;%%算法未使用，可屏蔽 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));%%计算xf,df,ef的功率谱密度，与上次结果平滑 % coherence 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(cohxd(echoBandRange));%%！！！到这里到这里不懂是用来干啥的 hnled = min(1 - cohxd, cohed);%%此处为了最大化回声抑制，二者均是相关性愈大，回声越小 [hnlSort2, hnlSortIdx2] = sort(hnled(echoBandRange)); cohedMean = mean(cohed(echoBandRange));%%计算所选点均值 NcodxdMean = mean(1 - cohxd(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 = 0; elseif cohedMean < 0.95 || NcodxdMean < 0.8 suppState = 1; end%%条件一下回声小，条件二下回声大，做抑制，其他情况不做处理 if NcodxdMean < cohxdLocalMin && NcodxdMean < 0.75 cohxdLocalMin = NcodxdMean;%%cohxdlocalmin，满足条件时赋值 end if cohxdLocalMin == 1%%此处表示未能更新coxdlocalmin，回声较小 ovrd = 3; hnled = 1-cohxd; hnlPrefAvg = NcodxdMean; hnlPrefAvgLow = NcodxdMean; end if suppState == 0%%不需要回声消除，为hnled，hnlprefavg等赋值 hnled = cohed; hnlPrefAvg = cohedMean; hnlPrefAvgLow = cohedMean; end if hnlPrefAvgLow < hnlLocalMin && hnlPrefAvgLow < 0.6 hnlLocalMin = hnlPrefAvgLow; hnlMin = hnlPrefAvgLow; ovrd = max(log(0.00000001)/(log(hnlMin + 1e-10) + 1e-10), 5);%%不懂更新规则，不知该值作用 end hnlLocalMin = min(hnlLocalMin + 0.0008/mult, 1);%%更新的意义不懂 cohxdLocalMin = min(cohxdLocalMin + 0.0004/mult, 1); if ovrd < ovrdSm ovrdSm = 0.99*ovrdSm + 0.01*ovrd; else ovrdSm = 0.9*ovrdSm + 0.1*ovrd; end ekEn = sum(Se); dkEn = sum(Sd); if divergeState == 0%%用于发散处理 if ekEn > dkEn ef = df; divergeState = 1; end else if ekEn*1.05 < dkEn divergeState = 0; else ef = df; end end if ekEn > dkEn*19.95%%异常，重置系数矩阵 WFb=zeros(N+1,M); % Block-based FD NLMS end ekEnV(kk) = ekEn; dkEnV(kk) = dkEn; hnlLocalMinV(kk) = hnlLocalMin; cohxdLocalMinV(kk) = cohxdLocalMin; hnlMinV(kk) = hnlMin; %%平滑滤波器系数及抑制指数 aggrFact = 0.3; wCurve = [0; aggrFact*sqrt(linspace(0,1,N))' + 0.1]; weight = wCurve; hnled = weight.*min(hnlPrefAvg, hnled) + (1 - weight).*hnled; od = ovrdSm*(sqrt(linspace(0,1,N+1))' + 1); hnled = hnled.^(od.*ones(N+1,1)); hnl = hnled; ef = ef.*(hnl); ovrdV(kk) = ovrdSm; hnledAvg(kk) = 1-mean(1-cohed(echoBandRange)); hnlxdAvg(kk) = 1-mean(cohxd(echoBandRange)); hnlS