오디오 필터(13)-IIR 버터워스 저역 통과 필터(Butterworth Low Pass Filter) 설계

---

저역통과 필터

오디오 필터로 가장 유용한 형태는 버터워스 필터입니다. 통과대역의 특성이 가장 평탄하기 때문입니다. 물론 이를 발전시킨 필터도 있습니다만, 사용성이나 실용적으로 구현하기에도 버터워스 필터는 오디오에 적합합니다. 

다른 글에서는 필터의 Q를 이용해서 저역 통과 필터의 계수를 구했습니다. 

https://medialink.tistory.com/74?category=958130 

이글에서는 필터의 일반적인 규격 조건을 통해 버터워스의 s-영역(s-domain)(라플라스 영역)에서 구하고 이를 z-영역(z-domain)으로 변환해서 필요한 디지털 계수를 구하는 방법을 알아보겠습니다.

그리고 필터의 특성을 GNU Octave에서 알아보도록 하겠습니다.

---

버터워스 필터를 구하는 방식

아날로그 신호 영역인 s-영역에서 디지털 신호 영역인 z-영역으로 변환할 때 사용하는 방법은 다음과 같습니다. 여기서$T_{s}$는 샘플 시간으로 샘플률 $f_{s}$의 역수입니다. 또 여기서 $\Omega$는 아날로그 각주파수, $\omega$는 디지털 각주파수입니다. 

• IIM (Impulse Invariance Method) : $z=e^{sT_{s}}$관계를 이용해서 두 영역을 변환합니다.  이때 $\omega =\Omega /T_{s}$의 관계를 가집니다. LPF, HPF를 매우 정교하게 만들 때 유용한 방식입니다.
• BLT(Binear Transform Method) : $s=\frac{2}{T_{s}}\left (  \frac{z-1}{z+1}\right )$관계를 이용해서 두영 역을 변환합니다. 그런데 이 관계식은 근사 변환이며, 주파수 또한 $\Omega =\frac{2}{T_{s}}tan\left ( \frac{\omega }{2} \right )$로 변환해야 합니다.

 

이 글에서는 BLT방식을 사용해 보겠습니다. 이 방식은 비선형이며, 다소의 필터 오류가 있지만 모든 필터를 만들기에 적합한 방식입니다.

---

BLT 방식으로 필터 만들기

BLT 방식으로 필터를 만들려면 다음 순서대로 하면 됩니다.

• 필터 규격을 결정합니다.
  • 샘플링 시간 : $T_{s}$
  • 통과 대역 각주파수 : $\Omega_{p}, \omega_{p}$
  • 저지 대역 각주파수 : $\Omega_{s}, \omega_{s}$
  • 통과 대역 이득 : $\delta_{p}$
  • 저지 대역 이득 : $\delta_{s}$

• 각 대역의 각주파수를 상호 변환해서 구해냅니다.
• 필터 차수 (N)을 위의 조건을 이용해서 구합니다. 여기서 ceil()은 실수를 정수로 올림 하는 함수입니다. 

$$N=ceil\left (\frac{1}{2}\frac{log\frac{\frac{1}{\delta_{s}^{2}}-1}{\delta_{p}^{2}-1}}{log\frac{\Omega_{s}}{\Omega_{p}}}  \right )$$

• 그리고 다음과 같이 필터 차단 주파수를 구합니다. (LPF 형태입니다.)

$$\Omega_{c}=\frac{\Omega_{p}}{\left(\delta_{p}^{2} -1 \right )^\frac{1}{2N}}$$

• 필터 차수와 차단 주파수를 이용해서 모든 극점을 구합니다. 여기서 K=0... N-1입니다. 

$$s_{K}=\pm \Omega_{c}e^{j(N+2K+1)\pi/2N}$$

여기서 s-domain의 좌측에 있는 극점만 필터용을 사용합니다. 이는 복소수는 $s_{K}$의 실수 부분이 음수인 조건을 사용하면 됩니다.

• 다음과 같이 s-영역에서 필터 전달 함수를 설정합니다. 이는 버터워스 필터의 형성 조건입니다. 즉 모든 극점을 조건을 분모로 하고 필터의 이득 부분을 분자로 정합니다.

$$H(s)=\frac{\Omega_{c}^{N}}{\prod_{k=0}^{N-1}(s-s_{k})}$$

• 그런 다음의 변환 관계를 이용해서 s-영역의 필터 전달 함수를 z-영역으로 만듭니다.

 $$s=\frac{2}{T_{s}}\left (  \frac{z-1}{z+1}\right )$$

• 마지막으로 전달 함수에서 필터 계수 부분을 추출하면 됩니다.

---

예제

필터 규격을 다음과 같이 정합니다.  

• 샘플링 시간 : $T_{s}=1.0/44100.0$
• 통과 대역 각주파수 : $\omega_{p}=0.2\pi$
• 저지 대역 각주파수 : $\omega_{s}=0.6\pi$
• 통과 대역 이득 : $\delta_{p}=0.8$
• 저지 대역 이득 : $\delta_{s}=0.2$

위의 조건으로 필터 계수를 구하는 GNU Octave용 스크립트의 내용은 다음과 같습니다.

조건에 의해 필터 차수 N=2이기 때문에 필터의 분모 부분은 coeff_num에 분자 부분은 coeff_dev에  1x3차원 행랼에 기록됩니다. 마지막 부분은 필터의 특성을 그래프로 알아보는 내용입니다.

###
# butterworth filter using BIT
###
## prepare memory
clc;
clear all;
## FFT
NFFT = 1024;
## load package
pkg load signal;
pkg load symbolic;
## symbol
# z-domain symbol
z = sym('z');
# s-domain symbol
s = sym('s');
## pole number symbole
K = sym('K');
## sampling rate
fs = 44100.0;
# sample time
Ts = 1/fs;
## linear filter attueation
delta_p=0.8;
delta_s= 0.2;
## digital filter radian frequency
w_p=0.01*pi;
w_s=0.05*pi;
## analog filter radian Frequency
# using BIT method
Omega_p = (2/Ts)*tan(w_p/2);
Omega_s = (2/Ts)*tan(w_s/2);
## filter order
N_num = log10(((1/(delta_s*delta_s))-1)/((1/(delta_p*delta_p))-1));
N_dev = log10((Omega_s)/(Omega_p));
#decision filter order
N=ceil((1/2)*N_num/N_dev);
## filter cuttoff frequency
Omega_c = Omega_p/((1/(delta_p*delta_p))-1)^(1/2/N);
## filter pole K= 0 ... N-1
M=double(N)-1
## s-domain trasfer function H_s
# numerator
H_num = Omega_c^N;
H_dev = 1;
for k =0 : M
  x = Omega_c*exp(j*(N+(2*k)+1)*pi/2/N);
  if (real(x) < 0)
    H_dev *= (s-x)
  endif
  x = -Omega_c*exp(j*(N+(2*k)+1)*pi/2/N);
  if real(x) < 0
    H_dev *= (s-x)
  endif
endfor
H_s = H_num/ H_dev;
# z-domain, transfer function, bilinear 
# substitution
H_z = subs(H_s,s,(2/Ts)*((z-1)/(z+1)));
# fuction expansion
H_z = simplify(expand(H_z));
# coefficient extraction
[num_z,den_z]= numden(H_z);
coeff_b = coeffs(num_z,z,'all');
coeff_a = coeffs(den_z,z,'all');
coeff_num = double(coeff_b);
coeff_dev = double(coeff_a);
# Coefficient normalization
a0 =  coeff_dev(1);
coeff_num /= a0
coeff_dev /= a0

[h,w] = freqz(coeff_num,coeff_dev,NFFT,fs);
m = 20.0*log10(abs(h));

figure(1);
clf;
zplane(coeff_num,coeff_dev);
figure(2);
clf;
subplot(2,1,1)
semilogx(w,m)
xlim([20 fs/2])
ylim([-40 10])
xlabel("frequency(Hz)")
ylabel("magitude(dB)")
grid()
subplot(2,1,2)
semilogx(w,arg(h)*180/pi)
xlim([20 fs/2])
ylim([-200 200])
xlabel("frequency(Hz)")
ylabel("phase(deg)")
grid()

butter_blt.m

0.00MB

---

필터 특성의 검토

z-영역에서 필터의 특성은 다음과 같습니다. 모든 극점이 단위원 안에 있어 안정적입니다. 

 

주파수 영역의 특성은 다음과 같습니다. 저주파 필터가 구성되어 있습니다.

---

광고좀 꾹 눌러주시면 고맙겠습니다. 

---

위의 내용을 참조용으로만 사용해주세요. 무단 도용이나 무단 복제는 불허합니다.

기타 문의 사항은 gigasound@naver.com에 남겨 주시면 고맙겠습니다.