超音波流体屋のプログラム備忘録
バターワースフィルタ定数
最終更新:
usapfrog
-
view
Matlab / Octave のbutter関数を手実装したかったので。
アルゴリズムのみ記述する。理論的背景は参考URLを参照すること。
アルゴリズムのみ記述する。理論的背景は参考URLを参照すること。
バターワース多項式
伝達関数$H(s)$は下記のバターワース多項式の逆数である[1]。(mは整数)
$\L B(s) = \prod_{k=1}^{n/2} \left[ s^2 -2s \cos\left(\frac{2k+n-1}{2n}\pi\right) +1 \right]$ (for n=2m)
$\L B(s) = (s+1)\prod_{k=1}^{(n-1)/2} \left[ s^2 -2s \cos\left(\frac{2k+n-1}{2n}\pi\right) +1 \right]$ (for n=2m+1)
例えば5次では下記の通り。
$B(s) = (s+1)\left[s^2 -2s\cos(6\pi/10)+1 \right] \left[ s^2 -2s\cos(8\pi/10)+1 \right]$
$\L B(s) = \prod_{k=1}^{n/2} \left[ s^2 -2s \cos\left(\frac{2k+n-1}{2n}\pi\right) +1 \right]$ (for n=2m)
$\L B(s) = (s+1)\prod_{k=1}^{(n-1)/2} \left[ s^2 -2s \cos\left(\frac{2k+n-1}{2n}\pi\right) +1 \right]$ (for n=2m+1)
例えば5次では下記の通り。
$B(s) = (s+1)\left[s^2 -2s\cos(6\pi/10)+1 \right] \left[ s^2 -2s\cos(8\pi/10)+1 \right]$
配列d,e→配列cの多項式の掛け算に関する係数演算を行う関数 poly_times(c,d,e) を実装し、
$\L B(s) = \sum_{k=0}^n c_k s^k$の係数配列を処理しておく。
$\L B(s) = \sum_{k=0}^n c_k s^k$の係数配列を処理しておく。
双一次変換・プリワーピング[2]
sにプリワーピングを考慮した$z^{-1}$の多項式を代入して、フィルタ係数配列a,bに整理する。
$\L s=\left(\frac{2}{T} \frac{1-z^{-1}}{1+z^{-1}} \right) \bigg/ \left( \frac{2}{T} \tan(2\pi f_c T/2) \right) $
$\L s=\left(\frac{2}{T} \frac{1-z^{-1}}{1+z^{-1}} \right) \bigg/ \left( \frac{2}{T} \tan(2\pi f_c T/2) \right) $
ここで、$y=z^{-1}$、および $w= \tan(2\pi f_c T/2) = \tan\left(\frac{\pi}{2} \frac{f_c}{F}\right)$ (Fはナイキスト周波数 1/2T)として
$\L s= \frac{1-y}{w(1+y)}$
$\L s= \frac{1-y}{w(1+y)}$
係数処理
$\L H=\frac{1}{B} = \frac{1}{\sum_{k=0}^n c_k s^k} = \frac{1}{\sum_{k=0}^n c_k/w^k \frac{(1-y)^k}{(1+y)^k}}$
部分分数を処理して、
$\L H= \frac{(1+y)^n}{\sum_{k=0}^n c_k/w^k (1-y)^k (1+y)^{n-k} } $
$a_0=1$とするために、この項を正規化する。
部分分数を処理して、
$\L H= \frac{(1+y)^n}{\sum_{k=0}^n c_k/w^k (1-y)^k (1+y)^{n-k} } $
$a_0=1$とするために、この項を正規化する。
$A_0 = \sum_{k=0}^n c_k/w^k$として、
$\L H= \frac{A_0^{-1} (1+y)^n}{A_0^{-1}\sum_{k=0}^n \frac{c_k}{w^k} (1-y)^k (1+y)^{n-k} } $
$\L H= \frac{A_0^{-1} (1+y)^n}{A_0^{-1}\sum_{k=0}^n \frac{c_k}{w^k} (1-y)^k (1+y)^{n-k} } $
まとめ
$\L b(y) = (1+y)^n \bigg/ \sum_{k=0}^n (c_k/w^k)$
$\L a(y) = \sum_{k=0}^n \frac{c_k}{w^k} (1-y)^k (1+y)^{n-k} \bigg/ \sum_{k=0}^n (c_k/w^k)$
where
$\L w = \tan\left(\frac{\pi}{2} \frac{f_c}{F}\right)$
$\L c_k(s) = B(s)$
$\L a(y) = \sum_{k=0}^n \frac{c_k}{w^k} (1-y)^k (1+y)^{n-k} \bigg/ \sum_{k=0}^n (c_k/w^k)$
where
$\L w = \tan\left(\frac{\pi}{2} \frac{f_c}{F}\right)$
$\L c_k(s) = B(s)$
これを
void poly_times(double *z, double *x, double *y)
void poly_pow(double *z, double *x, int *n)
あたりから実装すれば良い。
void poly_times(double *z, double *x, double *y)
void poly_pow(double *z, double *x, int *n)
あたりから実装すれば良い。
b,aが求まれば、下式からフィルタリングすればよい。
$\L y(t) = \sum_{k=0}^n b_kx(t-k) - \sum_{k=1}^n a_ky(t-k)$
$\L y(t) = \sum_{k=0}^n b_kx(t-k) - \sum_{k=1}^n a_ky(t-k)$
Fortranコード
main.f90
- program main
- use get_butter
- implicit none
-
- integer, parameter :: n = 5
- real(8), parameter :: fc = 0.1d0
- real(8) :: b(n+1), a(n+1)
-
- call butter(b,a,n,fc)
- print *, 'b', real(b)
- print *, 'a', real(a)
-
- end program main
-
get_butter.f90
- module get_butter
- implicit none
-
- contains
- subroutine butter(b,a,n,fc)
- implicit none
- integer, intent(in) :: n
- real(8), intent(out) :: b(n+1), a(n+1)
- real(8), intent(in) :: fc ! = freq_cutoff/freq_Nyquist
-
- integer :: k
- real(8) :: c(n+1), w, A0
- real(8) :: my(n+1), py(n+1)
-
- b(:) = 0d0; a(:) = 0d0;
-
- call butter_poly(c,n)
-
- A0 = 0
- do k=0,n
- A0 = A0 + c(k+1)/w**k
- end do
-
- !b_poly
- b(:) = 0d0; b(1) = 1d0; b(2)=1d0;
- call poly_pow(b, n, n+1)
- b(:) = b(:)/A0
-
- !a_poly
- a(:) = 0d0
- do k=0,n
- my(:) = 0d0; my(1) = 1d0; my(2) = -1d0;
- py(:) = 0d0; py(1) = 1d0; py(2) = 1d0;
- call poly_pow(my, k, n+1)
- call poly_pow(py, n-k, n+1)
- call poly_times(py,my,n+1)
-
- a(:) = a(:) + c(k+1)/w**k * py(:) / A0
- end do
- end subroutine butter
-
- subroutine butter_poly(c,n)
- implicit none
- integer, intent(in) :: n
- real(8), intent(out) :: c(n+1)
-
- integer :: k, kend
- real(8) :: s(n+1)
-
- c(:) = 0d0; c(1) = 1d0;
- if(mod(n,2) == 1) then
- c(2) = 1d0
- kend = (n-1)/2
- else
- kend = n/2
- end if
- do k=1,kend
- s(:) = 0d0
- s(1) = 1d0; s(3) = 1d0;
- call poly_times(c,s,n+1)
- end do
- end subroutine butter_poly
-
- !多項式関数群
- !桁あふれエラーチェックはしないので結果の最大次がmを超えないようにすること
- subroutine poly_times(x,y,m)
- implicit none
- integer, intent(in) :: m !結果の次数
- real(8), intent(inout) :: x(m)
- real(8), intent(in) :: y(m)
-
- integer :: i
- real(8) :: z(m)
- real(8) :: xy(m)
-
- z(:) = 0d0
- do i=1,m
- xy(:) = x(i)*y(:)
- z(i:m) = z(i:m) + xy(1:(m-i+1))
- end do
- x(:) = z(:)
- end subroutine
-
- subroutine poly_pow(x,p,m)
- implicit none
- integer, intent(in) :: p, m !p:べき乗
- real(8), intent(inout) :: x(m)
-
- integer :: i
- real(8) :: z(m)
-
- z(:) = 0d0; z(1) = 1d0
- do i=1,p
- call poly_times(z,x,m)
- end do
- x(:) = z(:)
- end subroutine poly_pow
-
- end module get_butter
-
結果
b 5.97957696E-05 2.98978848E-04 5.97957696E-04 5.97957696E-04 2.98978848E-04 5.97957696E-05 a 1.0000000 -3.9845433 6.4348674 -5.2536159 2.1651332 -0.35992831
Octave結果
octave-3.2.4.exe:39> [b,a]=butter(5,0.1) b = 5.9796e-005 2.9898e-004 5.9796e-004 5.9796e-004 2.9898e-004 5.9796e-005 a = 1.00000 -3.98454 6.43487 -5.25362 2.16513 -0.35993