% Partitioned block frequency domain adaptive filtering NLMS and % standard time-domain sample-based NLMS %fid=fopen('aecFar-samsung.pcm', 'rb'); % Load far end fid=fopen('aecFar.pcm', 'rb'); % Load far end %fid=fopen(farFile, 'rb'); % Load far end rrin=fread(fid,inf,'int16'); fclose(fid); %rrin=loadsl('data/far_me2.pcm'); % Load far end %fid=fopen('aecNear-samsung.pcm', 'rb'); % Load near end fid=fopen('aecNear.pcm', 'rb'); % Load near end %fid=fopen(nearFile, 'rb'); % Load near end ssin=fread(fid,inf,'int16'); %ssin = [zeros(1024,1) ; ssin(1:end-1024)]; fclose(fid); rand('state',13); fs=16000; mult=fs/8000; %rrin=rrin(fs*0+1:round(fs*120)); %ssin=ssin(fs*0+1:round(fs*120)); if fs == 8000 cohRange = 2:3; elseif fs==16000 cohRange = 2; end % Flags NLPon=1; % NLP CNon=1; % Comfort noise PLTon=1; % Plotting M = 16; % Number of partitions N = 64; % Partition length L = M*N; % Filter length if fs == 8000 mufb = 0.6; else mufb = 0.5; end %mufb=1; VADtd=48; alp = 0.1; % Power estimation factor alc = 0.1; % Coherence estimation factor beta = 0.9; % Plotting factor %% Changed a little %% step = 0.3;%0.1875; % Downward step size %% 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); %echoBandRange = ceil(1500*2/fs*N):floor(2500*2/fs*N); else echoBandRange = ceil(300*2/fs*N):floor(1800*2/fs*N); %echoBandRange = ceil(300*2/fs*N):floor(1800*2/fs*N); end %echoBandRange = ceil(1600*2/fs*N):floor(1900*2/fs*N); %echoBandRange = ceil(2000*2/fs*N):floor(4000*2/fs*N); suppState = 1; transCtr = 0; Nt=1; vt=1; ramp = 1.0003; % Upward ramp rampd = 0.999; % Downward ramp cvt = 20; % Subband VAD threshold; nnthres = 20; % Noise threshold shh=logspace(-1.3,-2.2,N+1)'; sh=[shh;flipud(shh(2:end-1))]; % Suppression profile len=length(ssin); w=zeros(L,1); % Sample-based TD NLMS WFb=zeros(N+1,M); % Block-based FD NLMS WFbOld=zeros(N+1,M); % Block-based FD NLMS YFb=zeros(N+1,M); erfb=zeros(len,1); erfb3=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; erifb=zeros(Nb+1,1)+0.1; erifb3=zeros(Nb+1,1)+0.1; ericn=zeros(Nb+1,1)+0.1; dri=zeros(Nb+1,1)+0.1; start=1; xo=zeros(N,1); do=xo; eo=xo; echoBands=zeros(Nb+1,1); cohxdAvg=zeros(Nb+1,1); cohxdSlow=zeros(Nb+1,N+1); cohedSlow=zeros(Nb+1,N+1); %overdriveM=zeros(Nb+1,N+1); cohxdFastAvg=zeros(Nb+1,1); cohxdAvgBad=zeros(Nb+1,1); cohedAvg=zeros(Nb+1,1); cohedFastAvg=zeros(Nb+1,1); hnledAvg=zeros(Nb+1,1); hnlxdAvg=zeros(Nb+1,1); ovrdV=zeros(Nb+1,1); dIdxV=zeros(Nb+1,1); SLxV=zeros(Nb+1,1); hnlSortQV=zeros(Nb+1,1); hnlPrefAvgV=zeros(Nb+1,1); mutInfAvg=zeros(Nb+1,1); %overdrive=zeros(Nb+1,1); hnled = zeros(N+1, 1); weight=zeros(N+1,1); hnlMax = zeros(N+1, 1); hnl = zeros(N+1, 1); overdrive = ones(1, N+1); xfwm=zeros(N+1,M); dfm=zeros(N+1,M); WFbD=ones(N+1,1); fbSupp = 0; 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; ovrdPos = floor((N+1)/4); ovrdSm = 2; hnlMin = 1; minCtr = 0; SeMin = 0; SdMin = 0; SeLocalAvg = 0; SeMinSm = 0; divergeFact = 1; dIdx = 1; hnlMinCtr = 0; hnlNewMin = 0; divergeState = 0; Sy=ones(N+1,1); Sym=1e7*ones(N+1,1); wins=[0;sqrt(hanning(2*N-1))]; ubufn=zeros(2*N,1); ebuf=zeros(2*N,1); ebuf2=zeros(2*N,1); ebuf4=zeros(2*N,1); mbuf=zeros(2*N,1); cohedFast = zeros(N+1,1); cohxdFast = zeros(N+1,1); cohxd = zeros(N+1,1); Se = zeros(N+1,1); Sd = zeros(N+1,1); Sx = zeros(N+1,1); SxBad = zeros(N+1,1); Sed = zeros(N+1,1); Sxd = zeros(N+1,1); SxdBad = zeros(N+1,1); hnledp=[]; cohxdMax = 0; %hh=waitbar(0,'Please wait...'); progressbar(0); %spaces = ' '; %spaces = repmat(spaces, 50, 1); %spaces = ['[' ; spaces ; ']']; %fprintf(1, spaces); %fprintf(1, '\n'); for kk=1:Nb pos = N * (kk-1) + start; % FD block method % ---------------------- Organize data xk = rrin(pos:pos+N-1); dk = ssin(pos:pos+N-1); xx = [xo;xk]; xo = xk; tmp = fft(xx); XX = tmp(1:N+1); dd = [do;dk]; % Overlap do = dk; tmp = fft(dd); % Frequency domain DD = tmp(1:N+1); % ------------------------ Power estimation pn0 = (1 - alp) * pn0 + alp * real(XX.* conj(XX)); pn = pn0; %pn = (1 - alp) * pn + alp * M * pn0; if (CNon) Yp = real(conj(DD).*DD); % Instantaneous power Sy = (1 - alp) * Sy + alp * Yp; % Averaged power mm = min(Sy,Sym); diff = Sym - mm; if (kk>50) Sym = (mm + step*diff) * ramp; % Estimated background noise power end end % ---------------------- Filtering XFm(:,1) = XX; for mm=0:(M-1) m=mm+1; YFb(:,m) = XFm(:,m) .* WFb(:,m); end yfk = sum(YFb,2); tmp = [yfk ; flipud(conj(yfk(2:N)))]; ykt = real(ifft(tmp)); ykfb = ykt(end-N+1:end); % ---------------------- Error estimation ekfb = dk - ykfb; %if sum(abs(ekfb)) < sum(abs(dk)) %ekfb = dk - ykfb; % erfb(pos:pos+N-1) = ekfb; %else %ekfb = dk; % erfb(pos:pos+N-1) = dk; %end %(kk-1)*(N*2)+1 erfb(pos:pos+N-1) = ekfb; tmp = fft([zm;ekfb]); % FD version for cancelling part (overlap-save) Ek = tmp(1:N+1); % ------------------------ Adaptation Ek2 = Ek ./(M*pn + 0.001); % Normalized error %Ek2 = Ek ./(pn + 0.001); % Normalized error %Ek2 = Ek ./(100*pn + 0.001); % Normalized error absEf = max(abs(Ek2), threshold); absEf = ones(N+1,1)*threshold./absEf; Ek2 = Ek2.*absEf; mEk = mufb.*Ek2; PP = conj(XFm).*(ones(M,1) * mEk')'; tmp = [PP ; flipud(conj(PP(2:N,:)))]; IFPP = real(ifft(tmp)); PH = IFPP(1:N,:); tmp = fft([PH;zeros(N,M)]); FPH = tmp(1:N+1,:); WFb = WFb + FPH; if mod(kk, 10*mult) == 0 WFbEn = sum(real(WFb.*conj(WFb))); %WFbEn = sum(abs(WFb)); [tmp, dIdx] = max(WFbEn); WFbD = sum(abs(WFb(:, dIdx)),2); %WFbD = WFbD / (mean(WFbD) + 1e-10); WFbD = min(max(WFbD, 0.5), 4); end dIdxV(kk) = dIdx; % NLP if (NLPon) ee = [eo;ekfb]; eo = ekfb; window = wins; if fs == 8000 %gamma = 0.88; gamma = 0.9; else %gamma = 0.92; gamma = 0.93; end %gamma = 0.9; 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); %fprintf(1,'%d: %f\n', kk, xf(4)); dfm(:,1) = df; 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)); %xRatio = real(xfwm(:,1).*conj(xfwm(:,1))) ./ ... % (real(xfwm(:,2).*conj(xfwm(:,2))) + 1e-10); %xRatio = Sx ./ (SxOld + 1e-10); %SLx = log(1/(N+1)*sum(xRatio)) - 1/(N+1)*sum(log(xRatio)); %SLxV(kk) = SLx; %freqSm = 0.9; %Sx = filter(freqSm, [1 -(1-freqSm)], Sx); %Sx(end:1) = filter(freqSm, [1 -(1-freqSm)], Sx(end:1)); %Se = filter(freqSm, [1 -(1-freqSm)], Se); %Se(end:1) = filter(freqSm, [1 -(1-freqSm)], Se(end:1)); %Sd = filter(freqSm, [1 -(1-freqSm)], Sd); %Sd(end:1) = filter(freqSm, [1 -(1-freqSm)], Sd(end:1)); %SeFast = ef.*conj(ef); %SdFast = df.*conj(df); %SxFast = xf.*conj(xf); %cohedFast = 0.9*cohedFast + 0.1*SeFast ./ (SdFast + 1e-10); %cohedFast(find(cohedFast > 1)) = 1; %cohedFast(find(cohedFast > 1)) = 1 ./ cohedFast(find(cohedFast>1)); %cohedFastAvg(kk) = mean(cohedFast(echoBandRange)); %cohedFastAvg(kk) = min(cohedFast); %cohxdFast = 0.8*cohxdFast + 0.2*log(SdFast ./ (SxFast + 1e-10)); %cohxdFastAvg(kk) = mean(cohxdFast(echoBandRange)); % coherence Sxd = gamma*Sxd + (1 - gamma)*xf.*conj(df); Sed = gamma*Sed + (1-gamma)*ef.*conj(df); %Sxd = filter(freqSm, [1 -(1-freqSm)], Sxd); %Sxd(end:1) = filter(freqSm, [1 -(1-freqSm)], Sxd(end:1)); %Sed = filter(freqSm, [1 -(1-freqSm)], Sed); %Sed(end:1) = filter(freqSm, [1 -(1-freqSm)], Sed(end:1)); cohed = real(Sed.*conj(Sed))./(Se.*Sd + 1e-10); %cohedAvg(kk) = mean(cohed(echoBandRange)); %cohedAvg(kk) = cohed(6); %cohedAvg(kk) = min(cohed); cohxd = real(Sxd.*conj(Sxd))./(Sx.*Sd + 1e-10); %freqSm = 0.5; %cohxd(3:end) = filter(freqSm, [1 -(1-freqSm)], cohxd(3:end)); %cohxd(end:3) = filter(freqSm, [1 -(1-freqSm)], cohxd(end:3)); %cohxdAvg(kk) = mean(cohxd(echoBandRange)); %cohxdAvg(kk) = (cohxd(32)); %cohxdAvg(kk) = max(cohxd); %xf = xfm(:,dIdx); %SxBad = gamma*SxBad + (1 - gamma)*real(xf.*conj(xf)); %SxdBad = gamma*SxdBad + (1 - gamma)*xf.*conj(df); %cohxdBad = real(SxdBad.*conj(SxdBad))./(SxBad.*Sd + 0.01); %cohxdAvgBad(kk) = mean(cohxdBad); %for j=1:N+1 % mutInf(j) = 0.9*mutInf(j) + 0.1*information(abs(xfm(j,:)), abs(dfm(j,:))); %end %mutInfAvg(kk) = mean(mutInf); %hnled = cohedFast; %xIdx = find(cohxd > 1 - cohed); %hnled(xIdx) = 1 - cohxd(xIdx); %hnled = 1 - max(cohxd, 1-cohedFast); hnled = min(1 - cohxd, cohed); %hnled = 1 - cohxd; %hnled = max(1 - (cohxd + (1-cohedFast)), 0); %hnled = 1 - max(cohxd, 1-cohed); if kk > 1 cohxdSlow(kk,:) = 0.99*cohxdSlow(kk-1,:) + 0.01*cohxd'; cohedSlow(kk,:) = 0.99*cohedSlow(kk-1,:) + 0.01*(1-cohed)'; end if 0 %if kk > 50 %idx = find(hnled > 0.3); hnlMax = hnlMax*0.9999; %hnlMax(idx) = max(hnlMax(idx), hnled(idx)); hnlMax = max(hnlMax, hnled); %overdrive(idx) = max(log(hnlMax(idx))/log(0.99), 1); avgHnl = mean(hnlMax(echoBandRange)); if avgHnl > 0.3 overdrive = max(log(avgHnl)/log(0.99), 1); end weight(4:end) = max(hnlMax) - hnlMax(4:end); end %[hg, gidx] = max(hnled); %fnrg = Sx(gidx) / (Sd(gidx) + 1e-10); %[tmp, bidx] = find((Sx / Sd + 1e-10) > fnrg); %hnled(bidx) = hg; %cohed1 = mean(cohed(cohRange)); % range depends on bandwidth %cohed1 = cohed1^2; %echoBands(kk) = length(find(cohed(echoBandRange) < 0.25))/length(echoBandRange); %if (fbSupp == 0) % if (echoBands(kk) > 0.8) % fbSupp = 1; % end %else % if (echoBands(kk) < 0.6) % fbSupp = 0; % end %end %overdrive(kk) = 7.5*echoBands(kk) + 0.5; % Factor by which to weight other bands %if (cohed1 < 0.1) % w = 0.8 - cohed1*10*0.4; %else % w = 0.4; %end % Weight coherence subbands %hnled = w*cohed1 + (1 - w)*cohed; %hnled = (hnled).^2; %cohed(floor(N/2):end) = cohed(floor(N/2):end).^2; %if fbSupp == 1 % cohed = zeros(size(cohed)); %end %cohed = cohed.^overdrive(kk); %hnled = gamma*hnled + (1 - gamma)*cohed; % Additional hf suppression %hnledp = [hnledp ; mean(hnled)]; %hnled(floor(N/2):end) = hnled(floor(N/2):end).^2; %ef = ef.*((weight*(min(1 - hnled)).^2 + (1 - weight).*(1 - hnled)).^2); cohedMean = mean(cohed(echoBandRange)); %aggrFact = 4*(1-mean(hnled(echoBandRange))) + 1; %[hnlSort, hnlSortIdx] = sort(hnled(echoBandRange)); [hnlSort, hnlSortIdx] = sort(1-cohxd(echoBandRange)); [xSort, xSortIdx] = sort(Sx); %aggrFact = (1-mean(hnled(echoBandRange))); %hnlSortQ = hnlSort(qIdx); hnlSortQ = mean(1 - cohxd(echoBandRange)); %hnlSortQ = mean(1 - cohxd); [hnlSort2, hnlSortIdx2] = sort(hnled(echoBandRange)); %[hnlSort2, hnlSortIdx2] = sort(hnled); hnlQuant = 0.75; hnlQuantLow = 0.5; qIdx = floor(hnlQuant*length(hnlSort2)); qIdxLow = floor(hnlQuantLow*length(hnlSort2)); hnlPrefAvg = hnlSort2(qIdx); hnlPrefAvgLow = hnlSort2(qIdxLow); %hnlPrefAvgLow = mean(hnled); %hnlPrefAvg = max(hnlSort2); %hnlPrefAvgLow = min(hnlSort2); %hnlPref = hnled(echoBandRange); %hnlPrefAvg = mean(hnlPref(xSortIdx((0.5*length(xSortIdx)):end))); %hnlPrefAvg = min(hnlPrefAvg, hnlSortQ); %hnlSortQIdx = hnlSortIdx(qIdx); %SeQ = Se(qIdx + echoBandRange(1) - 1); %SdQ = Sd(qIdx + echoBandRange(1) - 1); %SeQ = Se(qIdxLow + echoBandRange(1) - 1); %SdQ = Sd(qIdxLow + echoBandRange(1) - 1); %propLow = length(find(hnlSort < 0.1))/length(hnlSort); %aggrFact = min((1 - hnlSortQ)/2, 0.5); %aggrTerm = 1/aggrFact; %hnlg = mean(hnled(echoBandRange)); %hnlg = hnlSortQ; %if suppState == 0 % if hnlg < 0.05 % suppState = 2; % transCtr = 0; % elseif hnlg < 0.75 % suppState = 1; % transCtr = 0; % end %elseif suppState == 1 % if hnlg > 0.8 % suppState = 0; % transCtr = 0; % elseif hnlg < 0.05 % suppState = 2; % transCtr = 0; % end %else % if hnlg > 0.8 % suppState = 0; % transCtr = 0; % elseif hnlg > 0.25 % suppState = 1; % transCtr = 0; % end %end %if kk > 50 if cohedMean > 0.98 & hnlSortQ > 0.9 %if suppState == 1 % hnled = 0.5*hnled + 0.5*cohed; % %hnlSortQ = 0.5*hnlSortQ + 0.5*cohedMean; % hnlPrefAvg = 0.5*hnlPrefAvg + 0.5*cohedMean; %else % hnled = cohed; % %hnlSortQ = cohedMean; % hnlPrefAvg = cohedMean; %end suppState = 0; elseif cohedMean < 0.95 | hnlSortQ < 0.8 %if suppState == 0 % hnled = 0.5*hnled + 0.5*cohed; % %hnlSortQ = 0.5*hnlSortQ + 0.5*cohedMean; % hnlPrefAvg = 0.5*hnlPrefAvg + 0.5*cohedMean; %end suppState = 1; end if hnlSortQ < cohxdLocalMin & hnlSortQ < 0.75 cohxdLocalMin = hnlSortQ; end if cohxdLocalMin == 1 ovrd = 3; hnled = 1-cohxd; hnlPrefAvg = hnlSortQ; hnlPrefAvgLow = hnlSortQ; end if suppState == 0 hnled = cohed; hnlPrefAvg = cohedMean; hnlPrefAvgLow = cohedMean; end %if hnlPrefAvg < hnlLocalMin & hnlPrefAvg < 0.6 if hnlPrefAvgLow < hnlLocalMin & hnlPrefAvgLow < 0.6 %hnlLocalMin = hnlPrefAvg; %hnlMin = hnlPrefAvg; hnlLocalMin = hnlPrefAvgLow; hnlMin = hnlPrefAvgLow; hnlNewMin = 1; hnlMinCtr = 0; %if hnlMinCtr == 0 % hnlMinCtr = hnlMinCtr + 1; %else % hnlMinCtr = 0; % hnlMin = hnlLocalMin; %SeLocalMin = SeQ; %SdLocalMin = SdQ; %SeLocalAvg = 0; %minCtr = 0; % ovrd = max(log(0.0001)/log(hnlMin), 2); %divergeFact = hnlLocalMin; end if hnlNewMin == 1 hnlMinCtr = hnlMinCtr + 1; end if hnlMinCtr == 2 hnlNewMin = 0; hnlMinCtr = 0; %ovrd = max(log(0.0001)/log(hnlMin), 2); ovrd = max(log(0.00001)/(log(hnlMin + 1e-10) + 1e-10), 3); %ovrd = max(log(0.00000001)/(log(hnlMin + 1e-10) + 1e-10), 5); %ovrd = max(log(0.0001)/log(hnlPrefAvg), 2); %ovrd = max(log(0.001)/log(hnlMin), 2); end hnlLocalMin = min(hnlLocalMin + 0.0008/mult, 1); cohxdLocalMin = min(cohxdLocalMin + 0.0004/mult, 1); %divergeFact = hnlSortQ; %if minCtr > 0 & hnlLocalMin < 1 % hnlMin = hnlLocalMin; % %SeMin = 0.9*SeMin + 0.1*sqrt(SeLocalMin); % SdMin = sqrt(SdLocalMin); % %SeMin = sqrt(SeLocalMin)*hnlSortQ; % SeMin = sqrt(SeLocalMin); % %ovrd = log(100/SeMin)/log(hnlSortQ); % %ovrd = log(100/SeMin)/log(hnlSortQ); % ovrd = log(0.01)/log(hnlMin); % ovrd = max(ovrd, 2); % ovrdPos = hnlSortQIdx; % %ovrd = max(ovrd, 1); % %SeMin = sqrt(SeLocalAvg/5); % minCtr = 0; %else % %SeLocalMin = 0.9*SeLocalMin +0.1*SeQ; % SeLocalAvg = SeLocalAvg + SeQ; % minCtr = minCtr + 1; %end if ovrd < ovrdSm ovrdSm = 0.99*ovrdSm + 0.01*ovrd; else ovrdSm = 0.9*ovrdSm + 0.1*ovrd; end %end %ekEn = sum(real(ekfb.^2)); %dkEn = sum(real(dk.^2)); ekEn = sum(Se); dkEn = sum(Sd); if divergeState == 0 if ekEn > dkEn ef = df; divergeState = 1; %hnlPrefAvg = hnlSortQ; %hnled = (1 - cohxd); end else %if ekEn*1.1 < dkEn %if ekEn*1.26 < dkEn 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; %cohxdMaxLocal = max(cohxdSlow(kk,:)); %if kk > 50 %cohxdMaxLocal = 1-hnlSortQ; %if cohxdMaxLocal > 0.5 % %if cohxdMaxLocal > cohxdMax % odScale = max(log(cohxdMaxLocal)/log(0.95), 1); % %overdrive(7:end) = max(log(cohxdSlow(kk,7:end))/log(0.9), 1); % cohxdMax = cohxdMaxLocal; % end %end %end %cohxdMax = cohxdMax*0.999; %overdriveM(kk,:) = max(overdrive, 1); %aggrFact = 0.25; aggrFact = 0.3; %aggrFact = 0.5*propLow; %if fs == 8000 % wCurve = [0 ; 0 ; aggrFact*sqrt(linspace(0,1,N-1))' + 0.1]; %else % wCurve = [0; 0; 0; aggrFact*sqrt(linspace(0,1,N-2))' + 0.1]; %end wCurve = [0; aggrFact*sqrt(linspace(0,1,N))' + 0.1]; % For sync with C %if fs == 8000 % wCurve = wCurve(2:end); %else % wCurve = wCurve(1:end-1); %end %weight = aggrFact*(sqrt(linspace(0,1,N+1)')); %weight = aggrFact*wCurve; weight = wCurve; %weight = aggrFact*ones(N+1,1); %weight = zeros(N+1,1); %hnled = weight.*min(hnled) + (1 - weight).*hnled; %hnled = weight.*min(mean(hnled(echoBandRange)), hnled) + (1 - weight).*hnled; %hnled = weight.*min(hnlSortQ, hnled) + (1 - weight).*hnled; %hnlSortQV(kk) = mean(hnled); %hnlPrefAvgV(kk) = mean(hnled(echoBandRange)); hnled = weight.*min(hnlPrefAvg, hnled) + (1 - weight).*hnled; %od = aggrFact*(sqrt(linspace(0,1,N+1)') + aggrTerm); %od = 4*(sqrt(linspace(0,1,N+1)') + 1/4); %ovrdFact = (ovrdSm - 1) / sqrt(ovrdPos/(N+1)); %ovrdFact = ovrdSm / sqrt(echoBandRange(floor(length(echoBandRange)/2))/(N+1)); %od = ovrdFact*sqrt(linspace(0,1,N+1))' + 1; %od = ovrdSm*ones(N+1,1).*abs(WFb(:,dIdx))/(max(abs(WFb(:,dIdx)))+1e-10); %od = ovrdSm*ones(N+1,1); %od = ovrdSm*WFbD.*(sqrt(linspace(0,1,N+1))' + 1); od = ovrdSm*(sqrt(linspace(0,1,N+1))' + 1); %od = 4*(sqrt(linspace(0,1,N+1))' + 1); %od = 2*ones(N+1,1); %od = 2*ones(N+1,1); %sshift = ((1-hnled)*2-1).^3+1; sshift = ones(N+1,1); hnled = hnled.^(od.*sshift); %if hnlg > 0.75 %if (suppState ~= 0) % transCtr = 0; %end % suppState = 0; %elseif hnlg < 0.6 & hnlg > 0.2 % suppState = 1; %elseif hnlg < 0.1 %hnled = zeros(N+1, 1); %if (suppState ~= 2) % transCtr = 0; %end % suppState = 2; %else % if (suppState ~= 2) % transCtr = 0; % end % suppState = 2; %end %if suppState == 0 % hnled = ones(N+1, 1); %elseif suppState == 2 % hnled = zeros(N+1, 1); %end %hnled(find(hnled < 0.1)) = 0; %hnled = hnled.^2; %if transCtr < 5 %hnl = 0.75*hnl + 0.25*hnled; % transCtr = transCtr + 1; %else hnl = hnled; %end %hnled(find(hnled < 0.05)) = 0; ef = ef.*(hnl); %ef = ef.*(min(1 - cohxd, cohed).^2); %ef = ef.*((1-cohxd).^2); ovrdV(kk) = ovrdSm; %ovrdV(kk) = dIdx; %ovrdV(kk) = divergeFact; %hnledAvg(kk) = 1-mean(1-cohedFast(echoBandRange)); hnledAvg(kk) = 1-mean(1-cohed(echoBandRange)); hnlxdAvg(kk) = 1-mean(cohxd(echoBandRange)); %hnlxdAvg(kk) = cohxd(5); %hnlSortQV(kk) = mean(hnled); hnlSortQV(kk) = hnlPrefAvgLow; hnlPrefAvgV(kk) = hnlPrefAvg; %hnlAvg(kk) = propLow; %ef(N/2:end) = 0; %ner = (sum(Sd) ./ (sum(Se.*(hnl.^2)) + 1e-10)); % Comfort noise if (CNon) snn=sqrt(Sym); snn(1)=0; % Reject LF noise Un=snn.*exp(j*2*pi.*[0;rand(N-1,1);0]); % Weight comfort noise by suppression Un = sqrt(1-hnled.^2).*Un; Fmix = ef + Un; else Fmix = ef; end % Overlap and add in time domain for smoothness tmp = [Fmix ; flipud(conj(Fmix(2:N)))]; mixw = wins.*real(ifft(tmp)); mola = mbuf(end-N+1:end) + mixw(1:N); mbuf = mixw; ercn(pos:pos+N-1) = mola; end % NLPon % Filter update %Ek2 = Ek ./(12*pn + 0.001); % Normalized error %Ek2 = Ek2 * divergeFact; %Ek2 = Ek ./(pn + 0.001); % Normalized error %Ek2 = Ek ./(100*pn + 0.001); % Normalized error %divergeIdx = find(abs(Ek) > abs(DD)); %divergeIdx = find(Se > Sd); %threshMod = threshold*ones(N+1,1); %if length(divergeIdx) > 0 %if sum(abs(Ek)) > sum(abs(DD)) %WFb(divergeIdx,:) = WFb(divergeIdx,:) .* repmat(sqrt(Sd(divergeIdx)./(Se(divergeIdx)+1e-10))),1,M); %Ek2(divergeIdx) = Ek2(divergeIdx) .* sqrt(Sd(divergeIdx)./(Se(divergeIdx)+1e-10)); %Ek2(divergeIdx) = Ek2(divergeIdx) .* abs(DD(divergeIdx))./(abs(Ek(divergeIdx))+1e-10); %WFb(divergeIdx,:) = WFbOld(divergeIdx,:); %WFb = WFbOld; %threshMod(divergeIdx) = threshMod(divergeIdx) .* abs(DD(divergeIdx))./(abs(Ek(divergeIdx))+1e-10); % threshMod(divergeIdx) = threshMod(divergeIdx) .* sqrt(Sd(divergeIdx)./(Se(divergeIdx)+1e-10)); %end %absEf = max(abs(Ek2), threshold); %absEf = ones(N+1,1)*threshold./absEf; %absEf = max(abs(Ek2), threshMod); %absEf = threshMod./absEf; %Ek2 = Ek2.*absEf; %if sum(Se) <= sum(Sd) % mEk = mufb.*Ek2; % PP = conj(XFm).*(ones(M,1) * mEk')'; % tmp = [PP ; flipud(conj(PP(2:N,:)))]; % IFPP = real(ifft(tmp)); % PH = IFPP(1:N,:); % tmp = fft([PH;zeros(N,M)]); % FPH = tmp(1:N+1,:); % %WFbOld = WFb; % WFb = WFb + FPH; %else % WF = WFbOld; %end % Shift old FFTs %for m=M:-1:2 % XFm(:,m) = XFm(:,m-1); % YFm(:,m) = YFm(:,m-1); %end 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); %if mod(kk, floor(Nb/50)) == 0 % fprintf(1, '.'); %end if mod(kk, floor(Nb/100)) == 0 %if mod(kk, floor(Nb/500)) == 0 progressbar(kk/Nb); %figure(5) %plot(abs(WFb)); %legend('1','2','3','4','5','6','7','8','9','10','11','12'); %title(kk*N/fs); %figure(6) %plot(WFbD); %figure(6) %plot(threshMod) %if length(divergeIdx) > 0 % plot(abs(DD)) % hold on % plot(abs(Ek), 'r') % hold off %plot(min(sqrt(Sd./(Se+1e-10)),1)) %axis([0 N 0 1]); %end %figure(6) %plot(cohedFast); %axis([1 N+1 0 1]); %plot(WFbEn); %figure(7) %plot(weight); %plot([cohxd 1-cohed]); %plot([cohxd 1-cohed 1-cohedFast hnled]); %plot([cohxd cohxdFast/max(cohxdFast)]); %legend('cohxd', '1-cohed', '1-cohedFast'); %axis([1 65 0 1]); %pause(0.5); %overdrive end end progressbar(1); %figure(2); %plot([feat(:,1) feat(:,2)+1 feat(:,3)+2 mfeat+3]); %plot([feat(:,1) mfeat+1]); %figure(3); %plot(10*log10([dri erifb erifb3 ericn])); %legend('Near-end','Error','Post NLP','Final',4); % Compensate for delay %ercn=[ercn(N+1:end);zeros(N,1)]; %ercn_=[ercn_(N+1:end);zeros(N,1)]; %figure(11); %plot(cohxdSlow); %figure(12); %surf(cohxdSlow); %shading interp; %figure(13); %plot(overdriveM); %figure(14); %surf(overdriveM); %shading interp; figure(10); t = (0:Nb)*N/fs; rrinSubSamp = rrin(N*(1:(Nb+1))); plot(t, rrinSubSamp/max(abs(rrinSubSamp)),'b'); hold on plot(t, hnledAvg, 'r'); plot(t, hnlxdAvg, 'g'); plot(t, hnlSortQV, 'y'); plot(t, hnlLocalMinV, 'k'); plot(t, cohxdLocalMinV, 'c'); plot(t, hnlPrefAvgV, 'm'); %plot(t, cohxdAvg, 'r'); %plot(cohxdFastAvg, 'r'); %plot(cohxdAvgBad, 'k'); %plot(t, cohedAvg, 'k'); %plot(t, 1-cohedFastAvg, 'k'); %plot(ssin(N*(1:floor(length(ssin)/N)))/max(abs(ssin))); %plot(echoBands,'r'); %plot(overdrive, 'g'); %plot(erfb(N*(1:floor(length(erfb)/N)))/max(abs(erfb))); hold off tightx; figure(11) plot(t, ovrdV); tightx; %plot(mfeat,'r'); %plot(1-cohxyp_,'r'); %plot(Hnlxydp,'y'); %plot(hnledp,'k'); %plot(Hnlxydp, 'c'); %plot(ccohpd_,'k'); %plot(supplot_, 'g'); %plot(ones(length(mfeat),1)*rr1_, 'k'); %plot(ones(length(mfeat),1)*rr2_, 'k'); %plot(N*(1:length(feat)), feat); %plot(Sep_,'r'); %axis([1 floor(length(erfb)/N) -1 1]) %hold off %plot(10*log10([Se_, Sx_, Seu_, real(sf_.*conj(sf_))])); %legend('Se','Sx','Seu','S'); %figure(5) %plot([ercn ercn_]); figure(12) plot(t, dIdxV); %plot(t, SLxV); tightx; %figure(13) %plot(t, [ekEnV dkEnV]); %plot(t, dkEnV./(ekEnV+1e-10)); %tightx; %close(hh); %spclab(fs,ssin,erfb,ercn,'outxd.pcm'); %spclab(fs,rrin,ssin,erfb,1.78*ercn,'vqeOut-1.pcm'); %spclab(fs,erfb,'aecOutLp.pcm'); %spclab(fs,rrin,ssin,erfb,1.78*ercn,'aecOut25.pcm','vqeOut-1.pcm'); %spclab(fs,rrin,ssin,erfb,ercn,'aecOut-mba.pcm'); %spclab(fs,rrin,ssin,erfb,ercn,'aecOut.pcm'); %spclab(fs, ssin, erfb, ercn, 'out0.pcm');