Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Frequency Response Scipy.signal

I'm learning digital signal processing to implement filters and am using python to easily implement a test ideas. So I just started using the scipy.signal library to find the impulse response and frequency response of different filters.

Currently I am working through the book "Digital Signals, Processors and Noise by Paul A. Lynn (1992)" (and finding it an amazing resource for learning this stuff). In this book they have a filter with the transfer functions shown below:

I divided the numerator and denominator by in order to get the following equation:

I then implemented this with Scipy using:

NumeratorZcoefs = [1, -1, 1, -1]
DenominatorZcoefs = [1, 0.54048, -0.62519, -0.66354, 0.60317, 0.69341]

FreqResponse = scipy.signal.freqz(NumeratorZcoefs, DenominatorZcoefs)
fig = plt.figure(figsize = [8, 6])
ax = fig.add_subplot(111)
ax.plot(FreqResponse[0], abs(np.array(FreqResponse[1])))
ax.set_xlim(0, 2*np.pi)
ax.set_xlabel("$\Omega$")

and produce the plot shown below:

Plot showing Frequency response calculated by Scipy.sigal.freqz

However in the book the frequency response is shown to be the following:

Plot showing Frequency response from Book referenced above

They are the same shape but the ratio of the peaks at ~2.3 and 0.5 are very different for the 2 plots, could someone suggest why this is?

Edit:

To add to this, I've just implemented a function to calculate the frequency response by hand (by calculating the distance from the poles and zeros of the function) and I get a similar ratio to the plot generated by scipy.signal, however the numbers are not the same, does anyone know why this might by?

Implementation is as follows:

def H(omega):
    z1 = np.array([0,0]) # zero at 0, 0
    z2 = np.array([0,0]) # Another zero at 0, 0
    z3 = np.array([0, 1]) # zero at i
    z4 = np.array([0, -1]) # zero at -i
    z5 = np.array([1, 0]) # zero at 1

    z = np.array([z1, z2, z3, z4, z5])

    p1 = np.array([-0.8, 0])
    p = cmath.rect(0.98, np.pi/4)
    p2 = np.array([p.real, p.imag])
    p = cmath.rect(0.98, -np.pi/4)
    p3 = np.array([p.real, p.imag])
    p = cmath.rect(0.95, 5*np.pi/6)
    p4 = np.array([p.real, p.imag])
    p = cmath.rect(0.95, -5*np.pi/6)
    p5 = np.array([p.real, p.imag])

    p = np.array([p1, p2, p3, p4, p5])

    a = cmath.rect(1,omega)
    a_2dvector = np.array([a.real, a.imag])

    dz = z-a_2dvector
    dp = p-a_2dvector

    dzmag = []
    for dis in dz:
           dzmag.append(np.sqrt(dis.dot(dis)))

    dpmag = []
    for dis in dp:
           dpmag.append(np.sqrt(dis.dot(dis)))        

    return(np.product(dzmag)/np.product(dpmag))

I then plot the frequency response like so:

omegalist = np.linspace(0,2*np.pi,5000)
Hlist = []

for omega in omegalist:
    Hlist.append(H(omega))

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(omegalist, Hlist)
ax.set_xlabel("$\Omega$")
ax.set_ylabel("$|H(\Omega)|$")

and get the following plot:

Plot resulting from manual calculation of the frequency response.

like image 926
SomeRandomPhysicist Avatar asked Jun 15 '16 11:06

SomeRandomPhysicist


1 Answers

The SciPy generated frequency response is correct. In any case, I wouldn't trust the book's figure which appears to have been drawn by hand.

If you want to find the frequency response "manually", this can be simply done by defining a function returning the original Z-transform and evaluating it on the unit circle as follows

def H(z):
    num = z**5 - z**4 + z**3 - z**2
    denom = z**5 + 0.54048*z**4 - 0.62519*z**3 - 0.66354*z**2 + 0.60317*z + 0.69341
    return num/denom

import numpy as np
import matplotlib.pyplot as plt

w_range = np.linspace(0, 2*np.pi, 1000)
plt.plot(w_range, np.abs(H(np.exp(1j*w_range))))

The result is exactly the same as SciPy.

like image 137
Stelios Avatar answered Nov 15 '22 21:11

Stelios