2016-06-15 16 views
7

私はフィルターを実装するためにデジタル信号処理を学んでおり、簡単にテストアイデアを実装するためにPythonを使用しています。そこで私はscipy.signalライブラリを使用して、異なるフィルタのインパルス応答と周波数応答を見つけました。周波数応答Scipy.signal

現在、私は「Digital Signals、Processors and Noise by Paul A. Lynn(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が2つのプロットでは非常に異なるため、誰かがこれがなぜ示唆できるのでしょうか?

編集:

これに追加するには、私は(関数の極とゼロからの距離を計算して)手で周波数応答を計算する機能を実装していると私は同様の比率を取得します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のダウンロードを生成し、次のように

実装です周波数応答が正しい。いずれにせよ、私は手で描かれているような本の姿を信じていません。

あなたは「手動」周波数応答を検索したい場合は

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

を次のように、これは単純に、元のz変換を返す関数を定義し、単位円上で、それを評価することによって行うことができる結果がありますSciPyとまったく同じです。