[FIXED] Sign issues with firwin2 and signs

Issue

I suspected this wouldn't be a complicated issue, but, as is, I haven't found a solution:

• I want to design a (digital) filter using firwin2 from scipy.signal.
• The design works perfectly well...
• I then want to plot the results using freqz, and I'm interested in seeing the absolute gain. Using abs (or np.abs) is not possible, as, by definition, it returns only positive values (the filter contains a part which has negative gain (between f=1.0 and f=2.25):

• I plotted the phase of the output (phase), and as expected, it wraps around ±π. Using unwrap, I can get the phase 'straightened', with smaller jumps at the (gain) sign changes (π, indicating the sign change).

I could detect those jumps to change the sign, but surely there must be a more, hmm, scientific way do this? And, considering both remaining jumps are 'up', what sign do they represent?)

Here's the code to generate this:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
#  demo_sign_issue.py
#
#  Copyright 2022 John Coppens <[email protected]>
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
#  MA 02110-1301, USA.
#
#


import pylab as plt
import numpy as np
import scipy.signal as scs

def design_filter(f, g):
    taps = scs.firwin2(31, f, g)
    return scs.freqz(taps)


def plot(w, h):
    plt.subplot(2, 1, 1)
    plt.plot(w, np.abs(h), label = 'Gain')
    plt.ylabel('Gain')
    plt.grid(True)

    plt.subplot(2, 1, 2)
    plt.plot(w, np.angle(h), label = 'Angle')
    plt.plot(w, np.unwrap(np.angle(h)), label = 'Unwrapped angle')
    plt.ylabel('Angles (rad)')
    plt.grid(True)

    plt.legend()
    plt.show()


def main(args):
    w, h = design_filter(
                (0.0, 0.1, 0.5, 1.0),      # Node frequencies
                (0.0, 0.5, -0.4, 0.5))     # Node gains

    plot(w, h)
    return 0

if __name__ == '__main__':
    import sys
    sys.exit(main(sys.argv))

NOTE: In fact, to add to the confusion: The jump at f=1.0 seems to be special case: It is a combination of a +2π jump from the filter output and a -1π jump from the sign change.

Solution

Given your comment

I could detect those jumps to change the sign, but surely there must be a more, hmm, scientific way do this?

you might find this answer less than satisfactory (or redundant, if you've already implemented a similar calculation), but for what it's worth, here's a script that computes the "signed" magnitude response by estimating where the frequency response h(omega) goes through 0, and changes the sign of the graph of the magnitude at those points. The script generates this plot:

Here's the script. The function get_signs(h) estimate where h crosses the origin, and assigns either +1 or -1 to each value in h accordingly.

import numpy as np
from scipy.signal import freqz, remez, firwin2
import matplotlib.pyplot as plt


def get_signs(h):
    # This function returns an array of +/-1 with the same length as h.
    # Figure out the appropriate signs of the magnitude response by
    # finding where consecutive values in h jump either between the
    # first and third quadrants, or the second and fourth.  These
    # are the candidates for where the response h(omega) goes through
    # the origin, and those are where the phase changes by +/-pi,
    # and therefore where the sign of the "signed magnitude" changes.
    # For this approach to work, the frequencies should be sampled
    # relatively densely, to avoid missing the relevent changes in h.
    sr = np.sign(h.real)
    si = np.sign(h.imag)
    q = np.concatenate(([False], (sr[1:] != sr[:-1]) & (si[1:] != si[:-1])))
    s = 2*(1 - np.cumsum(q) & 1) - 1
    return s


taps = firwin2(31, [0, 0.1, 0.5, 1], [0.0, 0.5, -0.4, 0.5])
# taps = remez(35, [0, 0.3, 0.35, 1], [1, 0], fs=2)
# taps = remez(137, [0, 0.2, 0.3, 1], [0, 1], fs=2)

w, h = freqz(taps, worN=8000, fs=2)

s = get_signs(h)

ax = plt.subplot()
ax.plot(w, np.abs(h), alpha=0.6, linewidth=3.5,
        label='abs(h)')
ax.plot(w, np.abs(h)*s, 'k--', alpha=0.75, linewidth=1.5,
        label='"signed" abs(h)')

ax.set_xlabel('Frequency')
ax.set_ylabel('Gain')
ax.legend(framealpha=1, shadow=True)
ax.grid(True)

plt.show()
# plt.savefig('signed_amplitude_response.png')

The code requires a sufficiently dense sample of h to ensure that points where the sign changes are not missed. If you run this with filters that have many more taps, you might have to increase the argument worN to get a correct graph. It is also possible that there are degenerate cases where this method just won't work.

The code assume that the sign is positive at omega=0, which can result in an unexpected flip of the response; try the example taps = remez(137, [0, 0.2, 0.3, 1], [0, 1], fs=2) (commented out in the script) to see what I mean.

I was going to include a computation of the fully unwrapped phase response, but in fact there is no need. Once you unwrap the "usual" 2π jumps and then remove the π jumps that correspond to the sign changes, what is left is a straight line; i.e. you are left with the linear phase response that is such a great feature of symmetric FIR filters.

Answered By - Warren Weckesser