3
$\begingroup$

I need to design a short, symmetric low-pass FIR filter with a monotonic frequency response (no ripple). In other words, the differential of its amplitude response will be zero or less for all frequencies. How can I compute N coefficients such that I minimise stopband energy?

$\endgroup$
3
  • 1
    $\begingroup$ It seems to me that your filter is not fully specified. What I mean is that even a filter as simple as h=[1,1] would be an option: it is symmetric (has linear phase), has no ripple, and a zero at Nyquist. If that filter is not an option then your question is incomplete, and if it is in fact an option, then I wonder why you didn't find it yourself. $\endgroup$
    – Matt L.
    Feb 2, 2023 at 11:00
  • $\begingroup$ Yes, you're right, thank you, I haven't given enough information. Will edit. $\endgroup$
    – paperjam
    Feb 2, 2023 at 11:05
  • 1
    $\begingroup$ This can be formulated as a quadratic programming problem: it has a quadratic objective function (stopband energy) to be minimized, subject to linear constraints on the filter coefficients (derivative of amplitude $\le 0$). $\endgroup$
    – Matt L.
    Feb 2, 2023 at 13:09

2 Answers 2

8
$\begingroup$

This answer shows how to get coefficients for zero-phase "maxflat" filters. They do not minimize stop-band energy (as noted by @MattL. in a comment) generally, but may still be useful.

With $N$ unique coefficients $b_i$, $0 \le i < N$, the length of a zero-phase filter will be an odd number $2N - 1$. This is the delay-compensated version of the filter, with impulse response $[b_{N-1}, \ldots, b_1, b_0, b_1, \ldots b_{N-1}]$ symmetric around time zero. Its frequency response is:

$$H(\omega) = \sum_{i\ =\ -(N - 1)}^{N - 1} 2b_{|i|}e^{j i \omega} = b_0 + \sum_{i = 1}^{N - 1} 2b_i\cos(i \omega),\tag{1}\label{15347-1}$$

where $j$ is the imaginary unit and frequency being expressed such that $\omega = \pi$ represents half the sampling frequency.

To create a low-pass filter with a nonincreasing frequency response we can first constrain $H(0) = 1$ and $H(\pi) = 0$. For $N = 2$ that gives as many equations as unknowns ($b_0$ and $b_1$) and the filter coefficients can be solved unambiguously. For larger $N$ we can constrain long enough runs of early derivatives of $H(\omega)$ to zero at $\omega = 0$ and $\omega = \pi$ to get as many equations as unknowns. At which of the two frequencies you constrain to zero more derivatives will result in flattening more the pass or the stop band, allowing to tune the effective cutoff frequency in steps. For the given number of the coefficients the pass and stop band frequency responses are maximally flat for these filters in the sense that one cannot be made flatter without reducing flatness of the other.

The $n$th derivative of $H(\omega)$, denoted by $H^{(n)}(\omega)$, can be calculated by:

$$H^{(n)}(\omega) = \sum_{i = 1}^{N - 1}[2 \cos(i\omega), -2 \sin(i\omega), -2 \cos(i\omega), 2 \sin(i\omega)]_n\ i^n b_i,\quad n>0,\tag{2}\label{15347-2}$$

where $[\text{elements of a vector}]_n$ is a shorthand notation to get an element of a vector at zero-based index $n$ modulo the length of the vector. For example, $[0, 1, 2, 3]_4 = 0$.

To ensure maximum flatness, we are only interested in $H(0)$ and $H(\pi)$ and the derivatives at $\omega = 0$ and at $\omega = \pi$:

$$\begin{eqnarray}&&H(0) = b_0 + \sum_{i = 1}^{N - 1} 2 b_i\\ &&H(\pi) = b_0 + \sum_{i = 1}^{N - 1}[2, -2]_i b_i\\ &&H^{(n)}(0) = \sum_{i = 1}^{N - 1}[2, 0, -2, 0]_n i^n b_i,\quad n>0\\ &&H^{(n)}(\pi) = \sum_{i = 1}^{N - 1}[2, 0, -2, 0]_n [1, -1]_i i^n b_i,\quad n>0\end{eqnarray}\tag{3}\label{15347-3}$$

Note that at those frequencies the derivatives with odd $n$ are zero. Therefore we will only ever constrain even derivatives. Solving the resulting system will give rational number coefficients. Those can be calculated using the Python script at the end of the answer.

For a simple example, let's design a filter with $N = 3$ having coefficients $[b_2, b_1, b_0, b_1, b_2]$. That's 3 unknowns allowing us to have 3 equations. Two equations will be wasted by $H(0) = 1$ and $H(\pi) = 0$. That leaves one equation to constrain either $H^{(2)}(0) = 0$ or $H^{(2)}(\pi) = 0$. We choose the latter to further flatten the stop band. The system to solve will be:

$$\cases{b_0 + 2b_1 + 2b_2 = 1\\ b_0 - 2b_1 + 2b_2 = 0\\ 2b_1-8b_2 = 0,}\tag{4}\label{15347-4}$$

with a single solution $b_0 = \frac{3}{8} = 0.375$, $b_1 = \frac{1}{4} = 0.25$, $b_2 = \frac{1}{16} = 0.0625$. The frequency response is shown in Fig. 1 below.

A non-increasing low-pass frequency response over frequency 0 to pi
Figure 1. frequency response of the filter with coefficients [0.0625 0.25 0.375 0.25 0.0625]

If you just want to minimize stop band energy and don't care about the pass band response, then you can obtain the maximally flat stop band frequency response solution by starting from the unique $N=2$ solution impulse response [0.25 0.5 0.25] and by convolving the impulse response with [0.25 0.5 0.25] as many times as you please. Each convolution corresponds to frequency-domain multiplication, increasing the number of frequency response zeros at $\omega = \pi$ by two.

With $N \le 3$, a maxflat stopband minimizes stopband energy, but for $N > 3$ that is not generally true. Let's look at the $N = 3$ case. We use the weakest possible pass band specification $H(0) = 1$ to normalize the filter. We see from Eq. $\eqref{15347-3}$ that the frequency response is always flat at $\omega = 0$ and at $\omega = \pi$. We argue that $H(\pi) = 0$ because that is the only case for which we cannot bring the frequency response closer to zero everywhere in $0 < \omega \le \pi$ by "stretching" the response vertically by adjusting the coefficients in the response defined by Eq. $\eqref{15347-1}$. The operation does not affect the monotonicity or non-negativity of the frequency response and would reduce the pass band error for any cutoff frequency and any reasonable definition of pass band error, including least squares (that is, minimum stopband energy). It must then be that for the optimal filter, $H(\pi) = 0$.

That's already two constraints in total, so we have only one more degree of freedom left to adjust. For easier analysis, we can convert Eq. $\eqref{15347-1}$, which is a sum of cosine terms, into a polynomial by a monotonicity preserving change of variables $\omega = \operatorname{acos}(x)$, giving a mapping from $\omega \in [0\ldots\pi]$ to $x \in [1\ldots-1]$ (note the change in direction), and by using the relation $T_i\left(\cos\left(\omega\right)\right) = \cos\left(i\omega\right)$, where $T_i$ is a Chebyshev polynomial of the first kind. Each term for $i>0$ is rewritten using:

$$2b_i\cos(i\omega) = 2b_iT_i\left(\cos\left(\omega\right)\right) = 2b_iT_i\left(\cos\left(\operatorname{acos}\left(x\right)\right)\right) = 2b_iT_i(x)\tag{5}\label{15347-5}$$

This way we get a mapped zero-phase frequency response $G(x) = H(\omega)$ as function of $x$:

$$G(x) = b_0 + \sum_{i = 1}^{N - 1} 2b_iT_i(x),\tag{6}\label{15347-6}$$

For $N = 3$ it is simply:

$$G(x) = b_0 + 2b_1x + 2b_2\,(2x^2-1) = b_0 - 2b_2 + 2b_1x + 4b_2x^2\tag{7}\label{15347-7}$$

We reintroduce the earlier constraints as $G(1) = 1$ and $G(-1) = 0$. This resolves the two first coefficients as $b_0 = \frac{1}{2} - 2b_2$ and $b_1 = \frac{1}{4}$, giving a constrained mapped frequency response:

$$G(x) = \frac{1}{2} + \frac{1}{2}x + (4x^2 - 4)b_2\tag{8}\label{15347-8}$$

The factor $4x^2 - 4$ is zero at $x = -1$ and $x = 1$, and negative in $-1 < x < 1$. The only constrained monotone non-negative frequency response that cannot be brought closer to zero everywhere in $0 < \omega < \pi$ by increasing $b_2$ is the one for which $G(x) = 0$ at $x$ given by:

$$\begin{eqnarray}&G^{(1)}(x) = 0\\ \Rightarrow\quad&8b_2x + \frac{1}{2} = 0\\ \Rightarrow\quad&x = - \frac{1}{16b_2}\end{eqnarray}\tag{9}\label{15347-9}$$.

Equating $G(- \frac{1}{16b_2}) = 0$ gives:

$$\begin{eqnarray}&- 4b_2 - \frac{1}{64b_2} + \frac{1}{2} = 0\\ \Rightarrow\quad&b_2 = \frac{1}{16}\end{eqnarray}\tag{10}\label{15347-10}$$

That's the $N=3$ maxflat stopband filter.

To help analyze maxflat stopband filters with general $N$, I asked a question on Mathematics Stack Exchange, and an answer effectively showed that for $N > 3$ there are filters with a monotone non-negative response that, compared to the maxflat stopband filter, can give a lower response at some frequency even when both filters are normalized by $H(0) = 1$. The example given in the answer corresponds at $N=4$ to filter $\left[\frac{3}{64}, \frac{3}{32}, \frac{13}{64}, \frac{5}{16}, \frac{13}{64}. \frac{3}{32}, \frac{3}{64}\right]$ by substitution $x \rightarrow \frac{1}{2}x + \frac{1}{2}$, conversion to Chebyshev basis and replacing the terms with harmonic cosines to obtain the frequency response. For a stopband spanning the whole frequency range $0 \le \omega \le \pi$, that would give about 10 % lower stop band energy than the $N=4$ maxflat stopband filter $\left[\frac{1}{64}, \frac{3}{32}, \frac{15}{64}, \frac{5}{16}, \frac{15}{64}, \frac{3}{32}, \frac{1}{64}\right]$. The frequency responses are shown in Fig. 2.

enter image description here
Figure 2. Frequency responses of the filter $\left[\frac{3}{64}, \frac{3}{32}, \frac{13}{64}, \frac{5}{16}, \frac{13}{64}. \frac{3}{32}, \frac{3}{64}\right]$ (blue) and the $N=4$ maxflat stopband filter (red).

There are also Optimum "L" filters. A FIR filter with an Optimum "L" stopband would have the zeros of the derivative of the frequency response distributed over the stop band in a way that maximizes the slope at the cutoff frequency. Such filters may be better at reducing stopband energy compared to maxflat stopband filters.

Here's a literature reference about maxflat FIR filters that looks interesting: Chung, Daewon, Woon Cho, Inyeob Jeong, and Joonhyeon Jeon. 2021. "Design of Cut Off-Frequency Fixing Filters by Error Compensation of MAXFLAT FIR Filters" Electronics 10, no. 5: 553. https://doi.org/10.3390/electronics10050553

Python maxflat FIR low-pass filter coefficient calculation

This Python script calculates maxflat FIR low-pass filter coefficients for the given passband flatness P and stopband flatness S. The symbolic math package sympy is used to get exact rational number coefficients.

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
import mpmath as mp
import sympy

def maxflat(P, S):
  # P = Pass band flatness, positive integer
  # S = Stop band flatness, positive integer
  assert P > 0
  assert S > 0
  N = P + S # Number of unique coefficients. Filter will have 2N - 1 coefficients
  table4 = [2, 0, -2, 0]
  table2 = [1, -1]

  def H_0(n, i):
    if i == 0:
      if n == 0:
        return 1
      else:
        return 0
    else:
      return table4[n%4]*sympy.Pow(i, n)

  def H_pi(n, i):
    if i == 0:
      if n == 0:
        return 1
      else:
        return 0
    else:
      return table4[n%4]*table2[i%2]*sympy.Pow(i, n)

  a = sympy.zeros(N, N)
  for n in range(P):
    for i in range(N):
      a[n, i] = H_0(n*2, i)
  for n in range(S):
    for i in range(N):
      a[P + n, i] = H_pi(n*2, i)
  b = sympy.eye(N, 1)
  x = a.LUsolve(b)
  return list(sympy.Matrix.vstack(sympy.Matrix(x[-1:0:-1]), x))

P = 15
S = 35
N = P + S
b = maxflat(P, S)
b_float64 = list(np.array(b, dtype=np.float64))
print("P =", P)
print("S =", S)
print("N =", N)
print("Filter length =", 2*N - 1)
print("Exact:")
print(b)
print("float64:")
print(repr(b_float64))
w, h = signal.freqz(b_float64)
fig, ax1 = plt.subplots()
ax1.plot(w, np.real(h*np.exp(1j*w*(N-1))))
ax1.set_ylabel('Frequency response')
ax1.set_xlabel('$\omega$')
plt.show()

Example output

P = 15
S = 35
N = 50
Filter length = 99
Exact:
[60290077905/39614081257132168796771975168, 4923689695575/79228162514264337593543950336, 47204565464385/39614081257132168796771975168, 70002023063175/4951760157141521099596496896, 567668597301827/4951760157141521099596496896, 6508133039611035/9903520314283042199192993792, 12891158559189855/4951760157141521099596496896, 30126369125042085/4951760157141521099596496896, -2348311768807095/9903520314283042199192993792, -1300762204959251853/19807040628566084398385987584, -2646011521898726415/9903520314283042199192993792, -2167968001596345585/4951760157141521099596496896, 2804604525116597385/4951760157141521099596496896, 48368020984829412965/9903520314283042199192993792, 52858619408481968277/4951760157141521099596496896, -1970460616169115/4951760157141521099596496896, -1177325915813446793415/19807040628566084398385987584, -5713411640861450185065/39614081257132168796771975168, -843619214982331712295/19807040628566084398385987584, 2788813008762433384473/4951760157141521099596496896, 6671757799662009867225/4951760157141521099596496896, 3842694615105320832585/9903520314283042199192993792, -21783850434977250258275/4951760157141521099596496896, -46269202784045093873685/4951760157141521099596496896, -11398366009445266995507/9903520314283042199192993792, 563505082827768663706215/19807040628566084398385987584, 487574275710372412241925/9903520314283042199192993792, -37931076518167281266055/4951760157141521099596496896, -743351781673767967432725/4951760157141521099596496896, -1948852258912371721744521/9903520314283042199192993792, 549449500869281093733855/4951760157141521099596496896, 3158225231738691969966755/4951760157141521099596496896, 23000491348952162430469455/39614081257132168796771975168, -56175086812603900278484935/79228162514264337593543950336, -85848194261877596737454721/39614081257132168796771975168, -2844116251858741058375085/2475880078570760549798248448, 7602962723657608571703495/2475880078570760549798248448, 29108338515501648967444335/4951760157141521099596496896, 2052663369503261347831635/2475880078570760549798248448, -25001230805364534689842839/2475880078570760549798248448, -63084723033433146475033415/4951760157141521099596496896, 41408244257891574408551955/9903520314283042199192993792, 134542668383383518249441345/4951760157141521099596496896, 54710983442922686326925835/2475880078570760549798248448, -60180577167462385720128651/2475880078570760549798248448, -334341555218410707021236175/4951760157141521099596496896, -76076583929405887081624335/2475880078570760549798248448, 279904412570455622281448025/2475880078570760549798248448, 2855025008218647347270769855/9903520314283042199192993792, 7281563888583524448979432097/19807040628566084398385987584, 2855025008218647347270769855/9903520314283042199192993792, 279904412570455622281448025/2475880078570760549798248448, -76076583929405887081624335/2475880078570760549798248448, -334341555218410707021236175/4951760157141521099596496896, -60180577167462385720128651/2475880078570760549798248448, 54710983442922686326925835/2475880078570760549798248448, 134542668383383518249441345/4951760157141521099596496896, 41408244257891574408551955/9903520314283042199192993792, -63084723033433146475033415/4951760157141521099596496896, -25001230805364534689842839/2475880078570760549798248448, 2052663369503261347831635/2475880078570760549798248448, 29108338515501648967444335/4951760157141521099596496896, 7602962723657608571703495/2475880078570760549798248448, -2844116251858741058375085/2475880078570760549798248448, -85848194261877596737454721/39614081257132168796771975168, -56175086812603900278484935/79228162514264337593543950336, 23000491348952162430469455/39614081257132168796771975168, 3158225231738691969966755/4951760157141521099596496896, 549449500869281093733855/4951760157141521099596496896, -1948852258912371721744521/9903520314283042199192993792, -743351781673767967432725/4951760157141521099596496896, -37931076518167281266055/4951760157141521099596496896, 487574275710372412241925/9903520314283042199192993792, 563505082827768663706215/19807040628566084398385987584, -11398366009445266995507/9903520314283042199192993792, -46269202784045093873685/4951760157141521099596496896, -21783850434977250258275/4951760157141521099596496896, 3842694615105320832585/9903520314283042199192993792, 6671757799662009867225/4951760157141521099596496896, 2788813008762433384473/4951760157141521099596496896, -843619214982331712295/19807040628566084398385987584, -5713411640861450185065/39614081257132168796771975168, -1177325915813446793415/19807040628566084398385987584, -1970460616169115/4951760157141521099596496896, 52858619408481968277/4951760157141521099596496896, 48368020984829412965/9903520314283042199192993792, 2804604525116597385/4951760157141521099596496896, -2167968001596345585/4951760157141521099596496896, -2646011521898726415/9903520314283042199192993792, -1300762204959251853/19807040628566084398385987584, -2348311768807095/9903520314283042199192993792, 30126369125042085/4951760157141521099596496896, 12891158559189855/4951760157141521099596496896, 6508133039611035/9903520314283042199192993792, 567668597301827/4951760157141521099596496896, 70002023063175/4951760157141521099596496896, 47204565464385/39614081257132168796771975168, 4923689695575/79228162514264337593543950336, 60290077905/39614081257132168796771975168]
float64:
[1.521935533823476e-18, 6.21457009644586e-17, 1.1916107597695764e-15, 1.4136795975915104e-14, 1.1463976026446368e-13, 6.57153500278571e-13, 2.6033487386496266e-12, 6.083971793664778e-12, -2.37118892503337e-13, -6.567170883081183e-11, -2.671788856819528e-10, -4.378176512587472e-10, 5.663853733044288e-10, 4.883922024683702e-09, 1.0674713178958855e-08, -3.979313524156213e-13, -5.943976881208015e-08, -1.4422678652512737e-07, -4.25918859259393e-08, 5.631963019736235e-07, 1.34735075769772e-06, 3.8801299872766606e-07, -4.399213561173833e-06, -9.343991089171552e-06, -1.1509408420161798e-06, 2.844973630311391e-05, 4.9232420415918535e-05, -7.660119899681e-06, -0.00015011869680353805, -0.0001967837897097763, 0.0001109604430410174, 0.0006377985062915135, 0.0005806140296339985, -0.0007090292773417542, -0.0021671130955844492, -0.0011487294059494798, 0.003070812188951629, 0.005878382149329478, 0.0008290641324955418, -0.010097916705156768, -0.012739858359749346, 0.004181164166258318, 0.02717067549997218, 0.02209759023324968, -0.024306741545496235, -0.06751973936706467, -0.030727087546712783, 0.1130524919171508, 0.28828385438873455, 0.3676250291566483, 0.28828385438873455, 0.1130524919171508, -0.030727087546712783, -0.06751973936706467, -0.024306741545496235, 0.02209759023324968, 0.02717067549997218, 0.004181164166258318, -0.012739858359749346, -0.010097916705156768, 0.0008290641324955418, 0.005878382149329478, 0.003070812188951629, -0.0011487294059494798, -0.0021671130955844492, -0.0007090292773417542, 0.0005806140296339985, 0.0006377985062915135, 0.0001109604430410174, -0.0001967837897097763, -0.00015011869680353805, -7.660119899681e-06, 4.9232420415918535e-05, 2.844973630311391e-05, -1.1509408420161798e-06, -9.343991089171552e-06, -4.399213561173833e-06, 3.8801299872766606e-07, 1.34735075769772e-06, 5.631963019736235e-07, -4.25918859259393e-08, -1.4422678652512737e-07, -5.943976881208015e-08, -3.979313524156213e-13, 1.0674713178958855e-08, 4.883922024683702e-09, 5.663853733044288e-10, -4.378176512587472e-10, -2.671788856819528e-10, -6.567170883081183e-11, -2.37118892503337e-13, 6.083971793664778e-12, 2.6033487386496266e-12, 6.57153500278571e-13, 1.1463976026446368e-13, 1.4136795975915104e-14, 1.1916107597695764e-15, 6.21457009644586e-17, 1.521935533823476e-18]

The largest coefficient is exactly $\frac{7281563888583524448979432097}{19807040628566084398385987584}$ which rounds to 0.3676250291566483 in 64-bit floating point representation.

enter image description here

$\endgroup$
7
  • 1
    $\begingroup$ That's very neat +1, even though it doesn't actually minimize stopband energy. $\endgroup$
    – Matt L.
    Feb 3, 2023 at 10:45
  • 1
    $\begingroup$ Hay Olli, you remember this question (and answers? $\endgroup$ Feb 9, 2023 at 6:47
  • 1
    $\begingroup$ @robertbristow-johnson yes, it's very much related $\endgroup$ Feb 9, 2023 at 6:50
  • $\begingroup$ I have a MATLAB program written by me and another guy a few years back, do you want it Olli? $\endgroup$ Feb 9, 2023 at 6:53
  • 1
    $\begingroup$ I don't wanna just post some MATLAB code (that I can't even remember exactly how I and this other guy, Barry, derived it). It was about fitting the sum of cosines to that smoothstep thingie and getting a maxflat FIR. $\endgroup$ Feb 9, 2023 at 7:04
2
$\begingroup$

WTF, i decided to post it anyway.

N = 36;                             % 0 <= N 

f = zeros(1, 4*N+3);
f(2*N+2) = 0.5;                     % midpoint tap

h = zeros(1, 4*N+3);
h(2*N+2) = 0.5;                     % midpoint tap

g = zeros(1, 4*N+3);
g(2*N+2) = 0.5;                     % midpoint tap


H = 0.25 * factorial(2*N+1)/factorial(N) * 4^(-N);

G = 0.25*4^(-N);                    % two ways to do the same thing
for nn = N+1:2*N+1
    G = G * nn;
end


for k = 0:N

    ff = ((-1)^k * gamma(N + 3/2)^2) / (pi * (2*k+1) * gamma(k + N + 2) * gamma(-k + N + 1));

    hh = 0;
    gg = 0;
    for n = k:N

%           hh = hh + (-0.25)^n * (factorial(2*n)/factorial(n)/factorial(n-k)) / ( factorial(N-n)*factorial(n+k+1) );
            hh = hh + (-0.25)^n * (gamma(2*n+1)/(gamma(n+1))/gamma(n-k+1)) / ( gamma(N-n+1)*gamma(n+k+2) );

            B = 1;
            if n > 0
                for nn = n+1:2*n
                    B = B * nn;
                end
            end
            gg = gg + (-0.25)^n *  ( ( B/factorial(n-k) ) / ( factorial(N-n)*factorial(n+k+1) ) );          % two ways to do (almost) the same thing
    
    end

    f(2*N+2 + (2*k+1)) = ff;
    f(2*N+2 - (2*k+1)) = ff;

    h(2*N+2 + (2*k+1)) = H*hh;
    h(2*N+2 - (2*k+1)) = H*hh;

    g(2*N+2 + (2*k+1)) = G*gg;
    g(2*N+2 - (2*k+1)) = G*gg;
end

figure(1)
plot(g + 0.1, '.g')
hold on
plot(h + 0.05, '.r')
plot(f, '.b')
hold off

figure(2)
plot( abs(fft(g, 32*length(g))), 'g' )
hold on
plot( abs(fft(h, 32*length(h))), 'r' )
plot( abs(fft(f, 32*length(f))), 'b' )
hold off
axis([0 32*length(h) 0 1])

figure(3)
plot( 20*log10(abs(fft(g, 32*length(g))) + 1e-40), 'g' )
hold on
plot( 20*log10(abs(fft(h, 32*length(h))) + 1e-40), 'r' )
plot( 20*log10(abs(fft(f, 32*length(f))) + 1e-40), 'b' )
hold off
axis([0 32*length(h) -400 0])

linear frequency linear gain: maxflat

linear frequency dB gain: maxflat dB

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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