Can someone explain why signal.freqz() generates a different phase in the frequency fresponse than scipy.signal.dfreqresp()?

Question

I have been using scipy to analyze filter performance for a single-pole IIR filter, and I noticed a disagreement between the phase in the outputs of the freqz() function compared to the output of the dfreqresp() function for the same system.

The transfer function for this filter is: $$H(z) = \frac{b}{1-az^{-1}}$$

According to the documentation that is the format freqz() expects: $$H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2} + \dots + b_Nz^{-N}}{a_0 + a_1z^{-1} + a_2z^{-2} + \dots + a_Mz^{-M}}$$

but TransferFunction() expects the form with positive powers of z: $$H(z) = \frac{b_0z^{N} + b_1z^{N-1} + \dots + b_{N-1}z + b_N}{a_0z^{M} + a_1z^{M-1} + \dots + a_{M-1}z + a_M}$$

which I get by just multiplying by $\frac{z}{z}$: $$H(z) = \frac{b}{1-az^{-1}}\cdot\frac{z}{z} = \frac{bz}{z-a}$$

Luckily in this case the resulting coefficient arrays are identical ([b], [1, -a]), so that is not my problem.

Here is a short script to illustrate what I am talking about:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz, TransferFunction, dfreqresp, bode

def plot_freq_response(title, w, h, fs):
    "Plots the frequency response with magnitude in dB and phase in degrees"
    # Convert frequency, magnitude, and phase
    f = w * fs / (2*np.pi)
    h_db = 20 * np.log10(np.abs(h))
    h_phase = np.angle(h) * 180 / np.pi

    # Plot info about the cutoff frequency
    cutoff_idx = abs(h_db-(-3)).argmin()
    cutoff_freq = f[cutoff_idx]
    cutoff_mag = h_db[cutoff_idx]
    cutoff_phase = h_phase[cutoff_idx]
    cutoff_freq_str = f"Cutoff = {cutoff_freq:.2f} Hz"
    cutoff_phase_str = f"Phase = {cutoff_phase:.2f} degrees"


    # Magnitude
    plt.figure(figsize=(10,7))
    plt.subplot(2,1,1)
    plt.title(title)
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Gain [dB]")
    plt.semilogx(f, h_db)
    plt.scatter(cutoff_freq, cutoff_mag, color='red', label=cutoff_freq_str)
    plt.grid()
    plt.legend()

    # Phase
    plt.subplot(2,1,2)
    plt.ylabel("Phase [degrees]")
    plt.xlabel("Frequency [Hz]")
    plt.semilogx(f, h_phase)
    plt.scatter(cutoff_freq, cutoff_phase, color='red', label=cutoff_phase_str)
    plt.grid()
    plt.legend()

    plt.show()

# Filter parameters:
fs = 1000       # Sampling rate (Hz)
a = 0.75        # IIR filter decay value
b = 1 - a

# Use the same frequency range for both plots to make them comparable
w = np.linspace(0, np.pi, 1024)

# Generate the frequency response using freqz():
w, h_freqz = freqz([b], [1, -a], worN=w)

plot_freq_response(f"Low-Pass Frequency Response with decay={a:.2f} using freqz",w, h_freqz, fs)


# Now do the same with dfreqresp() which uses dlti() which uses TransferFunction():
s = TransferFunction([b], [1, -a], dt=1/fs)
w, h_dfreqresp = dfreqresp(s, w)

plot_freq_response(f"Low-Pass Frequency Response with decay={a:.2f} using dfreqresp", w, h_dfreqresp, fs)

# Plot Magnitude and phase errors
plt.figure(figsize=(10,4))
plt.xlabel("Frequency [Hz]")
plt.ylabel("Error")
f = w * fs / (2*np.pi)
mag_error = np.abs(np.abs(h_freqz) - np.abs(h_dfreqresp))
phase_error = np.abs(np.angle(h_freqz) - np.angle(h_dfreqresp)) * 180/np.pi
plt.plot(f, mag_error, label="Magnitude Error")
plt.plot(f, phase_error, label="Phase Error (degrees)")
plt.legend()
plt.show()

Here are the plots generated by the code above:

The phase of the freqz() plot clearly wraps back up to 0 degrees, while the dfreqresp() plot does not.

One interesting finding is that the magnitudes are the same (within error), but the difference between the phase in each plot is linear from 0 to 180 degrees as frequency goes from 0 to $\frac{fs}{2}$.

Can someone please explain why this is the case? I have run similar comparisons using both Octave and Matlab, and they behave the same way, so I think the error here is in my understanding of the differences between these two functions.

Answer 1

Let's use $H_1(z)$ to denote the transfer function with negative powers of $z$, and $H_2(z)$ is the one with positive powers of $z$. If $M=N$ we have $H_1(z)=H_2(z)$. However, in general we have

$$H_2(z)=z^{N-M}H_1(z)\tag{1}$$

For the frequency responses this means

$$H_2(e^{j\omega})=e^{j\omega (N-M)}H_1(e^{j\omega})\tag{2}$$

I.e., their magnitudes are identical, and their phases differ by a linear phase term. In your example, $N=0$ and $M=1$, so $H_2(z)$ has an additional linear phase term $e^{-j\omega}$.