7
$\begingroup$

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

firls vs least squares

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?

$\endgroup$
9
  • 1
    $\begingroup$ (+1) interesting question. Just curious, what sort of applications would need a stopband rejection of less than -212 dB? $\endgroup$
    – mhdadk
    Jan 5 at 9:17
  • 2
    $\begingroup$ The algorithm seems to encounter numeric issues (Octave complains about warning: matrix singular to machine precision, rcond = 2.66832e-17), and it's probably because of the very wide transition width relative to the order. If you try the frequency vector [0 .49 .51 1] you'll get a much better behaved response. Kaiser is better behaved because it's a simple multiplication (e.g. both the sinc() and the window are static computations, not the result of some algorithm). $\endgroup$ Jan 5 at 9:42
  • 2
    $\begingroup$ @mhdadk I work in the design of high precision atomic clocks. The requirement is driven from quantization noise effects in proprietary algorithms to ensure the noise added at critical nodes is less than the fractional frequency we are trying to measure. $\endgroup$ Jan 5 at 13:19
  • 3
    $\begingroup$ @Max I am not filtering signals captured from an ADC but working in applications where I need to maintain a very high SNR effectively. The analogy may be similar of how we require a 35 bit accumulator in an NCO when we want a precision of $f_clk/2^{35}$ -- we truncate later after the accumulation but within the accumulator node we require a very high number of bits. $\endgroup$ Jan 5 at 13:27
  • 2
    $\begingroup$ @Max maybe a similar question is why do we need the simultaneous dynamic range (as I define here: dsp.stackexchange.com/questions/38832/…) of double precision floating point? We do run into applications of numerical precision deeper within signal processing where we actually need greater than the 150 dB effective SNR that single precision floating point provides. But this has nothing to do with how many bits we grab from an ADC or send to a DAC. $\endgroup$ Jan 5 at 13:44

2 Answers 2

4
$\begingroup$

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:

enter image description here

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:

enter image description here

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

enter image description here

Zoomed in:

enter image description here

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

$\endgroup$
10
  • 1
    $\begingroup$ WOW, this is really impressive and interesting. Thanks Olli. Good to know if I ever hit a double precision floating point wall. Happy New Year to you and good to see you. $\endgroup$ Jan 6 at 20:04
  • 1
    $\begingroup$ Happy New Year! :-) I haven't had much time lately to relax in the DSP zone... $\endgroup$ Jan 6 at 20:42
  • 1
    $\begingroup$ Your result after rounding doesn't give as nice results but it is nicely consistent with the expected dynamic range of double precision floating point (as we would expect) $\endgroup$ Jan 6 at 21:54
  • 1
    $\begingroup$ That's interesting, thanks for adding! I wonder why the stopband ripple envelope of the Kaiser filter is non-monotonic. Could it be that is would need even more precision? $\endgroup$
    – Matt L.
    Jan 7 at 14:37
  • 1
    $\begingroup$ What surprised me in the designs shown in my answer, is that both filters are virtually identical. Since the computations are totally different in both cases I wouldn't expect this to be caused purely by numerical effects. $\endgroup$
    – Matt L.
    Jan 7 at 14:42
6
$\begingroup$

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}$.

enter image description here

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.

$\endgroup$
6
  • $\begingroup$ Very very nice Matt. I will try that again as I don’t recall seeing that error but if it is indeed there, that would make me feel a lot better. I guess second to that I can assume an ill-conditioned solution if I see chatter in the transition band so that is one indication I guess. I tried with a tighter transition to eliminate that but couldn’t get the lower rejection— I need to look closer at what you did and try that again. Thanks so much this addressed my primary question really well. $\endgroup$ Jan 5 at 14:37
  • $\begingroup$ @DanBoschen: I'm pretty sure that transition band artefacts and numerical ill-conditioning go hand in hand with LS designs, so seeing one you can be sure of the other. As I said in my answer, for such extreme designs, windowing is probably the best solution. $\endgroup$
    – Matt L.
    Jan 5 at 14:40
  • $\begingroup$ Thanks - yes I have no issue with doing the windowing in this case and more just wanted to be sure I understood the limits of least squares. Also I like the clear example where windowing is more straightforward and as you showed just as effective. But I also really like that you showed we can indeed still get there with firls (just with more aggravation) $\endgroup$ Jan 5 at 14:43
  • $\begingroup$ If the whole problem is that toeplitz matrix going ill conditioned thanks to an overly broad filter spec, it seems like the easiest solution would be to apply some slight regularization to push it towards a "natural" solution. For general filter design that would be a pain, but for a specific spec, you could probably just fiddle with it until a reasonable design fell out. $\endgroup$
    – Steve Cox
    Jan 5 at 19:27
  • $\begingroup$ @SteveCox yes indeed and as far as the algorithm itself I would prefer that it do that “under the hood” rather than make the user do the fiddling since no trade space was actually involved. Makes me think I could put a wrapper around the firls function to create “thefirls” $\endgroup$ Jan 6 at 1:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.