447 lines
12 KiB
Matlab
447 lines
12 KiB
Matlab
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); |