오디오 필터(3)-IIR 저역 통과 필터 (Low Pass Filter, LPF)
저역 통과 필터
저역통과 필터는 필터의 차단 주파수를 기준으로 이보다 낮은 주파수는 입력과 출력이 같도록 통과 시키고, 나머지 부분은의 주파수는 통과를 차단하는 필터입니다. 이 필터는 오디오용 필터의 기본이 되는 중요한 필터입니다. 이중에서 음성확성에 사용되어 좋은 음질을 제공하기도 합니다.
이 글에서는 저역 통과 필터의 필터 계수를 구하는 방법을 알아 보겠습니다. 필터 계수는 신호 처리 과정에서 IIR 필터에 의해 실제 입력되는 신호를 처리과정에 사용되기 때문입니다.
다른 글에서 IIR을 다루도록 하겠습니다.
Q를 이용한 저역 통과 필터의 설계
디지털 저역 통과 필터(low pass filter, LPF)를 구현하는 방법도 여러 가지가 있습니다. 여기서는 아날로그 필터의 Q값을 이용해서 필터를 구현하는 방법을 알아보겠습니다. 버터워스 2차의 경우 Q=0.71로 알려져 있습니다. 아날로그 회로에서 저역 통과 필터의 전달 함수는 다음과 같습니다.
$$H(s)=\frac{1}{s^{2}+\frac{s}{Q}+1}$$
필터 계수
위의 라플스 변환식은 z 변환 관계 테이블을 이용하면 디지털 필터용 전단 함수를 구할 수 있습니다. 중간 단계를 생략하고 바이쿼드의 필터 계수는 다음과 같이 정리됩니다. 여기서 $f=f_{0}/f_{s}$이며, $f_{0}$는 필터의 차단 주파수, $f_{s}$는 디지털 필터의 샘플링 주파수입니다. 이부분에 대햔 내용은 아래를 참조해 주세요
https://medialink.tistory.com/73?category=958130
https://medialink.tistory.com/77?category=958130
$$h=\frac{sin(2\pi f)}{2Q}$$
$$b_{0}=b_{1}=\frac{1-cos(2\pi f)}{2}$$
$$b_{1}=a_{1} = 1-cos(2\pi f)$$
$$a_{0} =1+h$$
$$a_{2} =1-h$$
최종적으로는 모든 필터 계수를 $a_{0}$로 나눈 값을 사용합니다. 그러면 $a_{0}=1$이기 때문에 총 6개의 필터 계수가 만들어 지지만 실제로 5개만 사용해도 되도록 연산량을 줄일 수 있습니다.
GNU Octave로 필터 검사
원하는 특성의 저역통과 필터로 구현되는지 그리고 이 필터가 안정적인지 검토할 필요가 있습니다. 이를 위해서 GNU Octave를 사용합니다. 이 프로그램은 matlab과 비슷하면서도 GNU로서 무료로 사용 가능합니다.
다음과 같이 lpf.m을 작성하고 실행하면 필터 주파수 f0를 기준으로 필터에 필요한 계수를 생성하고 이를 z영역에서 피터의 안정성을, 주파수 특성을 그래프로 볼 수 있습니다. 이 코드는 다른 필터를 검사할 때도 사용합니다.
clear all;
pkg load signal;
pkg load symbolic;
fs = 44100.0;
f0 =1000.0;
NFFT = 1024;
Q=0.71; # butterworth, 2 order
function [coef_b, coef_a] = LPF(f0,fs,Q)
coef_b = zeros(1,3);
coef_a = zeros(1,3);
w0 = 2 * pi * f0 / fs;
c0 = cos(w0);
s0 = sin(w0);
h = s0 /(2 * Q);
b0 = (1 - c0)/2;
b1 = 1 - c0;
b2 = b0;
a0 = 1 + h;
a1 = -2 * c0;
a2 = 1 - h;
coef_b(1) = b0 / a0;
coef_b(2) = b1 / a0;
coef_b(3) = b0 / a0;
coef_a(1) = 1;
coef_a(2) = a1 / a0;
coef_a(3) = a2 / a0;
endfunction
[b,a]=LPF(f0,fs,Q);
[h,w]=freqz(b,a,NFFT,fs);
m = 20.0*log10(abs(h));
figure(1);
clf;
zplane(b,a);
figure(2);
clf;
subplot(2,1,1)
semilogx(w,m)
xlim([20 44100/2])
ylim([-40 10])
xlabel("frequency(Hz)")
ylabel("magitude(dB)")
grid()
subplot(2,1,2)
semilogx(w,arg(h)*180/pi)
xlim([20 44100/2])
ylim([-200 200])
xlabel("frequency(Hz)")
ylabel("phase(deg)")
grid()필터의 특성 검토
z 영역에서 필터의 특성을 볼 수 있습니다 모든 극점(x)이 단위원 안에 있어 필터는 안정적인 상태입니다.
주파수 특성을 보면 1kHz를 기준으로 저역 통과 필터가 형성됨을 확인할 수 있습니다.
필터의 합성
고차의 필터는 아래 그림과 같이 1차와 2차 필터를 합성해서 만듭니다. 이때 하나의 필터를 스테이지(stage)라고 부릅니다.
필터의 차단 주파수와 필터의 형식 그리고 필터의 차수를 결정하면 다음과 같이 다수의 Q 값을 이용해서 최종 필터를 구합니다. Q를 쉽게 구하는 방법으로는 FilterPro를 이용하는 것입니다.
아래 표는 필터 차수와 필터 형태에 따라 몇개의 필터 스테이지와 Q가 필요한지 표시한 내용입니다.
| 필터차수 | 필터 형식 | Q[0] | Q[1] | Q[2] | Q[3] |
| 1 | 0.5 | - | |||
| 2 | Bessel | 0.58 | - | ||
| Butterworth | 0.71 | - | |||
| L-R | 0.5 | - | |||
| 3 | Bessel | 0.5 | 0.69 | - | |
| Butterworth | 0.5 | 1.0 | - | ||
| L-R | - | ||||
| 4 | Bessel | 0.52 | 0.81 | - | |
| Butterworth | 0.54 | 1.31 | - | ||
| L-R | 0.707 | 0.707 | - | ||
| 5 | Bessel | 0.5 | 0.56 | 0.92 | - |
| Butterworth | 0.5 | 0.62 | 1,62 | - | |
| L-R | - | ||||
| 6 | Bessel | 0.51 | 0.61 | 1.02 | - |
| Butterworth | 0.52 | 0.71 | 1.93 | - | |
| L-R | - | ||||
| 7 | Bessel | 0.5 | 0.53 | 0.66 | 1.13 |
| Butterworth | 0.5 | 0.55 | 0.8 | 2.25 | |
| L-R | - | ||||
| 8 | Bessel | 0.51 | 0.56 | 0.71 | 1.23 |
| Butterworth | 0.51 | 0.6 | 0.9 | 2.56 | |
| L-R | 0.54 | 0.54 | 1.31 | 1.31 | |
이를 이용해서 저역 통과 필터를 구하는 스크립트는 다음과 같습니다.
- 필터의 합성을 위해서 signal과 symbolic의 패키지를 사용합니다.
- 필터의 형태와 차수와 관련되어 Q값을 배열로 만듭니다.
- 필터의 차수에 따라서 필터 계수를 합성합니다.
- 필터의 특성을 그래프로 표시합니다.
## 필요한 패키지를 불러온다
pkg load signal;
pkg load symbolic;
## 기존의 변수를 모두 지운다.
clear all
## 심볼릭 변수를 선언하여, 전달함수를 수식처럼 구한다.
syms z Hd Hn dev num N;
Hd = 1
Hn = 1
## 샘플링 주파수
SAMPLING_RATE = 44100.0
## fft 크기
NFFT=1024
# Bessel filter의 Q
BESSEL_Q = [0.5, 0.0, 0.0, 0.0;
0.58, 0.0, 0.0, 0.0;
0.5, 0.69, 0.0, 0.0;
0.52, 0.81, 0.0, 0.0;
0.5, 0.56, 0.92, 0.0;
0.51, 0.61, 1.02, 0.0;
0.5, 0.53, 0.66, 1.13;
0.51, 0.56, 0.71, 1.23];
# Butterworth의 Q
BW_Q = [0.5, 0.0, 0.0, 0.0;
0.71, 0.0, 0.0, 0.0;
0.5, 1.0, 0.0, 0.0;
0.54, 1.31, 0.0, 0.0;
0.50, 0.62, 1.62, 0.0;
0.52, 0.71, 1.93, 0.0;
0.5, 0.55, 0.80, 2.25;
0.51, 0.60, 0.90, 2.56];
# L_R의 Q
LR_Q = [0.0, 0.0, 0.0, 0.0;
0.5, 0.0, 0.0, 0.0;
0.0, 0.0, 0.0, 0.0;
0.707, 0.707, 0.0, 0.0;
0.0, 0.0, 0.0, 0.0;
0.0, 0.0, 0.0, 0.0;
0.0, 0.0, 0.0, 0.0;
0.54, 0.54, 1.31, 1.31];
#################################################
# LPF 구하는 함수
#################################################
function [coef_a, coef_b]=lpf(f0, fc, Q)
w0 = 2.0 * pi * f0 / fc;
c0 = cos(w0);
s0 = sin(w0);
h = s0 /(2.0 * Q);
b0 = (1.0 - c0)/2.0;
b1 = 1.0 - c0;
b2 = b0;
a0 = 1.0 + h;
a1 = -2.0 * c0;
a2 = 1.0 - h;
coef_b(1) = b0 / a0;
coef_b(2) = b1 / a0;
coef_b(3) = b0;
coef_a(1) = 1.0;
coef_a(2) = a1 / a0;
coef_a(3) = a2 / a0;
endfunction
## 필터 차수
order = 8;
## 필터 차단 주파수
f0 = 500.0;
#################################################
# 필터의 전달함수를 구한다.
#################################################
# 무조건 4개의 Q에 대해 반복연산한다.
for i = 1 : 4
# 필터 차수와, 필터 종류에 따라 Q값을 가져온다.
Q = BW_Q(order,i);
#Q = BW_Q(order,i);
#Q = LR_Q(order,i);
# Q가 0 이면 필터 전달함수를 계산하지 않는다.
if (Q != 0)
# Q에 의해 LPF의 계수를 구한다.
[coef_a, coef_b] = lpf(f0,SAMPLING_RATE,Q);
# 필터 계수를 다항식으로 변경한다.
dev = poly2sym(coef_a,z);
num = poly2sym(coef_b,z);
# 기존의 필터 다항식과 현재 다항식을 곱해서 전달함수를구한다.
Hd *= dev;
Hn *= num;
endif
endfor
#################################################
# 필터 결과를 구한다.
#################################################
# 필터 전달함수에서 필터 계수만 추출
coef_a = sym2poly(Hd);
coef_b = sym2poly(Hn);
# 필터의 주파수 응답을 구함
[h,w] = freqz(coef_b,coef_a,NFFT,SAMPLING_RATE);
#################################################
# 극점 영점을 구한다.
#################################################
figure(1)
# 기존 그림 지우기
clf
# 극점 영점 표시
zplane(coef_b,coef_a)
#################################################
# 주파수 응답을 구한다.
#################################################
# 진폭응답을 로그로 처리
m =20.0*log10(abs(h));
figure(2)
# 기존 그림 지우기
clf
# 진폭응답 그리기
subplot(2,1,1)
semilogx(w,m)
xlim([20 44100/2])
ylim([-100 10])
xlabel("frequency(Hz)")
ylabel("magitude(dB)")
grid()
# 위상 응답 그리기
subplot(2,1,2)
semilogx(w,arg(h)*180/pi)
xlim([20 44100/2])
ylim([-200 200])
xlabel("frequency(Hz)")
ylabel("phase(deg)")
grid()광고좀 꾹 눌러주시면 고맙겠습니다.
참조 :
https://www.musicdsp.org/en/latest/Analysis/186-frequency-response-from-biquad-coefficients.html
https://www.musicdsp.org/en/latest/Filters/197-rbj-audio-eq-cookbook.html
위의 내용을 참조용으로만 사용해주세요, 무단 도용이나 무단 복제는 불허합니다.
기타 문의 사항은 gigasound@naver.com에 남겨 주시면 고맙겠습니다.