fullaec.m

% Partitioned block frequency domain adaptive filtering NLMS and
% standard time-domain sample-based NLMS
%fid=fopen('aecFar-samsung.pcm', 'rb'); % Load far end
addpath('./aec_record');
fid=fopen('ref001.raw', '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('echo001.raw', '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=8000;
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 = 12; % 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');