apulSoft Blog

Feb 19th, 2022 - filters math dsp

All-Pass FIR Approximation


I decided to create all-pass filters to apply arbitrary phase changes using the frequency sampling method. This immediately raises red flags because by definition it is impossible to calculate the required no-gain impulse responses using inverse FFT (iFFT). This set me on a long journey and after many attempts, I found a relation between phase and magnitude change when doing an iFFT followed by windowing.

A better way to create all-pass approximating FIR filters would be a fitting procedure optimizing the IR coefficients to stay close to magnitude 1 while approximating the desired phase response. For me, that was not an option because performance is of the essence. I needed a solution that allows calculating changing FIR filters in real-time.

If complex frequency response bins with magnitude 1 and random phase values are transformed using the real iFFT the result is a random series of IR samples. The iFFT is circular - its beginning and end are connected. Such a signal is unusable as an FIR filter as it has no beginning and end. IR windowing does enforce this by multiplying each sample with a factor defined by a windowing function. These usually have 0 at the endings and 1 at the center.

To keep the math manageable I'm using one of the simplest signal windows, the Hann window:

$w[n] = 0.5 - 0.5 * cos((2pin)/N)$

A few years ago I read an eye-opening article on embedded.com about applying a window in the frequency domain.. allowing to change the order of windowing and FFT.

DSP Tricks: Frequency Domain Windowing

It boils down to this: The multiplication of the IR with a signal window is the same as applying a convolution to the frequency domain bins. If the window is a sum of cosines, the convolution kernel gets very small. The Hann window is just two cosines overlayed (one with freq 0..) and the convolution for each bin $b[n]$ is:

$b'[n] = 0.25*b[n-1] + 0.5*b[n] + 0.25*b[n+1]$

In the common case of doing linear-phase things with FFTs and iFFT, this just means applying the convolution to the bin magnitudes, as the phases are all 0 anyway. However, I realized the process also applies to complex bins. If we want to create an all-pass filter, all bins have magnitude 1, and only the phase angle changes. Looking at three such bins $a, b$ & $c$, the convolution can be visualized on the complex plane:

In polar form: ${(a = |a|*e^(ialpha),=,0.25*e^(ialpha)), (b = |b|*e^(ibeta),=,0.5*e^(ibeta)), (c = |c|*e^(igamma),=,0.25*e^(igamma)):}$

We are interested in the magnitude of the convolution result $x$ and how its size varies with the angles between the three source bins.

$x = a + b + c$

Looking at these complex numbers as vectors, the scalar product can be used to get the magnitude squared.

$|x|^2 = vecx*vecx = (veca+vecb+vecc)*(veca+vecb+vecc)$

$=veca*veca + vecb*vecb + vecc*vecc + 2*veca*vecb + 2*veca*vecc + 2*vecb*vecc$

$=|a|^2+|b|^2+|c|^2+2cos(alpha-beta)*|a|*|b|+2cos(beta-gamma)*|b|*|c|+2cos(alpha-gamma)*|a|*|c|$

using angles between the bins/vectors

${(/_vecavecb=alpha-beta=theta_0), (/_vecbvecc=beta-gamma=theta_1), (/_vecavecc=alpha-gamma=theta_0+theta_1):}$

$|x|^2=1/16+1/4+1/16+2cos theta_0*1/4*1/2+2cos theta_1*1/2*1/4+2cos(theta_0+theta_1)*1/4*1/4$

$|x|=sqrt((3+2cos theta_0 +2cos theta_1 +cos(theta_0+theta_1))/8)$

To figure out a limit for the phase angle, let's set $theta=theta_0=theta_1$. After simplification, the result is:

$|x|=1/2(cos theta+1)$

and solved for theta:

$theta = acos(2|x|-1)$

This means neighboring complex FFT filter bins with magnitude 1 that have a phase difference below theta won't change their magnitude to anything less than $|x|$ when the Hann window is applied after the iFFT.

The maximum possible gain change is not dependent on the bin center frequency or the FFT resolution.

Example: If all FFT bin magnitudes should stay between -0.1dB and 0dB, the phase cannot change more than $acos(2*10^(-0.1/20) - 1) = 0.21439...$ per bin.

Additionally, if any changes are larger than the limit, only the gains of the neighboring bins will be affected and it will be lowered. So if there are two bins in a row where we don't care about gain, an arbitrary phase jump can be executed.

All gain changes caused by the phase changes will lower gain. The expected change can be approximated and applied in inverse before doing the windowing to end up with gains staying closer to 0dB. With $theta$ being the average phase change on both sides we get for the input gain $g$ that results in gain 1 (0dB) after the windowing:

$1 = g*1/2(cos(theta) + 1)$

$g = 2/(cos(theta) + 1)$

as $theta$ is an approximation and quite small, a short taylor approximation can be used:

$g ~~ 1 + 0.25*theta^2$

Setting all gain values to their respective $g$ values before applying the window halves the gain error.

This knowledge can be used to smooth arbitrary phase changes into a complex frequency response approximating an all-pass fir filter.

Such an approach is the basis for apulSoft apQualizr2's mixed phase mode - IIR filter phase changes are compensated to get linear phase filtering where phase changes are small.

 
Comments