function [emicrophone,aaa]=compsup(microphone,TheFarEnd,avtime,samplingfreq);
% microphone = microphone signal
% aaa = nonlinearity input variable
% TheFarEnd = far end signal
% avtime = interval to compute suppression from (seconds)
% samplingfreq = sampling frequency
%if(nargin==6)
% fprintf(1,'suppress has received a delay sequence\n');
%end
Ap500=[ 1.00, -4.95, 9.801, -9.70299, 4.80298005, -0.9509900499];
Bp500=[ 0.662743088639636, -2.5841655608125, 3.77668102146288, -2.45182477425154, 0.596566274575251, 0.0];
Ap200=[ 1.00, -4.875, 9.50625, -9.26859375, 4.518439453125, -0.881095693359375];
Bp200=[ 0.862545460994275, -3.2832804496114, 4.67892032308828, -2.95798023879133, 0.699796870041299, 0.0];
maxDelay=0.4; %[s]
histLen=1; %[s]
% CONSTANTS THAT YOU CAN EXPERIMENT WITH
A_GAIN=10.0; % for the suppress case
oversampling = 2; % must be power of 2; minimum is 2; 4 works
% fine for support=64, but for support=128,
% 8 gives better results.
support=64; %512 % fft support (frequency resolution; at low
% settings you can hear more distortion
% (e.g. pitch that is left-over from far-end))
% 128 works well, 64 is ok)
lowlevel = mean(abs(microphone))*0.0001;
G_ol = 0; % Use overlapping sets of estimates
% ECHO SUPPRESSION SPECIFIC PARAMETERS
suppress_overdrive=1.0; % overdrive factor for suppression 1.4 is good
gamma_echo=1.0; % same as suppress_overdrive but at different place
de_echo_bound=0.0;
mLim=10; % rank of matrix G
%limBW = 1; % use bandwidth-limited response for G
if mLim > (support/2+1)
error('mLim in suppress.m too large\n');
end
dynrange=1.0000e-004;
% other, constants
hsupport = support/2;
hsupport1 = hsupport+1;
factor = 2 / oversampling;
updatel = support/oversampling;
win=sqrt(designwindow(0,support));
estLen = round(avtime * samplingfreq/updatel)
runningfmean =0.0;
mLim = floor(hsupport1/2);
V = sqrt(2/hsupport1)*cos(pi/hsupport1*(repmat((0:hsupport1-1) + 0.5, mLim, 1).* ...
repmat((0:mLim-1)' + 0.5, 1, hsupport1)));
fprintf(1,'updatel is %5.3f s\n', updatel/samplingfreq);
bandfirst=8; bandlast=25;
dosmooth=0; % to get rid of wavy bin counts (can be worse or better)
% compute some constants
blockLen = support/oversampling;
maxDelayb = floor(samplingfreq*maxDelay/updatel); % in blocks
histLenb = floor(samplingfreq*histLen/updatel); % in blocks
x0=TheFarEnd;
y0=microphone;
%input
tlength=min([length(microphone),length(TheFarEnd)]);
updateno=floor(tlength/updatel);
tlength=updatel*updateno;
updateno = updateno - oversampling + 1;
TheFarEnd =TheFarEnd(1:tlength);
microphone =microphone(1:tlength);
TheFarEnd =[zeros(hsupport,1);TheFarEnd(1:tlength)];
microphone =[zeros(hsupport,1);microphone(1:tlength)];
% signal length
n = min([floor(length(x0)/support)*support,floor(length(y0)/support)*support]);
nb = n/blockLen - oversampling + 1; % in blocks
% initialize space
win = sqrt([0 ; hanning(support-1)]);
sxAll2 = zeros(hsupport1,nb);
syAll2 = zeros(hsupport1,nb);
z500=zeros(5,maxDelayb+1);
z200=zeros(5,hsupport1);
bxspectrum=uint32(zeros(nb,1));
bxhist=uint32(zeros(maxDelayb+1,1));
byspectrum=uint32(zeros(nb,1));
bcount=zeros(1+maxDelayb,nb);
fcount=zeros(1+maxDelayb,nb);
fout=zeros(1+maxDelayb,nb);
delay=zeros(nb,1);
tdelay=zeros(nb,1);
nlgains=zeros(nb,1);
% create space (mainly for debugging)
emicrophone=zeros(tlength,1);
femicrophone=complex(zeros(hsupport1,updateno));
thefilter=zeros(hsupport1,updateno);
thelimiter=ones(hsupport1,updateno);
fTheFarEnd=complex(zeros(hsupport1,updateno));
afTheFarEnd=zeros(hsupport1,updateno);
fmicrophone=complex(zeros(hsupport1,updateno));
afmicrophone=zeros(hsupport1,updateno);
G = zeros(hsupport1, hsupport1);
zerovec = zeros(hsupport1,1);
zeromat = zeros(hsupport1);
% Reset sums
mmxs_a = zerovec;
mmys_a = zerovec;
s2xs_a = zerovec;
s2ys_a = zerovec;
Rxxs_a = zeromat;
Ryxs_a = zeromat;
count_a = 1;
mmxs_b = zerovec;
mmys_b = zerovec;
s2xs_b = zerovec;
s2ys_b = zerovec;
Rxxs_b = zeromat;
Ryxs_b = zeromat;
count_b = 1;
nog=0;
aaa=zeros(size(TheFarEnd));
% loop over signal blocks
fprintf(1,'.. Suppression; averaging G over %5.1f seconds; file length %5.1f seconds ..\n',avtime, length(microphone)/samplingfreq);
fprintf(1,'.. SUPPRESSING ONLY AFTER %5.1f SECONDS! ..\n',avtime);
fprintf(1,'.. 20 seconds is good ..\n');
hh = waitbar_j(0,'Please wait...');
for i=1:updateno
sb = (i-1)*updatel + 1;
se=sb+support-1;
% analysis FFTs
temp=fft(win .* TheFarEnd(sb:se));
fTheFarEnd(:,i)=temp(1:hsupport1);
xf=fTheFarEnd(:,i);
afTheFarEnd(:,i)= abs(fTheFarEnd(:,i));
temp=win .* microphone(sb:se);
temp=fft(win .* microphone(sb:se));
fmicrophone(:,i)=temp(1:hsupport1);
yf=fmicrophone(:,i);
afmicrophone(:,i)= abs(fmicrophone(:,i));
ener_orig = afmicrophone(:,i)'*afmicrophone(:,i);
if( ener_orig == 0)
afmicrophone(:,i)=lowlevel*ones(size(afmicrophone(:,i)));
end
% use log domain (showed improved performance)
xxf= sqrt(real(xf.*conj(xf))+1e-20);
yyf= sqrt(real(yf.*conj(yf))+1e-20);
sxAll2(:,i) = 20*log10(xxf);
syAll2(:,i) = 20*log10(yyf);
mD=min(i-1,maxDelayb);
xthreshold = sum(sxAll2(:,i-mD:i),2)/(maxDelayb+1);
[yout, z200] = filter(Bp200,Ap200,syAll2(:,i),z200,2);
yout=yout/(maxDelayb+1);
ythreshold = mean(syAll2(:,i-mD:i),2);
bxspectrum(i)=getBspectrum(sxAll2(:,i),xthreshold,bandfirst,bandlast);
byspectrum(i)=getBspectrum(syAll2(:,i),yout,bandfirst,bandlast);
bxhist(end-mD:end)=bxspectrum(i-mD:i);
bcount(:,i)=hisser2( ...
byspectrum(i),flipud(bxhist),bandfirst,bandlast);
[fout(:,i), z500] = filter(Bp500,Ap500,bcount(:,i),z500,2);
fcount(:,i)=sum(bcount(:,max(1,i-histLenb+1):i),2); % using the history range
fout(:,i)=round(fout(:,i));
[value,delay(i)]=min(fout(:,i),[],1);
tdelay(i)=(delay(i)-1)*support/(samplingfreq*oversampling);
% compensate
idel = max(i - delay(i) + 1,1);
% echo suppression
noisyspec = afmicrophone(:,i);
% Estimate G using covariance matrices
% Cumulative estimates
xx = afTheFarEnd(:,idel);
yy = afmicrophone(:,i);
% Means
mmxs_a = mmxs_a + xx;
mmys_a = mmys_a + yy;
if (G_ol)
mmxs_b = mmxs_b + xx;
mmys_b = mmys_b + yy;
mmy = mean([mmys_a/count_a mmys_b/count_b],2);
mmx = mean([mmxs_a/count_a mmxs_b/count_b],2);
else
mmx = mmxs_a/count_a;
mmy = mmys_a/count_a;
end
count_a = count_a + 1;
count_b = count_b + 1;
% Mean removal
xxm = xx - mmx;
yym = yy - mmy;
% Variances
s2xs_a = s2xs_a + xxm .* xxm;
s2ys_a = s2ys_a + yym .* yym;
s2xs_b = s2xs_b + xxm .* xxm;
s2ys_b = s2ys_b + yym .* yym;
% Correlation matrices
Rxxs_a = Rxxs_a + xxm * xxm';
Ryxs_a = Ryxs_a + yym * xxm';
Rxxs_b = Rxxs_b + xxm * xxm';
Ryxs_b = Ryxs_b + yym * xxm';
% Gain matrix A
if mod(i, estLen) == 0
% Cumulative based estimates
Rxxf = Rxxs_a / (estLen - 1);
Ryxf = Ryxs_a / (estLen - 1);
% Variance normalization
s2x2 = s2xs_a / (estLen - 1);
s2x2 = sqrt(s2x2);
% Sx = diag(max(s2x2,dynrange*max(s2x2)));
Sx = diag(s2x2);
if (sum(s2x2) > 0)
iSx = inv(Sx);
else
iSx= Sx + 0.01;
end
s2y2 = s2ys_a / (estLen - 1);
s2y2 = sqrt(s2y2);
% Sy = diag(max(s2y2,dynrange*max(s2y2)));
Sy = diag(s2y2);
iSy = inv(Sy);
rx = iSx * Rxxf * iSx;
ryx = iSy * Ryxf * iSx;
dbd= 7; % Us less than the full matrix
% k x m
% Bandlimited structure on G
LSEon = 0; % Default is using MMSE
if (LSEon)
ryx = ryx*rx;
rx = rx*rx;
end
p = dbd-1;
gaj = min(min(hsupport1,2*p+1),min([p+(1:hsupport1); hsupport1+p+1-(1:hsupport1)]));
cgaj = [0 cumsum(gaj)];
G3 = zeros(hsupport1);
for kk=1:hsupport1
ki = max(0,kk-p-1);
if (sum(sum(rx(ki+1:ki+gaj(kk),ki+1:ki+gaj(kk))))>0)
G3(kk,ki+1:ki+gaj(kk)) = ryx(kk,ki+1:ki+gaj(kk))/rx(ki+1:ki+gaj(kk),ki+1:ki+gaj(kk));
else
G3(kk,ki+1:ki+gaj(kk)) = ryx(kk,ki+1:ki+gaj(kk));
end
end
% End Bandlimited structure
G = G3;
G(abs(G)<0.01)=0;
G = suppress_overdrive * Sy * G * iSx;
if 1
figure(32); mi=2;
surf(max(min(G,mi),-mi)); view(2)
title('Unscaled Masked Limited-bandwidth G');
end
pause(0.05);
% Reset sums
mmxs_a = zerovec;
mmys_a = zerovec;
s2xs_a = zerovec;
s2ys_a = zerovec;
Rxxs_a = zeromat;
Ryxs_a = zeromat;
count_a = 1;
end
if (G_ol)
% Gain matrix B
if ((mod((i-estLen/2), estLen) == 0) & i>estLen)
% Cumulative based estimates
Rxxf = Rxxs_b / (estLen - 1);
Ryxf = Ryxs_b / (estLen - 1);
% Variance normalization
s2x2 = s2xs_b / (estLen - 1);
s2x2 = sqrt(s2x2);
Sx = diag(max(s2x2,dynrange*max(s2x2)));
iSx = inv(Sx);
s2y2 = s2ys_b / (estLen - 1);
s2y2 = sqrt(s2y2);
Sy = diag(max(s2y2,dynrange*max(s2y2)));
iSy = inv(Sy);
rx = iSx * Rxxf * iSx;
ryx = iSy * Ryxf * iSx;
% Bandlimited structure on G
LSEon = 0; % Default is using MMSE
if (LSEon)
ryx = ryx*rx;
rx = rx*rx;
end
p = dbd-1;
gaj = min(min(hsupport1,2*p+1),min([p+(1:hsupport1); hsupport1+p+1-(1:hsupport1)]));
cgaj = [0 cumsum(gaj)];
G3 = zeros(hsupport1);
for kk=1:hsupport1
ki = max(0,kk-p-1);
G3(kk,ki+1:ki+gaj(kk)) = ryx(kk,ki+1:ki+gaj(kk))/rx(ki+1:ki+gaj(kk),ki+1:ki+gaj(kk));
end
% End Bandlimited structure
G = G3;
G(abs(G)<0.01)=0;
G = suppress_overdrive * Sy * G * iSx;
if 1
figure(32); mi=2;
surf(max(min(G,mi),-mi)); view(2)
title('Unscaled Masked Limited-bandwidth G');
end
pause(0.05);
% Reset sums
mmxs_b = zerovec;
mmys_b = zerovec;
s2xs_b = zerovec;
s2ys_b = zerovec;
Rxxs_b = zeromat;
Ryxs_b = zeromat;
count_b = 1;
end
end
FECestimate2 = G*afTheFarEnd(:,idel);
% compute Wiener filter and suppressor function
thefilter(:,i) = (noisyspec - gamma_echo*FECestimate2) ./ noisyspec;
ix0 = find(thefilter(:,i)<de_echo_bound); % bounding trick 1
thefilter(ix0,i) = de_echo_bound; % bounding trick 2
ix0 = find(thefilter(:,i)>1); % bounding in reasonable range
thefilter(ix0,i) = 1;
% NONLINEARITY
nl_alpha=0.8; % memory; seems not very critical
nlSeverity=0.3; % nonlinearity severity: 0 does nothing; 1 suppresses all
thefmean=mean(thefilter(8:16,i));
if (thefmean<1)
disp('');
end
runningfmean = nl_alpha*runningfmean + (1-nl_alpha)*thefmean;
aaa(sb+20+1:sb+20+updatel)=10000*runningfmean* ones(updatel,1); % debug
slope0=1.0/(1.0-nlSeverity); %
thegain = max(0.0,min(1.0,slope0*(runningfmean-nlSeverity)));
% END NONLINEARITY
thefilter(:,i) = thegain*thefilter(:,i);
% Wiener filtering
femicrophone(:,i) = fmicrophone(:,i) .* thefilter(:,i);
thelimiter(:,i) = (noisyspec - A_GAIN*FECestimate2) ./ noisyspec;
index = find(thelimiter(:,i)>1.0);
thelimiter(index,i) = 1.0;
index = find(thelimiter(:,i)<0.0);
thelimiter(index,i) = 0.0;
if (rem(i,floor(updateno/20))==0)
fprintf(1,'.');
end
if mod(i,50)==0
waitbar_j(i/updateno,hh);
end
% reconstruction; first make spectrum odd
temp=[femicrophone(:,i);flipud(conj(femicrophone(2:hsupport,i)))];
emicrophone(sb:se) = emicrophone(sb:se) + factor * win .* real(ifft(temp));
end
fprintf(1,'\n');
close(hh);