超音波流体屋のプログラム備忘録

バターワースフィルタ定数

最終更新:

usapfrog

- view
管理者のみ編集可
Matlab / Octave のbutter関数を手実装したかったので。
アルゴリズムのみ記述する。理論的背景は参考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]$

配列d,e→配列cの多項式の掛け算に関する係数演算を行う関数 poly_times(c,d,e) を実装し、
$\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) $

ここで、$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 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$とするために、この項を正規化する。

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

これを
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)$



Fortranコード

main.f90
  1. program main
  2. use get_butter
  3. implicit none
  4.  
  5. integer, parameter :: n = 5
  6. real(8), parameter :: fc = 0.1d0
  7. real(8) :: b(n+1), a(n+1)
  8.  
  9. call butter(b,a,n,fc)
  10. print *, 'b', real(b)
  11. print *, 'a', real(a)
  12.  
  13. end program main
  14.  

get_butter.f90
  1. module get_butter
  2. implicit none
  3.  
  4. real(8), parameter :: PI=3.1415926
  5. contains
  6. subroutine butter(b,a,n,fc)
  7. implicit none
  8. integer, intent(in) :: n
  9. real(8), intent(out) :: b(n+1), a(n+1)
  10. real(8), intent(in) :: fc ! = freq_cutoff/freq_Nyquist
  11.  
  12. integer :: k
  13. real(8) :: c(n+1), w, A0
  14. real(8) :: my(n+1), py(n+1)
  15.  
  16. b(:) = 0d0; a(:) = 0d0;
  17.  
  18. call butter_poly(c,n)
  19. w = tan(PI/2d0*fc)
  20.  
  21. A0 = 0
  22. do k=0,n
  23. A0 = A0 + c(k+1)/w**k
  24. end do
  25.  
  26. !b_poly
  27. b(:) = 0d0; b(1) = 1d0; b(2)=1d0;
  28. call poly_pow(b, n, n+1)
  29. b(:) = b(:)/A0
  30.  
  31. !a_poly
  32. a(:) = 0d0
  33. do k=0,n
  34. my(:) = 0d0; my(1) = 1d0; my(2) = -1d0;
  35. py(:) = 0d0; py(1) = 1d0; py(2) = 1d0;
  36. call poly_pow(my, k, n+1)
  37. call poly_pow(py, n-k, n+1)
  38. call poly_times(py,my,n+1)
  39.  
  40. a(:) = a(:) + c(k+1)/w**k * py(:) / A0
  41. end do
  42. end subroutine butter
  43.  
  44. subroutine butter_poly(c,n)
  45. implicit none
  46. integer, intent(in) :: n
  47. real(8), intent(out) :: c(n+1)
  48.  
  49. integer :: k, kend
  50. real(8) :: s(n+1)
  51.  
  52. c(:) = 0d0; c(1) = 1d0;
  53. if(mod(n,2) == 1) then
  54. c(2) = 1d0
  55. kend = (n-1)/2
  56. else
  57. kend = n/2
  58. end if
  59. do k=1,kend
  60. s(:) = 0d0
  61. s(1) = 1d0; s(3) = 1d0;
  62. s(2) = -2d0 * cos(PI * (2d0*k+n-1d0)/(2d0*n))
  63. call poly_times(c,s,n+1)
  64. end do
  65. end subroutine butter_poly
  66.  
  67. !多項式関数群
  68. !桁あふれエラーチェックはしないので結果の最大次がmを超えないようにすること
  69. subroutine poly_times(x,y,m)
  70. implicit none
  71. integer, intent(in) :: m !結果の次数
  72. real(8), intent(inout) :: x(m)
  73. real(8), intent(in) :: y(m)
  74.  
  75. integer :: i
  76. real(8) :: z(m)
  77. real(8) :: xy(m)
  78.  
  79. z(:) = 0d0
  80. do i=1,m
  81. xy(:) = x(i)*y(:)
  82. z(i:m) = z(i:m) + xy(1:(m-i+1))
  83. end do
  84. x(:) = z(:)
  85. end subroutine
  86.  
  87. subroutine poly_pow(x,p,m)
  88. implicit none
  89. integer, intent(in) :: p, m !p:べき乗
  90. real(8), intent(inout) :: x(m)
  91.  
  92. integer :: i
  93. real(8) :: z(m)
  94.  
  95. z(:) = 0d0; z(1) = 1d0
  96. do i=1,p
  97. call poly_times(z,x,m)
  98. end do
  99. x(:) = z(:)
  100. end subroutine poly_pow
  101.  
  102. end module get_butter
  103.  

結果
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


目安箱バナー