13 asptrcpbfdaf

Purpose

Block filtering and coefficient update in frequency domain using the Reduced Complexity Partitioned Block Frequency Domain (RCPBFDAF) algorithm.

Syntax

[W,X,x,y,e,Px,ci,w]=asptrcpbfdaf(M,x,xn,dn,X,W,mu,n,c,b,Px,ci)

Description

asptrcpbfdaf() is a reduced complexity and extended version of asptpbfdaf(). The computational complexity reduction is achieved by constraining one or more partition each call to asptrcpbfdaf() instead of constraining all partitions as in the case of asptpbfdaf(). This saves two FFT operations for each skipped partition while keeping the performance almost unaffected. No reduction in complexity is achieved for unconstrained filters. asptrcpbfdaf() is also designed to accommodate the general case of $M = g * L$, where $M$ is the partition length, $L$ is the block length, and $g$ is an integer, which allows choosing $L$ less than $M$ to further reduce the processing delay. Similar to asptpbfdaf(), the adaptive filter is splitted into $P$ partitions, from which only the first partition introduces a delay of $L$ samples, and therefore, reducing the processing delay by a factor of $P$ compared to the BFDAF for the same filter length. asptrcpbfdaf() performs filtering and coefficient update in the frequency domain using the overlap-save method, and therefore, provides an efficient implementation for long adaptive filters in applications such as acoustic echo cancelers where the adaptive filter can be a few thousand coefficients long. asptrcpbfdaf() is a block processing algorithm, every call processes $L$ input and $L$ desired samples ($L$ is the block length), to produce $L$ filter output samples and $L$ error samples, besides updating all filter coefficients in the frequency domain. The parameters of asptrcpbfdaf() are summarized below ( see Fig. 4.2 ).

Input Parameters [Size]:: 
   M  : partition length
   x  : previous overlap-save vector [B x 1]
   xn : new input block [L x 1]
   dn : new desired block [L x 1]
   X  : previous matrix of F-domain input samples [B x P]
   W  : previous matrix of F-domain filter coef. [B x P]
   mu : adaptation constant
   n  : normalization flag, 0 means no normalization
   c  : if 0, unconstrained PBFDAF is used, 
        if -r (r +ve int), all partitions are constrained
        if +r (r +ve int), only r partitions are constrained
   b  : forgetting factor for input power estimation
   Px : previous estimate of the power of X [B x 1]
   ci : next partition to constrain 
Output parameters::
   W  : updated filter coefficients (F-domain)
   X  : updated matrix of past frequency input samples 
   x  : updated overlap-save input vector 
   y  : filter output block
   e  : error vector block
   Px : updated estimate of the power of X
   ci : next partition to constrain 
   w  : time domain filter (calculated only if required)

Example

iter = 5000;           % Number of samples to process
% Complex unknown impulse response
h    = [.9 + i*.4; 0.7+ i*.2; .5; .3+i*.1; .1];     
xt   = 2*(rand(iter,1)-0.5); % Input signal
% although xn is real, dn will be complex 
dt   = osfilter(h,xt);       % Unknown filter output 
en   = zeros(iter,1);        % estimation error
% Initialize RCPBFDAF with a filter of 2*4 coef.
P = 2; M = 4; L = M/2;
[W,x,d,e,y,Px,X,ci,w]=init_rcpbfdaf(L,M,P); 
%% Processing Loop
for (m=1:L:iter-L)
   xn = xt(m:m+L-1,:);             % input block
   dn = dt(m:m+L-1,:)+ 1e-3*rand;  % desired block    
   % call RCPBFDAF to calculate the filter output, 
   % estimation error and update the filter coef. 
   [W,X,x,y,e,Px,ci,w] = asptrcpbfdaf(M,x,xn,dn,...
                         X,W,0.06,1,1,0.98,Px,ci);

   % save the last error block to plot later
   en(m:m+L-1,:) = e;  
end;
% display the results
subplot(2,2,1);stem([real(w) imag((w))]); grid;
subplot(2,2,2);
eb = filter(.1, [1 -.9], en(1:m) .* conj(en(1:m)));
plot(10*log10(eb ));grid

Figure 4.14: The adaptive filter coefficients after convergence and the learning curve for the complex FIR system identification problem using the RCPBFDAF algorithm.

Running the above script will produce the graph shown in Fig. 4.14. The left side graph of the figure shows the adaptive filter coefficients after convergence which are almost identical to the unknown filter h. The right side graph shows the square error in dB versus time during the adaptation process, which is usually called the learning curve. The lower limit of the error signal power in the learning curve is defined here by the additive white noise added at the filter output (-60 dB).

Algorithm

asptrcpbfdaf() performs the following operations (see Fig. 4.2).

• buffers $L$ input samples, composes an overlap-save input vector $\vx (n)$ and computes its FFT, $\vX (f)$ and adds it to the input samples matrix $\mX (f)$
• element-wise multiplies $\mX (f)$ by the adaptive filter coefficients matrix $\mW (f)$ and evaluates the sum at each frequency bin (circular convolution in time domain for the $P$ partitions). The result is converted to time domain using IFFT and the linear convolution samples are extracted to produce the filter-output vector $\vy (n)$
• buffers $L$ desired samples and evaluates the current error block $\ve (n) = \vd (n) - \vy (n)$. The error vector is padded with zeros and transformed to frequency domain giving $\vE (f)$
• estimates the input signal power at each frequency bin and normalizes the step size at each bin.
• evaluates the cross-correlation between $\mX (f)$ and $\vE (f)$ to produce the block gradient vector. This vector is used to update the frequency domain filter coefficients.
• constrains only the partitioned given by the input argument $c$. This is performed by first taking the IFFT of the corresponding columns of $\mW (f)$, applying a rectangular window on the time domain coefficients, and taking the FFT of the windowed coefficients.

Remarks

• Supports both real and complex signals.
• asptrcpbfdaf() constrains the filter coefficients $\mW (f)$ rather than the gradient vector as in the official PBFDAF algorithms since this has been proven to result in a more stable update.
• The unconstrained RCPBFDAF (c = 0) saves $2P$ FFT operations each block on the cost of accuracy. Partially constrained filters $(c = r; r<P)$ save $2(P-r)$ FFT operations each block.
• The time domain filter coefficients $\vw (n)$ will be calculated and returned by asptrcpbfdaf() only if the output variable $w$ is given.
• The convergence properties of the constrained RCPBFDAF algorithm are superior to time domain algorithms since normalization is performed at each frequency bin which eliminates the eigen value spread problem. Unconstrained RCPBFDAF suffers from eigenvalue spread within each bin due to the overlap in each FFT block which can be reduced by decreasing the overlap.
• asptrcpbfdaf() introduces a processing delay between its input $x(n)$ and output $y(n)$ equals to the block length $L$, since the algorithm has to collect $L$ samples before processing a block. This delay is however $P$ times smaller than that introduced by BFDAF for the same total filter length.
• Very efficient for long adaptive filters since convolution and correlation are performed in frequency domain. Maximum efficiency is obtained when the block length is chosen to be equal to the filter length $L = M$.

Resources

The resources required for direct implement of the RCPBFDAF algorithm in real time is given in the table below. The computations given are those required to process L samples using the fully constrained RCPBFDAF. Partially constrained RCPBFDAF saves 2 FFT operations for each skipped partition, for instance calling asptrcpbfdaf with P=16 and c=1 saves 30 FFT operations every L samples, which can be huge saving for filters composed of many partitions. In the table below $C(FFT_B)$ is used to indicate the number of operations required to implement an FFT or IFFT of length $B=2^{nextpow2(M+L-1)}$

MEMORY	$B(G+P+4) + 3L + 4$
MULTIPLY	$6BP + (3+2P) C(FFT_B)$
ADD	$2BP + L + (3+2P) C(FFT_B)$
DIVIDE	$B + (3+2P) C(FFT_B)$

See Also

INIT_ RCPBFDAF, ECHO_ RCPBFDAF, ASPTPBFDAF, ASPTBFDAF.

Reference

[1] and [9] for detailed description of frequency domain adaptive filters.