The scipy library of matplotlib in Python 3 contains a function \"scipy.signal.f

Question

The scipy library of matplotlib in Python 3 contains a function "scipy.signal.freqz" which computes the frequency response of a digital filter. I am trying to write my own version of this function, to relay the same output as this imbedded function. I plan to use the inputs: freqz(b, a = 1, worN=None, whole=None): [b=numerator of linear filter; a=denominator of linear filter; worN=at none, compute at 512 frequencies equally spaced around the unit circle; whole=at none, frequency computed around upper half of unit circle, from 0 to pi]. This function should output two outputs: w and h

After implementing the function, it should pass all the tests provided in the last part of code.

DON'T USE FFT FUNCTION

I am going to provide the code that I am given.

Code:

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sig
import matplotlib as mpl
mpl.rcParams['figure.figsize'] = (5,3)
%matplotlib inline

# Implementation of freqz using different way.

def my_freqz(b, a = 1, worN=None, whole=None):
#replace with your answer
pass

Testing code:

from collections import namedtuple
Test = namedtuple('Test',['b','a','worN','whole'])
tests = [ Test([1,0.5,0.1], [1,1], 256, None),
Test(1,1,None, True),
Test(1,1,12, False),
Test([1,2],1,[0,1,2,3], None) ]

b,a = sig.iirdesign(0.2, 0.3, 1, 20)
tests.append(Test(b,a,None,None))

b = sig.firwin(63, 0.25)
tests.append(Test(b,1,1024,None))

for test in tests:
w, h = sig.freqz(test.b, test.a, test.worN, test.whole)
my_w, my_h = my_freqz(test.b, test.a, test.worN, test.whole)
print(np.allclose(w, my_w), np.allclose(h, my_h))

Tests should output TRUE for all.

Explanation / Answer

freqz(b, a = 1, worN=None, whole=None):
[b=numerator of linear filter;
a=denominator of linear filter;
worN=at none, compute at 512 frequencies equally spaced around the unit circle;
whole=at none, frequency computed around upper half of unit circle, from 0 to pi]
w and h
[w=ndarray, the normalized frequencies at which h was computed, in radians/sample;
h=ndarray, the frequency response]

def freqz(b, a=1, worN=None, whole=False):

#Initialize a
a = atleast_1d(a)

#Initialize b
b = atleast_1d(b)

#Check condition
if worN is None:

#Initialize worN
worN = 512

#Initialize
h = None

#Try block
try:

#Update
worN = operator.index(worN)

#Exception
except TypeError:

#Update
w = atleast_1d(worN)

#Otherwise
else:

#Check condition
if worN < 0:

#Throw an exception
raise ValueError('worN should be non negative, we got %s' % (worN,))

#Update
finalpoint = 2 * pi if whole else pi

#Update
w = np.linspace(0, finalpoint, worN, endpoint=False)

#Check condition
if (a.size == 1 and worN >= b.shape[0] and fftpack.next_fast_len(worN) == worN and (b.ndim == 1 or (b.shape[-1] == 1))):

#Update
nfft = worN if whole else worN * 2

#Check condition
if np.isrealobj(b) and np.isrealobj(a):

#Update
fftfunc = np.fft.rfft

#Otherwise
else:

#Update
fftfunc = fftpack.fft

#Update
h = fftfunc(b, n=nfft, axis=0)[:worN]

#Update
h /= a

#Check condition
if fftfunc is np.fft.rfft and whole:

#Update
stop = -1 if nfft % 2 == 1 else -2

#Update
hVal_flip = slice(stop, 0, -1)

#Update
h = np.concatenate((h, h[hVal_flip].conj()))

#Check condition
if b.ndim > 1:

#Update
h = h[..., 0]

#Update
h = np.rollaxis(h, 0, h.ndim)
#Delete
del worN

#Check condition
if h is None:

#Update
zmVal = exp(-1j * w)

#Update
h = (npp_polyval(zmVal, b, tensor=False)/npp_polyval(zmVal, a, tensor=False))

#Check condition
return w, h

def mfreqz(b,a, wp, ws):
w,h = freqz(b,a)
h_dB = 20 * np.log10 (abs(h))
ax1 = plt.subplot(211)
plt.plot(w/max(w),h_dB)
plt.vlines(wp, np.min(h_dB), np.max(h_dB), 'g', linewidth = 2)
plt.vlines(ws, np.min(h_dB), np.max(h_dB), 'r', linewidth = 2)
plt.ylim(np.min(h_dB), np.max(h_dB))
plt.ylabel('Magnitude (db)')
plt.xlabel(r'Normalized Frequency (x$pi$rad/sample)')
plt.title(r'Frequency response; order A: '+str(len(a))+' order B: '+str(len(b)))
plt.subplot(212, sharex = ax1)
h_Phase = np.unwrap(np.arctan2(np.imag(h),np.real(h)))
plt.plot(w/max(w),h_Phase)
plt.vlines(wp, np.min(h_Phase), np.max(h_Phase), 'g', linewidth = 2)
plt.vlines(ws, np.min(h_Phase), np.max(h_Phase), 'r', linewidth = 2)
plt.ylabel('Phase (radians)')
plt.xlabel(r'Normalized Frequency (x$pi$rad/sample)')
plt.title(r'Phase response')
plt.subplots_adjust(hspace=0.5)

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))

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)|$")

Related Questions

Navigate

• Browse (All)
• Browse T
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.