(Library) Butterworthローパスフィルター

作成者: 神代 剛

GPhysオブジェクトの任意の次元にButterworthローパスフィルターを施すインスタンスメソッドです。

フィルターをかける軸方向に欠損値があると使えませんが、その方向にすべて欠損値という場合が含まれているのはOKです(具体的には、陸面がマスクされた海上観測データ時系列、とかを想定しています)。

butterworth.rb (2016-10-18版)

GPhys (電脳Rubyホームページを参照のこと)

上記 butterworth.rb を置いたディレクトリで、

require "./butterworth"

とするとモジュールを読み込めます。

私自身は、~/lib/ruby/filter に butterworth.rb を置き (環境変数 RUBYLIB=$HOME/lib/ruby を設定しています)、

require "filter/butterworth"

として使っています。

メソッドの引数については、ソースのコメントを参照してください。 ソースの後半はサンプルプログラムになっていますので、そちらも参考になると思います。

% ruby butterworth.rb

と上記ファイルそのものを実行すると、サンプルプログラムが実行されます。

※サンプルプログラムの実行には、納多哲史さんの gphys_io_util.rb が必要です。

Butterworthフィルターは再帰型フィルター (recursive filter) なので、実用的には、時系列の最後に0を適当な長さ加えてからフィルターをかける必要があります (いわゆる zero-padding)。 当初はいくつ0を加えるのが適当なのかがよくわからなかったのですが、いろいろ調べてみて、以下のサイトを見つけました:

<URL:http://mailman.ucar.edu/pipermail/ncl-talk/2015-June/003008.html>

これはNCLのユーザーメーリングリスト ncl-talk のアーカイブで、分野も我々と同じですので、ここの記述にしたがってデフォルト値を決めています。

もちろん、陽に指定すれば任意の個数の0を時系列の末尾に加えてフィルターをかけることができます。

# -*- coding: utf-8 -*-
require "numru/gphys/gphys"

module NumRu
  class GPhys

    ###########################################################################
    #
    # Butterworth low-pass filter
    #
    # [in]
    #   dim  (Integer or String) : dimension
    #   q    (Float)             : order
    #   fc   (Float)             : cut-off frequency
    #   npad (Integer)           : number of zeros padded at the end of the
    #                              series (default: 1.5*q/fc)
    # [out]
    #   gp_ret (GPhys) : filtered data
    #
    # References:
    # Hamming, R. W., 1998: Digital Filters, 3rd ed. Dover Publications, Inc.,
    #     284 pp.
    # 伊藤久徳, 見延庄士郎, 2010: 気象学と海洋物理学で用いられるデータ解析法.
    #     気象研究ノート 第221号, 日本気象学会, 253 pp.
    # Roberts, J. and T. D. Roberts, 1978: Use of Butterworth low-pass filter
    #     for oceanographic data. J. Geophys. Res., 83(C11), 5510-5514.
    # [ncl-talk] question regarding Butterworth filter
    #     http://mailman.ucar.edu/pipermail/ncl-talk/2015-June/003008.html
    #
    def butterworth(dim, q, fc, npad=nil)

      npad = (1.5*q/fc).to_i unless npad

      dim = dim_index(dim) if dim.is_a?(String)
      dim += rank if dim < 0
      if dim < 0 || dim >= rank
        raise ArgumentError, "dim #{dim} does not exist"
      end
      if q%2 == 1
        q_is_odd = true
      else
        q_is_odd = false
      end

      x = self.data.val.to_type("float")

      is_nam = false
      if x.class == NArrayMiss
        is_nam = true
        if x.count_invalid == 0
          x = x.to_na
          missing = false
        else
          missing_mask = x.get_mask
          count_dim = missing_mask.to_type("int").sum(dim)
          count_dim_0 = count_dim.eq(0).to_type("int").sum
          count_dim_ndim = count_dim.eq(x.shape[dim]).to_type("int").sum

          if count_dim_0 + count_dim_ndim == x.total/x.shape[dim]
            missing = true
            x = x.to_na(0.0)
          else
            raise "GPhys data include missing values along the filtered dimention."
          end
        end
      end

      x2_shape = x.shape
      x2_shape[dim] += npad
      x2 = NArray.float(*x2_shape).fill(0.0)
      idx = [true]*x.rank
      idx[dim] = 0..(x.shape[dim] - 1)
      x2[*idx] = x

      mm = q/2
      y2 = nil
      # forward
      mm.times do |m|
        y2 = bw_2nd_order(x2, dim, fc, m, q, true)
        x2 = y2
      end
      y2 = bw_1st_order(x2, dim, fc, true) if q_is_odd
      x2 = y2
      # backward
      mm.times do |m|
        y2 = bw_2nd_order(x2, dim, fc, m, q, false)
        x2 = y2
      end
      y2 = bw_1st_order(x2, dim, fc, false) if q_is_odd

      y = y2[*idx].to_type(typecode)
      if is_nam
        if missing
          y = NArrayMiss.to_nam(y,missing_mask)
        else
          y = NArrayMiss.to_nam(y)
        end
      end
      gp_ret = self.copy
      gp_ret.data.replace_val(y)

      return gp_ret
    end

    private

    def bw_1st_order(x, dim, fc, dir_forward)

      x_shape = x.shape
      x_rank = x.rank
      x_shape_dim = x.shape[dim]

      wc = Math::tan(PI*fc)
      c0 = wc/(1 + wc)
      c1 = wc/(1 + wc)
      b1 = (1 - wc)/(1 + wc)

      y = NArray.float(*x_shape).fill(0.0)
      x_shape_dim.times do |i|
        t = [true]*x_rank
        if dir_forward
          t[dim] = i
        else # backward
          t[dim] = (x_shape_dim - 1) - i
        end
        if i == 0
          y[*t] = c0*x[*t]
        else
          t_1 = [true]*x_rank
          if dir_forward
            t_1[dim] = i - 1
          else # backward
            t_1[dim] = (x_shape_dim - 1) - (i - 1)
          end
          y[*t] = c0*x[*t] + c1*x[*t_1] + b1*y[*t_1]
        end
      end

      return y
    end

    def bw_2nd_order(x, dim, fc, k, n, dir_forward)

      x_shape = x.shape
      x_rank = x.rank
      x_shape_dim = x.shape[dim]

      wc = Math::tan(PI*fc)
      b0   = wc**2 + 2*wc*Math::sin(PI*(2*k+1)/(2*n)) + 1
      b1b0 = 2 - 2*wc**2
      b2b0 = -(wc**2 - 2*wc*Math::sin(PI*(2*k+1)/(2*n)) + 1)
      c0b0 = wc**2
      c1b0 = 2*wc**2
      c2b0 = wc**2

      b1 = b1b0/b0
      b2 = b2b0/b0
      c0 = c0b0/b0
      c1 = c1b0/b0
      c2 = c2b0/b0

      y = NArray.float(*x_shape).fill(0.0)
      x_shape_dim.times do |i|
        t = [true]*x_rank
        if dir_forward
          t[dim] = i
        else # backward
          t[dim] = (x_shape_dim - 1) - i
        end
        if i == 0
          y[*t] = c0*x[*t]
        else
          t_1 = [true]*x_rank
          if dir_forward
            t_1[dim] = i - 1
          else # backward
            t_1[dim] = (x_shape_dim - 1) - (i - 1)
          end
          if i == 1
            y[*t] = c0*x[*t] + c1*x[*t_1] + b1*y[*t_1]
          else
            t_2 = [true]*x_rank
            if dir_forward
              t_2[dim] = i - 2
            else # backward
              t_2[dim] = (x_shape_dim - 1) - (i - 2)
            end
            y[*t] = c0*x[*t] + c1*x[*t_1] + c2*x[*t_2] + b1*y[*t_1] + b2*y[*t_2]
          end
        end
      end

      return y
    end

  end
end

if $0 == __FILE__
##################################################
# 気象研究ノート No. 221
# 「気象学と海洋物理学で用いられるデータ解析法」
# 第7章 時間フィルタリング を参考に作図
##################################################
require "numru/ggraph"
require "gphys_io_util"
include NumRu

iws = (ARGV[0] ? ARGV[0].to_i : 4)
DCL.gropn(iws)
DCL.sldiv("y", 2, 1)
DCL.sgpset("lclip", true)

### 図7.8 (p.108) ################################
# 元の図はハイパスフィルターの例なので、ローパス
# フィルターの場合は「フィルター時系列」と「元時
# 系列との差」が入れ替わった結果になることに注意。

na = NArray.float(257).fill(0.0)
(124..128).each do |t|
  na[t] = 2.0*(t.to_f - 123.0)/5.0
end
(129..133).each do |t|
  na[t] = 2.0*(133 - t.to_f)/5.0
end
gp = na2gphys(na,[NArray.float(na.size).indgen])
GGraph.next_fig("viewport"=>[0.15,0.85,0.15,0.85],"window"=>[0.0,256.0,-0.8,2.2])
GGraph.next_axes("xunits"=>"","xtitle"=>"","yunits"=>"","ytitle"=>"")
GGraph.line(gp,true,"title"=>"")
gp2 = gp.butterworth(0, 8, 1.0/30.0)
GGraph.line(gp2,false,"index"=>3)
gp3 = gp - gp2
GGraph.line(gp3,false,"type"=>3)
DCL.uulinz([170,200], [1.8,1.8], 1, 1)
DCL.uulinz([170,200], [1.6,1.6], 1, 3)
DCL.uulinz([170,200], [1.4,1.4], 3, 1)
DCL.sgtxzu(210.0, 1.8, "original", 0.02, 0, -1, 3)
DCL.sgtxzu(210.0, 1.6, "Butter",   0.02, 0, -1, 3)
DCL.sgtxzu(210.0, 1.4, "diff.",    0.02, 0, -1, 3)

### 次数を高くするとリンギングが強く出る #########
# これは図にはないが、本文に説明あり。

GGraph.next_fig("viewport"=>[0.15,0.85,0.15,0.85],"window"=>[0.0,256.0,-0.2,0.9])
GGraph.next_axes("xunits"=>"","xtitle"=>"","yunits"=>"","ytitle"=>"")
GGraph.line(gp,true,"title"=>"")

gp2 = gp.butterworth(0, 2, 1.0/30.0)
GGraph.line(gp2,false,"index"=>23)
gp4 = gp.butterworth(0, 4, 1.0/30.0)
GGraph.line(gp4,false,"index"=>33)
gp8 = gp.butterworth(0, 8, 1.0/30.0)
GGraph.line(gp8,false,"index"=>43)
DCL.sgtxzu(240.0, 0.77, "2nd", 0.02, 0, 1, 3)
DCL.sgtxzu(240.0, 0.70, "4th", 0.02, 0, 1, 3)
DCL.sgtxzu(240.0, 0.63, "8th", 0.02, 0, 1, 3)
DCL.uulinz([180,210], [0.77,0.77], 1, 23)
DCL.uulinz([180,210], [0.70,0.70], 1, 33)
DCL.uulinz([180,210], [0.63,0.63], 1, 43)

gp1 = gp.butterworth(0, 1, 1.0/30.0)
GGraph.line(gp1,false,"type"=>2,"index"=>1)
gp3 = gp.butterworth(0, 3, 1.0/30.0)
GGraph.line(gp3,false,"type"=>3,"index"=>1)
DCL.sgtxzu(240.0, 0.51, "1st", 0.02, 0, 1, 3)
DCL.sgtxzu(240.0, 0.44, "3rd", 0.02, 0, 1, 3)
DCL.uulinz([180,210], [0.51,0.51], 2, 1)
DCL.uulinz([180,210], [0.44,0.44], 3, 1)

### 図7.9 (p.109) ################################
# 元の図左の破線は、キャプションには陽に書かれて
# いないが、オリジナルデータにそのままフィルター
# をかけた結果ではなく、データの最後にゼロを30個
# 加えてフィルターをかけた結果だと思われる(その
# ように作図すると、元の図と一致した)。本文とは
# 齟齬があるのだが、左の破線は両端も含めて低周波
# 成分をうまく取り出せているように見える(つまり、
# 右の実線のほうに対応している)。

t = NArray.float(150).indgen
f  = 2.0*sin(2.0*PI*t/12.0) + 1.5*cos(2.0*PI*t/5.5)
f2 = 2.0*sin(2.0*PI*t/12.0)

# 左図
gp = na2gphys(f,[t])
GGraph.next_fig("viewport"=>[0.15,0.85,0.15,0.85],"window"=>[0.0,149.0,-4.0,4.0])
GGraph.next_axes("xtickint"=>5,"xlabelint"=>20,"xunits"=>"day","xtitle"=>"Time","ytickint"=>0.2,"ylabelint"=>1,"yunits"=>"","ytitle"=>"Amplitude")
GGraph.line(gp, true, "title"=>"")
gp_lp = gp.butterworth(0, 8, 1.0/8.0, 30)
GGraph.line(gp_lp, false, "type"=>2)

# 右図
gp2 = na2gphys(f2,[t])
GGraph.next_fig("viewport"=>[0.15,0.85,0.15,0.85],"window"=>[0.0,149.0,-1.0,1.0])
GGraph.next_axes("xtickint"=>5,"xlabelint"=>20,"xunits"=>"day","xtitle"=>"Time","ytickint"=>0.1,"ylabelint"=>0.2,"yunits"=>"","ytitle"=>"")
GGraph.line(gp2 - gp_lp, true, "title"=>"")
gp_lp2 = gp.butterworth(0, 8, 1.0/8.0, 0)
GGraph.line(gp2 - gp_lp2, false, "type"=>2)

DCL.grcls
end

実装にあたっては、西本絵梨子さんの (Library) GPhys（加重）移動平均メソッド および (Library) Lanczosフィルター が大変参考になりました。ありがとうございます。

Hamming, R. W., 1998: Digital Filters, 3rd ed. Dover Publications, Inc., 284 pp.

伊藤久徳, 見延庄士郎, 2010: 気象学と海洋物理学で用いられるデータ解析法. 気象研究ノート 第221号, 日本気象学会, 253 pp.

Roberts, J. and T. D. Roberts, 1978: Use of Butterworth low-pass filter for oceanographic data. J. Geophys. Res., 83(C11), 5510-5514.

[ncl-talk] question regarding Butterworth filter <URL:http://mailman.ucar.edu/pipermail/ncl-talk/2015-June/003008.html>