%% Acoustic Echo Cancellation in Real-Time
% Develop a real-time Acoustic Echo Canceller for real-time full-band communications.
%% Introduction
% Hands-free full-duplex devices have one or many microphones and loudspeakers that are always active, which can produce an undesirable feedback. Acoustic Echo Cancellation (AEC) is an algorithm that suppresses the portion of the microphone signal that is produced by the loudspeaker and it's echo in the room (shown as the yellow paths in the figure below).
%
% This example uses subband filtering to processes full band signals from real recordings and can also be executed in real-time. For a first introduction to AEC or for a demonstration of the frequency domain adaptive filter method (FDAF) using a simulated room response, you may also look at our other AEC example. This example can be used as a tutorial on how to design a real-world AEC system.
%
%% Room and Device Characterization
% To design a real-world AEC system, it may be useful to first measure the impulse response (IR) of the device in a representative room. This IR can be used to help us design our adaptive filter, as it provides information such as the initial delay (hardware, software and acoustic latency) and decay length. Having these measured IR's also make it easier to test and refine our algorithm offline, while reducing the risk that our algorithm would not perform as well when executing live, in real-time.
% If the loudspeaker and microphone are connected to the computer through a full-duplex audio device, we can simply use the impulseResponseMeasurer app to measure and export (save) an IR. It is also possible to measure the IR from an offline recording using impzest as shown below. But first, the following figure shows how the impulseResponseMeasurer app would be used.
%
% If the system cannot present itself as a full-duplex audio interface, but can play and record audio, the same result can be achieved using impzest. To do so, produce test signal using sweeptone, play and record it on your system, and then provide these two signals to impzest.

fs = 48e3;
swtone = sweeptone(1.5,1.5,fs,'SweepFrequencyRange',[10 fs/2],'ExcitationLevel',-24.1);
tonerep3 = repmat(swtone,3,1);
audiowrite('tonerep3.wav',tonerep3,fs);

% Here, play tonerep3.wav on the system and record the tonerec.wav file. Load it and measure the impulse response of the system.

tonerec = audioread('tonerec.wav');
ir = impzest(swtone, tonerec);

% We can plot the IR in the same way the app does.

plot((1:size(ir,1))/fs, ir)
title('Impulse Response of Device in Office Room')
xlabel('Time (seconds)')
ylabel('Amplitude')
ylim([-.45 .7])
grid on

% In addition to providing latency and decay-length information, the measured IR can be used to produce a large amount of test files for the AEC algorithm. In addition to the loudspeaker to microphone IR, multiple IR's can also be measured for several interlocuter positions in the room. Alternatively, we can play and record the near and far-end interlocuter recordings to obtain our test files. This is demonstrated in the next section.
%
%% Recording Audio for Offline Use
% When designing an AEC system, much of the work of designing and refining the algorithms can be achieved offline, i.e. without the hardware in the loop. Although the desired result is to demonstrate the AEC algorithm live in real-time, that is not a practical way to design and measure the performance of our system. Here, using one microphone, we will record two loudspeakers. One loudspeaker is used to play the far-side interlocuter and corresponds to the AEC device loudspeaker, and the other loudspeaker is playing the near-end interlocuter (corresponding to an actual interlocuter in the room). This method also allows us to produce test recordings with a controlled amount of double-talk and a known target for our double-cross detector. This way, we can design our adaptive filters independently of the double-talk detector, which makes both tasks that much easier.
% The end goal of the AEC algorithm will be to remove the far-end speaker from these recording, with only the knowledge of the far-end signal before it went through the device's loudspeaker and reverberated in the room.
% Our first step is to load two short sentences, one for the far-end speaker and one for the near-end speaker. We will then produce a target "detector" and longer near- and far-end sequences that both include period of silence at different times to control the amount of double-talk. For our test, the talkers will not overlap at first, and then will overlap towards the end. This is illustrated in the figure below.

% Load speech file recordings (without reverb)
[xa,spfs] = audioread('FH47_07.wav');
[xb,smfs] = audioread('MG39_01.wav');
assert(spfs==smfs);
fs = spfs;
xz = zeros(1e5,1);

% Define the 'oracle' that matches the above overlaps
% (close to ideal detection of when to adapt filter)
da = ones(size(xa),'logical');
db = ones(size(xb),'logical');
dz = zeros(size(xz),'logical');
detectA = [da;dz;da;da;da;dz];
detectB = [dz;db;dz;dz;db];
detector = and(detectA, not([detectB;dz(1:length(detectA)-length(detectB))]));

% Create some near and far-end talking with more or less overlapping
xa = [xa;xz;xa;xa;xa;xz];
xb = [xz;xb;xz;xz;xb];
xb = [xb;xz(1:length(xa)-length(xb))];

aecref = [xa xb];
% We can save aecref with the following code:
% audiowrite('aecref.wav',aecref,fs);

% Plot the near and far-end talkers
t = (1:size(xa,1))/fs;
subplot(2,1,1)
plot(t,xa)
title('far-end')
ylim([-1 1])
grid on
subplot(2,1,2)
plot(t,xb)
title('near-end')
xlabel('Time (seconds)')
ylim([-1 1])
grid on

% Next, using audioPlayerRecorder, we shall play this and record the microphone. We check for dropped samples as we may want to redo the recording if there were any loss frames. Finally, we save the recording for offline use.

frameSize = 512; % 512/48 = 10.7 ms latency
abOut = dsp.AsyncBuffer(size(aecref,1)+frameSize);
abInp = dsp.AsyncBuffer(size(aecref,1));
write(abInp,aecref);

apr = audioPlayerRecorder('SampleRate',fs,'Device','ASIO MADIface USB','BitDepth','16-bit integer',...
                          'PlayerChannelMapping',[1 2],'RecorderChannelMapping',9);
                      
droppedFlag = false;
while abInp.NumUnreadSamples > 0
    % Play from file and store record audio
    audioToBePlayed = read(abInp,frameSize);
    [recorded, nUnderruns, nOverruns] = apr(audioToBePlayed);
    write(abOut,recorded);
    if ~droppedFlag
        droppedFlag = nUnderruns+nOverruns > 0;
    end
end

% Report any loss of synchronization
if droppedFlag
    fprintf('Audio queue dropped samples.\n');
end

% Save the audio recording
aecrec = read(abOut);
% We can save aecrec with the following code:
% audiowrite('aecrec.wav',aecrec,fs);

release(abInp);
release(abOut);
release(apr);

figure
plot((1:size(aecrec,1))/fs, aecrec)
title('In-Room Recording of Near and Far-End Interlocuters')
xlabel('Time (seconds)')
ylabel('Amplitude')
grid on

%
%% Algorithm Design
% Here, we will first design our algorithm offline before attempting to test it in real-time, using the actual recording we just made. Alternatively, we could have used the room impulse responses that we've measured instead.
% The AEC algorithm here will use an adaptive filter (NLMS) to characterize the response of the loudspeaker-room-microphone system. Thus, the error signal of the NLMS filter correspond to the sounds in the room that did not come from the loudspeaker. Since we are using full-band (48 kHz) signals and we want to minimize the lengths of the individual NLMS filters, we design a four-band filter bank so that we now need four NLMS filters, but with a less coefficients each. This filter bank design is encapsulated in the QuadFilterBank object and it is based on firpr2chfb.
% For this first step, we are using an "oracle" for the filter adaptation, meaning that we are using our knowledge of when the near-end talker is active instead of developing a double-talk detector while we are trying to fine-tune our adaptive filters. Once our adaptive filters are well behaved, we will replace the oracle by our own detector and we will be able to compare it to the oracle.

% Load the reference and recording
[aecref,~]  = audioread('aecref.wav');
[aecrec,fs] = audioread('aecmic.wav');

% Initalize buffers
maxlen = max(size(aecref,1),size(aecrec,1));
mic = dsp.AsyncBuffer(maxlen);   
spk = dsp.AsyncBuffer(maxlen+frameSize); 
yad = dsp.AsyncBuffer(maxlen+frameSize);
det = dsp.AsyncBuffer(size(aecrec,1));

write(mic,aecrec);
write(spk,aecref(:,1));
write(det,detector);

% Four-band filterbank
L = 96;             % QMF Filter length
B = 4;              % Number of filterbands
stepsize = 0.3;     % LMS filter step size

qfb1 = QuadFilterBank(L,L*2/3);
qfb2 = QuadFilterBank(L,L*2/3);

% Create LMS Filters (NLMS)
NLms = 2000; % Lms filter lengths
lms = cell(1,B);
for ii = 1:B
    lms{ii} = dsp.LMSFilter(NLms, 'Method', 'Normalized LMS', 'StepSizeSource', 'Input port');
end

%% Processing Loop

while spk.NumUnreadSamples > 0
   
    xa = read(spk,frameSize);
    xm = read(mic,frameSize);
    
    xfa = analyze(qfb1,xa);
    xfm = analyze(qfb2,xm);
    
    % Oracle detection
    detect = read(det,frameSize);
    outPower = sum(xfa(:,1).^2)/frameSize*B;
    adaptFrame = min(detect) && outPower>1e-6;
    thisstep = adaptFrame*stepsize;
        
    % Adaptive filters
    efb = zeros(frameSize/B,B); 
    for ii = 1:B
        % Filter from output signal (far end) to microphone (echo + near end).
        % The output error (after adaptation) corresponds to the near end speech.
        [~,efb(:,ii)] = lms{ii}(xfa(:,ii),xfm(:,ii),thisstep);
    end

    ylms = synthesize(qfb1,efb);
    write(yad, ylms);
    
end

release(apr);

% Display results

ylmsout = read(yad);

plot(aecref(:,1))
hold on
plot(aecref(:,2))
plot(aecrec,':')
plot(ylmsout)
plot(detector*stepsize)
hold off
legend('far','near','recording','output','step size')

% Compare the output obtained from the AEC algorithm to the desired signal. We observe that after approximately half a second, the amplitude of the undesired signal is very low. This can be further improved by adding so-called "non-linear processing" (NLP). NLP consists of attenuating the residual output signal when we detect that there is no active near-side interlocuter talking.

plot(aecref(:,2))
title('Desired Output')

% In this next step, we will replace the "oracle" detector by our own algorithm. We will also add a basic NLP that will reduce the output level if no near-end speech is detected while there is far-end speech. Also, since adaptive filters perform better with "flat" (i.e. "white") signals, we will apply a frequency pre-emphasis to the signals that will enter the NLMS filter, and we apply the corresponding frequency de-emphasis to the error (i.e. output) signal. Another small optimization added here is to use smaller filter lengths for the higher-frequency bands. This was done because we can pass the IR through the QuadFilterBank and observe the outputs and observe that the decay length is different in each band.
% To design a basic double-talk detector, we filter the far-end signal with the NLMS filter and measure the similarity with the microphone signal. We optimize the threshold in the first moments because we know for a fact that the NLMS filter is not adapted yet.
% With all these changes, we observe that the output signal is visibly improved (in particular between the 2 and 4 second marks).

% Initialization
reset(qfb1);
reset(qfb2);
reset(mic);
reset(spk);

write(mic,aecrec);
write(spk,aecref(:,1));

% Pre-Emphasis on xa and xir (whitening filter)
spkEmph = dsp.IIRFilter('Numerator',[1 -.8],'Denominator',[1 .97]);
lmsDemp = dsp.IIRFilter('Numerator',[1 .97],'Denominator',[1 -.8]);
micEmph = clone(spkEmph);

myDetect = dsp.AsyncBuffer(maxlen+frameSize);

lmsFilt = dsp.FIRFilter('NumeratorSource','Input port');

% Create LMS Filters (NLMS)
% Use different LMS filter lengths for each band
NLms = [2000, 1250, 1250, 1000];
lms = cell(1,B);
for ii = 1:B
    lms{ii} = dsp.LMSFilter(NLms(ii), 'Method', 'Normalized LMS', 'StepSizeSource', 'Input port');
end

%% Processing Loop

wts = zeros(NLms(1),B);
yfb = zeros(frameSize/B,B);
efb = zeros(frameSize/B,B); 

thres = 1;
oldGain = 1;
crossLength = 50;
init = 10;
ca = false;

while spk.NumUnreadSamples > 0
   
    xa = spkEmph(read(spk,frameSize));
    xm = micEmph(read(mic,frameSize));
    
    xfa = analyze(qfb1,xa);
    xfm = analyze(qfb2,xm);
    
    % Step size control (only ON/OFF for now)
    % Use prediction gain to detect double-talk
    pred = lmsFilt(xfa(:,1),wts(:,1).');
    outPower = sum(xfa(:,1).^2)/frameSize*B;
    perr = xfm(:,1) - pred;
    predGain = sum(xfm(:,1).^2)/sum(perr.^2);
    prevca = ca;
    th2 = 1e-8;
    if init > 0
        ca = outPower>th2;
        init = init - 1;
    else
        ca = (predGain>thres && outPower>th2);
    end
    if ca % change the threshold while adapting
        thres = min(2,thres*1.0025);
    end
    % Ratio between predicted feedback and microphone input
    thisstep = stepsize*ca;
    write(myDetect,repmat(thisstep,frameSize,1));
    % Adaptive filters
    for ii = 1:B
        % Filter from output signal (far end) to microphone (echo + near end).
        % The output error (after adaptation) corresponds to the near end speech.
        [yfb(:,ii),efb(:,ii),wts(1:NLms(ii),ii)] = lms{ii}(xfa(:,ii),xfm(:,ii),thisstep);
        if ca && prevca % non-linear processing
            newGain = (1-min(1,predGain/200));
            efb(:,ii) = efb(:,ii).*[linspace(oldGain,newGain,crossLength).' ; newGain*ones(frameSize/B-crossLength,1)];
            oldGain = newGain;
        end
    end

    ylms = synthesize(qfb1,efb);

    write(yad, lmsDemp(ylms));
    
end

% Display results
ylmsout = read(yad);
myAdapt = read(myDetect);

plot(aecref(:,1))
hold on
plot(aecref(:,2))
plot(aecrec,':')
plot(ylmsout)
plot(myAdapt,'LineWidth',2)
hold off
legend('far','near','recording','output','step size')

% Finally, we could replace mic by an audio device to try the algorithm in real-time.
% Now, we can listen to the microphone recording and the AEC algorithm output.

sound(ylmsout,fs)

% We could also listen to the microphone capture before AEC processing by listening to aecrec.
%
%% Real-Time Demonstration
% Once satisfied with the AEC algorithm, we can make slight modifications to make it a real-time demonstration. We simply need to replace the "near-end" signal by silence, and we use the full-duplex audio device to play the far-end audio and record the microphone. The helperAcousticEchoCancelRealTime function implements this, and can be found in this example's directory.
% Copyright 2019 The MathWorks, Inc.