webrtc AEC 线性滤波 PBFDAF（均匀分块频域自适应滤波）介绍

计算一个脉冲响应和输入信号的卷积，除了使用原始的时域卷积以外，还有如下方法：

1. FFT卷积的方法：对输入信号（长度M）和脉冲响应（长度N）分别补零到K（K>M+N-1)，然后分别计算FFT，然后相乘，最后反FFT，取前M+N-1个元素作为最终的卷积结果
2. 输入信号很长时，将输入信号分成一帧一帧，使用overlap-add或者overlap-save的方法
3. 当脉冲信号和输入信号都很长时，可使用均匀分块卷积

均匀分块卷积

        均匀分块卷积与频域自适应滤波(FDAF)结合，就是WebRTC AEC中线性滤波处理中的算法核心。

在介绍PBFDAF之前，我们来看一下均匀分块卷积的计算流程图：

我们分几个部分讲解上图的计算流程：

1、脉冲响应分块

        如上图红色矩形部分，将脉冲响应分成P个等长的不重叠的小块，每小块的长度为B，我们把这些小块叫做子滤波器（filter part 1,2...P)，将每个小块后面补B个零，分别构成2B长度的序列，然后进行实数FFT。我们知道实数序列的FFT结果具有对称性，因此实数FFT返回B+1个点（类似numpy的rfft.fft）

2、将输入信号分块

        如上图红色线框内的部分，将输入信号分成不重叠的等长的分块（帧），分块长度为B，通过一个buffer构造重叠长度为B，这样输入buffer的长度为2B，将输入buffer进行实数FFT，得到B+1个复数结果。然后将fft结果存入一个长度为P的队列，队列进口处存储最新的输入buffer fft结果，旧的输入buffer的fft结果从队列的出口扔掉。这个队列有个名字叫Frequency-domain delay line。

3、频域相乘相加和反变换

        第三部分如上图红色矩形内，将第一部分准备的P个分块脉冲响应的FFT结果分别与第二部分形成的输入buffer fft结果的队列分别相乘，然后沿着P的方向求和。得到B+1长度的FFT结果，然后在进行复数到实数的IFFT（如numpy.rfft.ifft)，输出是2B个实数序列，取后B个元素作为输入block对于的输出。

WebRTC AEC中的分块卷积

    % 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);

    
    % ----------------------   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); 

xk是近端麦克风输入信号，要对近端信号进行线性滤波，得到估计的回声信号。

xx就是输入buffer，xk是输入的N个样本点，xo是上一次的输入N个样本点。对输入buffer进行傅里叶变换的到XX，将XX存入XFm,XFm就是频域的那个队列

然后将队列中各个输入buffer的fft结果与WFb进行相乘相加。WFb就是存放脉冲响应分块傅里叶变换的结果，因为这是自适应滤波，所以WFb矩阵的初始值的元素全部是0。M与上文中的P对应，N与上文中的B对应。

WebRTC AEC中的PBFDAF

% 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));

% Flags
NLPon=0;  % NLP
CNon=0; % 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;  

alp = 0.1; % Power estimation factor alc = 0.1; % Coherence estimation factor



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);

zm=zeros(N,1);
XFm=zeros(N+1,M);

pn0=10*ones(N+1,1);
pn=zeros(N+1,1);
NN=len;
Nb=floor(NN/N)-M;

start=1;
xo=zeros(N,1);
do=xo;
eo=xo;

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);
	
    
    % ------------------------  Power estimation
    pn0 = (1 - alp) * pn0 + alp * real(XX.* conj(XX));
    pn = pn0;
    %pn = (1 - alp) * pn + alp * M * pn0;
    
    % ----------------------   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; 
    
    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
	
	
	absEf = max(abs(Ek2), threshold);
	absEf = ones(N+1,1)*threshold./absEf;
	Ek2 = Ek2.*absEf; % 让EK2限定到threshold
	
	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;
    
    
	
    XFm(:,2:end) = XFm(:,1:end-1);
    
	
end

audiowrite('aec_out.wav', erfb/32768, fs);

图片：

Frank Wefers. Partitioned convolution algorithms for real-time auralization