2016-06-15 3 views
7

Я изучаю обработку цифровых сигналов для реализации фильтров и использую python для простого внедрения тестовых идей. Поэтому я просто начал использовать библиотеку scipy.signal, чтобы найти импульсную характеристику и частотную характеристику различных фильтров.Частотный отклик Scipy.signal

В настоящее время я работаю над книгой «Цифровые сигналы, процессоры и шум» Пола А. Линна (1992) »(и нахождение этого удивительного ресурса для изучения этого материала). В этой книге они имеют фильтр с передаточными функциями показано ниже:

Я разделил числитель и знаменатель на , чтобы получить следующее уравнение:

Я тогда осуществил это с помощью Scipy, используя:

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

и произвести Участок показано ниже:

Plot showing Frequency response calculated by Scipy.sigal.freqz

Однако в книге показана частотная характеристика будет следующее:

Plot showing Frequency response from Book referenced above

Они являются такой же формы, но отношение пиков при ~ 2,3 и 0.5 очень разные для двух сюжетов, может кто-нибудь предположить, почему это так?

Edit:

Чтобы добавить к этому, я просто реализовал функцию для вычисления частотной характеристики вручную (путем вычисления расстояния от полюсов и нулей функции), и я получаю подобное отношение к сюжет, созданный scipy.signal, однако цифры не совпадают, кто-нибудь знает, почему это возможно?

Реализация выглядит следующим образом:

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

и получить следующий сюжет:

Plot resulting from manual calculation of the frequency response.

ответ

5

SciPy генерируемый частотный отклик правильный. В любом случае, я бы не стал доверять фигуре книги, которая, как представляется, была взята вручную.

Если вы хотите найти частотную характеристику «вручную», это может быть просто сделано путем определения функции возвращения оригинального Z-преобразованием и его оценки на единичной окружности следующего

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

Результат точно так же, как SciPy.

Смежные вопросы