1
$\begingroup$

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:

plot of frequency response using freqz()

plot of frequency response using dfreqresp()

plot of magnitude and phase error between freqz and dfreqresp

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.

$\endgroup$
5
  • $\begingroup$ Saluton DrEsperanto! Mi estas denaske Esperantisto. Bonvenon al DSP.SE! $\endgroup$ Feb 10, 2022 at 12:17
  • $\begingroup$ @DanBoschen Ĉu vere? Dankon pro la bonveno! Mi neniam divenus ke mi hazarde renkontus denaskulon en forumo, kiu ne temas pri Esperanto mem ;) Mojose! $\endgroup$ Feb 10, 2022 at 19:27
  • $\begingroup$ Mia patro: scua.library.umass.edu/boschen-allan-c Li parolis al mi nur en Esperanto gxis mi estis 13-jara, sed nun mi ne multe parolas gxin. $\endgroup$ Feb 10, 2022 at 19:39
  • $\begingroup$ @Dan Boschen La mondo vere estas malgranda. Mi pensas ke vi renkontis mian koramikinon antaŭ kelkaj jaroj ĉe antaŭfunebra ceremonio en Beverly, MA. Ŝi rekonis vian verdan stelon. Kia eta mondo. $\endgroup$ Feb 13, 2022 at 3:34
  • $\begingroup$ Ĉu vere! Jes ja, mi memoras tion kaj ŝin. La filino de Rikardo. Tio estas mirinda. $\endgroup$ Feb 13, 2022 at 6:32

1 Answer 1

1
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ Interesting, I was thinking that because these two functions are mathematically equivalent they would have to have the same response (since I multiplied by z/z = 1). Does this boil down to a change in filter implementation (causal vs non-causal)? $\endgroup$ Feb 10, 2022 at 19:18

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.