SciPy를 활용하여 필터 계수 설계하는 법

물론입니다. 아래는 scipy.signal 라이브러리의 cheb2ord()와 cheby2() 함수를 사용해서, passband와 stopband 조건을 충족하는 Chebyshev-II type low-pass 필터를 설계하는 예입니다. 이 함수들의 보다 자세한 설명은 SciPy 문서를 참고하세요 (https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.cheb2ord.html).

#!/usr/bin/env xmodelpy

import scipy.signal
import numpy as np

# design of low-pass chebyshev-II filter
wp = 2*np.pi * 0.9e6    # passband edge frequency (rad/s)
ws = 2*np.pi * 1.1e6    # stopband edge frequency (rad/s)
rp = 3                  # maximum attenuation within pass band in dB
rs = 40                 # minimum attenuation within stop band in dB

N, wn = scipy.signal.cheb2ord(wp, ws, rp, rs, analog=True)
z,p,k = scipy.signal.cheby2(N, rs, wn, 'lowpass', analog=True, output='zpk')

그리고 아래는 Chebyshev-II type의 band-pass 필터를 설계하는 예입니다.

# design of bass-pass chebyshev-II filter
wp = [2*np.pi*1.1e6, 2*np.pi*2.9e6]     # passband edge frequency (rad/s)
ws = [2*np.pi*0.9e6, 2*np.pi*3.1e6]     # stopband edge frequency (rad/s)
rp = 3                                  # maximum attenuation within pass band in dB
rs = 40                                 # minimum attenuation within stop band in dB

N, wn = scipy.signal.cheb2ord(wp, ws, rp, rs, analog=True)
z,p,k = scipy.signal.cheby2(N, rs, wn, 'bandpass', analog=True, output='zpk')

비슷한 방법으로 SciPy 라이브러리의 cheb1ord()와 cheby1() 함수를 사용해서 Chebyshev-I 타잎의 필터를, buttord()와 butter() 함수를 사용해서 Butterworth 타잎 필터를, 그리고 ellipord()와 ellip() 함수를 사용해서 elliptic 타잎 필터를 각각 설계할 수 있습니다. 위키피디아 문서에는 이들 선형 필터들의 차이가 잘 설명되어 있으니 참고하시기 바랍니다 (https://en.wikipedia.org/wiki/Chebyshev_filter).

위의 예제에서 z, p, k 변수의 결과값은 각각 zero들의 주파수, pole들의 주파수, 그리고 gain값에 해당합니다. 다만, 이들 값을 XMODEL의 filter primitive에서 사용하는 'gain', 'poles', 'zeros' 파라메터들의 값으로 사용하기 위해서는 아래의 몇가지 사항들을 고려하셔야 합니다.

1. XMODEL filter primitive의 poles와 zeros 파라메터에서는 pole과 zero의 주파수를 Hertz 단위로 표현합니다. SciPy의 결과값은 rad/s 단위로 표현되어 있습니다.
2. XMODEL filter primitive에서는 pole과 zero의 주파수의 실수부 부호가 반대입니다. 즉, XMODEL은 left-half plane (LHP) pole/zero의 실수부를 양수로, right-half plane (RHP) pole/zero의 실수부를 음수로 표현하는데 반해, SciPy는 LHP pole/zero의 실수부를 음수로, RHP pole/zero의 실수부를 양수로 표현합니다.
3. XMODEL filter primitive의 poles와 zeros 파라메터에서는 각 pole 또는 zero 주파수의 실수부와 허수부를 아래처럼 별도 항목으로 나열합니다. 이는 SystemVerilog가 복소수를 직접 표현하지 못하기 때문입니다.

poles = '{real(fp1), imag(fp1), real(fp2), imag(fp2), ..., real(fpM), imag(fpM)}
zeros = '{real(fz1), imag(fz1), real(fz2), imag(fz2), ..., real(fzN), imag(fzN)}

4. XMODEL filter primitive의 poles와 zeros 파라메터에 나열된 복소수 값의 pole/zero들은 반드시 정확한 conjugate 쌍을 이루어야 합니다.
5. XMODEL filter primitive에서는 zero가 전혀 없는 전달함수를 표현하기 위해서 zeros 파라메터를 '{0.0}으로 정의합니다.
6. XMODEL filter primitive는 strictly proper한 시스템만 구현할 수 있기 때문에, pole의 수가 항상 zero의 수보다 많아야 합니다. 이 조건을 만족시키기 위해, 충분히 높은 주파수를 가져 필터 특성에 영향을 주지 않는 dummy pole을 추가할 수도 있습니다 (예를 들면 1THz pole).

예를 들면, 아래의 z, p, k 값들이 구해지면:

z = [ -3.72529030e-09+16418754.44763247j  -5.23868948e-09-16418754.44763249j
        1.92618690e-09 +6800870.76948295j   7.78721962e-11 -6800870.76948295j]
p = [ -1075430.76302337+2991438.64240786j -3170099.69183802+1512931.2492777j
      -3170099.69183801-1512931.2492777j  -1075430.76302337-2991438.64240787j]
k = (0.01+0j)

이를 아래의 gain, poles, zeros 파라메터 값을 갖는 filter primitive 인스턴스로 변환할 수 있습니다:

filter #(.gain(0.01),
    .poles('{171160, 476102, 504537, 240790, 504537, -240790, 171160, -476102, 1e12, 0.0}),
    .zeros('{0.0, 2613125.9, 0.0, -2613125.9, 0.0, 1082392.2, 0.0, -1082392.2})) 
    filter_inst (.in(in), .out(out));

복소수 pole 또는 zero들이 conjugate 쌍끼리 모이도록 배치를 좀 바꾸었습니다. 그리고, strictly proper한 시스템을 위해서, 다른 pole/zero 주파수들보다 월등히 높은 주파수인 1THz에 dummy pole을 추가하였습니다.

아래 Python 스크립트는 z, p, k 값을 XMODEL filter primitive의 gain, poles, zeros 파라메터 값으로 변환하는 한가지 예입니다.

    def sort_pz(p):
        res = []
        for w in sorted(p, key=np.real):
            f = w / (2*np.pi)
            if f.imag == 0.0:
                res += [f.real, 0.0]
            elif f.imag > 0.0:
                res += [f.real, f.imag, f.real, -f.imag]
        return res

    gain = (k / np.prod(-p) * np.prod(-z)).real
    poles = sort_pz(-p)
    zeros = sort_pz(-z)
    if len(p) <= len(z): poles += [1e12, 0.0] * (len(z)-len(p)+1)
    if len(z) == 0: zeros += [0.0]

    print("gain  : %g" % gain)
    print("poles : %s" % ', '.join(['%g' % v for v in poles]))
    print("zeros : %s" % ', '.join(['%g' % v for v in zeros]))