High Dynamic Range FIR Filters

Question

Related to this question on using the Windowing method for FIR filter design versus the optimized algorithms such as least squares (firls in MATLAB, Octave and Python scipy.signal, my usual "go-to" for filter design), I ran into another interesting case where, as far as I could tell, could not be resolved with the least squares or Parks-McClellan algorithms. Specifically, I am in need of a set of high dynamic range filters, with stopband rejection (and passband distortion) less than -212 dB (and in general as a personal challenge, I would like to know how to approach the noise floor of any available numerical precision if ever needed).

Note since I am being asked: This is not for an application of processing analog signals captured with an ADC or for signals that will be converted to analog with a DAC. This is for an dsp algorithm where a very small quantization noise is needed or equivalently very high SNR.

My question is as follows: The least squares algorithm "converged" and provided results, but clearly in comparison to what I achieved with a Kaiser window with the same number of taps, the least squares solution is not "optimum in the least squares sense" as advertised. This was a surprise to me since there wasn't any indication otherwise of a convergence issue (as did occur with Parks-McClellan, here with least squares I got a result with no complaint), so then even in lower SNR cases, under which conditions can we no longer "trust" that the solution provided by the algorithm is actually optimum in the least squares sense?

Is what I am finding a theoretical limit/ constraint as part of the least squares algorithm and if so what is that constraint; or am I approaching it incorrectly somehow? Can the least squares algorithm achieve rejections below 212 dB for my specific example case given below? Are there other implementations of the least squares algorithm that can do this, and the specific ones in Python and Octave are somehow inferior?

Example

With the following in either Octave or Python:

coeff = firls(287, [0, .4, .6, 1], [1, 1, 0, 0])

I get the following results for this half-band filter shown below, which is directly compared to Kaiser windowing a Sinc impulse response (window method of FIR design) with the same number of coefficients. Note that increasing the number of coefficients further with the least squares algorithm was not productive in reducing the stopband further. I was also unsuccessful with using the weighting functions, but in this case I do want equal performance in passband and stopband: the plots show the passband magnitude error as well in orange as $20\log_{10}(1-|H(\omega)|)$

I do notice the wasted "chatter" in the wider transition band for the least squares case, which I can eliminate by tightening that transition band (such as [0, .49, .51, 1]), but doing that also reduces the stopband rejection (as expected). Increasing the number of taps from there eventually just adds the "chatter" back in to that "don't care region" rather than providing more rejection.

The least squares algorithm is detailed here, but I do not see any mention as to constraints on the ultimate dynamic range that can be achieved, nor what would be limiting it. It does concern me that the coefficients were returned without any issue, yet they are not actually optimum in the least squares sense, a point which up to now I assumed was assured by the algorithm, as long as it successfully converged- so at which point does that assurance no longer hold?

Answer 1

Like @MattL. and @aconcernedcitizen say, the issue is numerical.

Python's scipy.signal.firls uses internally the solver scipy.linalg.solve. For your input, the solver throws a "matrix singular" error, but firls suppresses the error and falls back to another solver scipy.linalg.lstsq which doesn't throw an error but also doesn't get the problem solved properly.

Increasing the numerical resolution helps.

Arbitrary precision firls implementation

Here is my modified version of scipy.signal.firls. It uses mpmath for an arbitrary precision floating point data type and for arbitrary precision solvers. I removed all comments for brevity. The code is licensed under the BSD 3-Clause license of the original Scipy 1.7 source and is copyright 2001-2002 Enthought, Inc. and 2003-2019, SciPy Developers; all rights reserved.

import numpy as np
from scipy.linalg import (hankel, toeplitz)
import mpmath as mp

def _get_fs(fs, nyq):
    if nyq is None and fs is None:
        fs = 2
    elif nyq is not None:
        if fs is not None:
            raise ValueError("Values cannot be given for both 'nyq' and 'fs'.")
        fs = 2*nyq
    return fs

def firls(numtaps, bands, desired, weight=None, nyq=None, fs=None):
    nyq = 0.5 * _get_fs(fs, nyq)

    numtaps = int(numtaps)
    if numtaps % 2 == 0 or numtaps < 1:
        raise ValueError("numtaps must be odd and >= 1")
    M = (numtaps-1) // 2

    nyq = mp.mpf(nyq)
    if nyq <= 0:
        raise ValueError('nyq must be positive, got %s <= 0.' % nyq)
    bands = np.asarray(bands, dtype=mp.mpf).flatten() / nyq
    if len(bands) % 2 != 0:
        raise ValueError("bands must contain frequency pairs.")
    if (bands < 0).any() or (bands > 1).any():
        raise ValueError("bands must be between 0 and 1 relative to Nyquist")
    bands.shape = (-1, 2)

    desired = np.asarray(desired, dtype=mp.mpf).flatten()
    if bands.size != desired.size:
        raise ValueError("desired must have one entry per frequency, got %s "
                         "gains for %s frequencies."
                         % (desired.size, bands.size))
    desired.shape = (-1, 2)
    if (np.diff(bands) <= 0).any() or (np.diff(bands[:, 0]) < 0).any():
        raise ValueError("bands must be monotonically nondecreasing and have "
                         "width > 0.")
    if (bands[:-1, 1] > bands[1:, 0]).any():
        raise ValueError("bands must not overlap.")
    if (desired < 0).any():
        raise ValueError("desired must be non-negative.")
    if weight is None:
        weight = np.ones(len(desired), dtype=mp.mpf)
    weight = np.asarray(weight, dtype=mp.mpf).flatten()
    if len(weight) != len(desired):
        raise ValueError("weight must be the same size as the number of "
                         "band pairs (%s)." % (len(bands),))
    if (weight < 0).any():
        raise ValueError("weight must be non-negative.")

    n = np.arange(numtaps)[:, np.newaxis, np.newaxis]
    q = np.dot(np.diff(np.vectorize(mp.sincpi)(bands * n) * bands, axis=2)[:, :, 0], weight)

    Q1 = toeplitz(q[:M+1])
    Q2 = hankel(q[:M+1], q[M:])
    Q = Q1 + Q2

    n = n[:M + 1]
    m = (np.diff(desired, axis=1) / np.diff(bands, axis=1))
    c = desired[:, [0]] - bands[:, [0]] * m
    b = bands * (m*bands + c) * np.vectorize(mp.sincpi)(bands * n)

    b[0] -= m * bands * bands / 2.
    b[1:] += m * np.vectorize(mp.cospi)(n[1:] * bands) / (mp.pi * n[1:]) ** 2
    b = np.dot(np.diff(b, axis=2)[:, :, 0], weight)

    a = mp.lu_solve(Q, b)

    coeffs = np.hstack((a[:0:-1], 2 * a[0], a[1:]))
    return coeffs

With enough precision (set through variable mp.mp.prec) this version of firls works fine on your inputs. To properly show the, in this case, sub -400 dB magnitude values, the magnitude frequency response calculation must be done in arbitrary precision, which is quite slow without an arbitrary precision Fast Fourier Transform (FFT) implementation.

import mpmath as mp
import matplotlib.pyplot as plt

mp.mp.prec = 150

N = 287
coeff = firls(N, [0, .4, .6, 1], [1, 1, 0, 0])

num_freqs = 256 + 1
f = np.arange(num_freqs, dtype=mp.mpf)/(num_freqs - 1)
phases = np.arange(N)*f.reshape((num_freqs, 1))
phasors = np.vectorize(mp.cospi)(phases) - 1j*np.vectorize(mp.sinpi)(phases)
response = np.dot(phasors, coeff)
plt.plot(0.5*f, 20*np.vectorize(mp.log10)(abs(response)))
plt.grid(True)
plt.ylabel("dB")
plt.xlabel("Frequency (cycles/sample)")
plt.show()

Output:

So yeah, with enough precision, least squares filter design works fine!

Rounding the arbitrary precision coefficients to double precision floating point by:

coeff = coeff.astype(np.float64)

doesn't give as nice results:

Arbitrary precision Kaiser windowed FIR filter design

Here is also an arbitrary precision implementation (not extensively tested) of Kaiser windowed lowpass finite impulse response (FIR) filter design. In the test, I adjusted the Kaiser parameter beta so that the first frequency response zero matches that of the least square filter with your input parameters:

def mpf_kaiser_lowpass_fir(N, cutoff, beta):
  n = np.arange(N, dtype=mp.mpf)
  alpha = mp.mpf((N-1)/2)
  win = np.vectorize(lambda x: mp.besseli(0, x))(beta*np.vectorize(mp.sqrt)(1 - ((n-alpha)/alpha)**2))/mp.besseli(0, beta)
  coeffs = cutoff*win*np.vectorize(mp.sincpi)((n-(N-1)/2)*cutoff)
  return coeffs

coeff2 = mpf_kaiser_lowpass_fir(N, 0.5, 44.9535)

It looks like the least squares design (blue) is winning against the Kaiser windowed design (orange):

Zoomed in:

Adjusting beta to equate the two frequency responses at 0.3 cycles / sample gives a qualitatively similar comparison result (not shown).

Answer 2

The problem lies in the formulation of the desired response, and especially in the "don't care" region, which is extremely wide for the chosen filter length. Even though I can't give any exact relation between transition band width and filter length, I know that in the case of a least squares design, the matrix of the system of linear equations that must be solved becomes ill-conditioned if the "don't care" region is too wide for a given filter length.

A similar effect is known for the design of equiripple (Parks-McClellan) FIR filters, where the frequency response exhibits "bumps" in the transition band(s) if the filter length is too large for the chosen transition band width. This phenomenon occurs for bandpass and multi-band filters.

In the case of your specific example, using Octave's firls function, I get a warning that the matrix is ill-conditioned. I'm almost sure that you should get the same warning. When I try to use the same specifications to design a Parks-McClellan filter (using remez.m), the algorithm stops with the error message "too many extremals - cannot continue".

So indeed both mainstream algorithms fail with that specification. However, both algorithms can be used to achieve extremely high stopband attenuations. But in that case the specification needs to be changed. There are two approaches: 1. choose (by trial and error) a transition band that "matches" the given filter length, such that there are no more artefacts of the frequency response in the transition band. 2. instead of a "don't care" region without any specification, choose a desired response in the transition band.

With the first approach we will end up with a very high filter order and a very narrow transition band, which is probably much narrower than needed. The difficulty of the second approach lies in the choice of the desired response in the transition band. If we choose some response that is not "natural" for an FIR filter of the given length, we will end up with a large approximation error. I remember a paper discussing that problem:

Least squared error FIR filter design with transition bands by Burrus, Soewito, and Gopinath.

Just in order to show that a least squares method can achieve an extremely small error, assuming that we somehow know how to specify an ideal response in the transition band, I designed a linear-phase least squares filter with the transition band specification taken from the filter designed by the Kaiser window. Of course, we might as well just use the Kaiser windowed filter in that case, but the experiment was just to show that there is no reason why a least squares filter can't achieve a very small approximation error, given that we know how to choose an appropriate filter spec.

The design cannot be done with firls because that function allows only a piecewise linear desired response. I used the function lslevin, which allows arbitrary complex desired responses specified on a frequency grid.

The figure below shows the result. The responses of both filters (Kaiser window and least squares) are virtually identical; the maximum difference between their magnitudes is in the order of $10^{-12}$.

In sum, a least squares design can achieve extremely small errors. However, for frequency-selective filters, the transition band needs to be relatively wide for achieving such small errors, which causes numerical ill-conditioning of the corresponding system of linear equations. In order to avoid that, one needs to specify a desired response in the transition band. The remaining problem is how to optimally choose such a response. The above-mentioned paper could be a good starting point. However, I believe that that problem is mainly of academic interest, because if such a filter is really needed, the (Kaiser) window method is very likely the best choice.