(Library) Lanczosフィルター

作成者: 西本絵梨子

Lanczosフィルターに関するライブラリです。 時空間領域に対応した重み関数、周波数空間に対応した応答関数、を返します。

Lanczos.rb ライブラリ（2013/11/08版）。

github Lanczos

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

上記ファイルを同一ディレクトリもしくは、パスの通っているディレクトリに格納し、

require 'Lanczos'

でモジュールを読み込んで下さい。

Tips

データに時空間フィルターをかけるには、 時空間領域でのフィルタリング と 周波数空間でのフィルタリング の二通りがある。 詳細は、Duchon(1979)、伊藤・見延(2010)などを参照してください。

• 時空間領域でのフィルタリングの仕方

重み関数（Lanczos::weight_fなど）を計算し、時空間データに対して加重移動平均 （GPhys加重移動平均メソッド）を施す。

• 周波数空間でのフィルタリングの仕方

まず、フーリエ変換（GPhysのFFTメソッド）などによって、時空間データを周波数領域データに変換する。次にその周波数領域データに、応答関数（Lanczos::ref_fなど）をかけてフィルターをかける。最後にフーリエ逆変換などによって、時空間データに戻す。

主な関数は次のようになっています。

Lanczos::weight_f( fc , n ) =>（ローパスの）重み関数を返します

• fc : カットオフ周波数
• n : 重み

Lanczos::weight_f＿bp( fc1, fc2 , n ) => バンドパスの重み関数を返します

• fc1: カットイン周波数
• fc2: カットオフ周波数
• n : 重み

Lanczos::res_f( wgt, f ) => 周波数空間での応答関数を返します

• wgt: 重み関数
• f : 連続した周波数

Lanczos::bp_filter( f, fc1, fc2, n )

Lanczos::low_filter( f, fc2, n )

Lanczos::high_filter( f, fc1, n )

それぞれ、周波数空間での応答関数を返します。 中で、Lanczos::weight_f と Lanczos::res_f を呼んでいます。

# -*- coding: utf-8 -*-

require 'narray'
include NMath
module Lanczos

  #############
  # NArray版 虚数単位 i (長さ1の配列)
  I = NArray.complex(1).fill( Complex(0.0,1.0) )

  module_function

  # sinc function
  #
  # 入力
  #  - x (NArray)
  # 出力
  #   - sinc (NArray)
  def sinc(x)
    n = x.length
    n_2 = (n-1)/2

    sinc=NMath::sin(PI*x)/(PI*x)
    sinc[n_2]=1 if x[n_2]==0
    return sinc
  end

  # Lanczos weight function
  #
  # 入力
  #   - fc (Float) : cut-off(or cut-in) frequency
  #   - n (Integer) : number of weights
  # 出力
  #   - w (NArray) : weight function
  def weight_f(fc,n)
    n = n.to_i
    if n%2 == 0
      raise "weights number need to be odd number!\n"
    end
    n_2 = (n-1)/2
    k = NArray.float(n).indgen!-n_2

    sigma = sinc(k/n_2)
    w = 2.0 * fc * sinc(2.0*fc*k) * sigma

    return w
  end

  # Lanczos band-pass weight function
  #
  # 入力
  #   - fc1 (Float) : cut-in frequency
  #   - fc2 (Float) : cut-off frequency
  #   - n (Integer) : number of weights
  # 出力
  #   - w (NArray) : weight function
  def weight_f_bp(fc1,fc2,n)
    n = n.to_i
    if n%2 == 0
      raise "weights number need to be odd number!\n"
    end
    n_2 = (n-1)/2
    if n_2 < (1.3/(fc2-fc1)).ceil
      raise "weights number need to be more than 2*1.3/(fc2-fc1)+1\n"
    end

    k = NArray.float(n).indgen!-n_2
    sigma = sinc(k/n_2)
    w = 2.0 * ( fc2*sinc(2.0*fc2*k) - fc1*sinc(2.0*fc1*k) ) * sigma

    return w
  end

  # Lanczos Response function
  #
  # 入力
  #   - w (NArray) : weight function
  #   - f (NArray) : sequence of frequencies
  # 出力
  #   - r (NArray) : response function
  def res_f(w,f)
    n = w.length
    n_2 = (n-1)/2
    k = NArray.complex(1,n).indgen!-n_2 
    w = w.reshape(1,n).to_type(NArray::COMPLEX)
    r = (w*NMath::exp(I*2*PI*f*k)).sum(-1)
    #p "r"
    #p r

    return r
  end

#########

  # Lanczos Band-Pass Response Function
  #
  # 入力
  #   - f (NArray) : sequence of frequencies
  #   - fc1 (Float) : cut-in frequency
  #   - fc2 (Float) : cut-off frequency
  #   - n (Integer) : number of weights
  # 出力
  #   - r (NArray) : band-pass response function
  def bp_filter( f, fc1, fc2, n )
    w = Lanczos::weight_f_bp( fc1, fc2, n )
    r = Lanczos::res_f( w, f )
    return r
  end

  # Lanczos Low-Pass Response Function
  #
  # 入力
  #   - f (NArray) : sequence of frequencies
  #   - fc2 (Float) : cut-off frequency
  #   - n (Integer) : number of weights
  # 出力
  #   - r (NArray) : low-pass response function
  def low_filter( f, fc2, n )
    w = Lanczos::weight_f( fc2, n )
    r = Lanczos::res_f( w, f )
    return r
  end

  # Lanczos High-Pass Response Function
  #
  # 入力
  #   - f (NArray) : sequence of frequencies
  #   - fc1 (Float) : cut-in frequency
  #   - n (Integer) : number of weights
  # 出力
  #   - r (NArray) : high-pass response function
  def high_filter( f, fc1, n )
    w = Lanczos::weight_f( fc1, n )
    r = 1-Lanczos::res_f( w, f )
    return r
  end

#########

end


##########################################################
# メイン
if $0 == __FILE__
  require 'numru/dcl'
  include NumRu

  #< 引数解釈 >
  fc1 = ARGV[0] ? ARGV[0].to_f : 0.2
  fc2 = ARGV[1] ? ARGV[1].to_f : 0.3
  n   = ARGV[2] ? ARGV[2].to_i : 43

  #< f >
  m = 365*2
  f = NArray.sfloat(m).indgen/m.to_f
  fN = 1.0/2.0  # Nyquist frequency

  DCL.gropn 1
  DCL.grfrm

  DCL.grsvpt 0.2,0.8,0.2,0.8
  DCL.grswnd 0.0,1.0,-0.2,1.1
  DCL.grstrn 1
  DCL.grstrf

  DCL.usdaxs

  # idealized filter response
  DCL.uulinz [0.0,fc1],[0.0,0.0],1,23
  DCL.uulinz [fc1,fc2],[1.0,1.0],1,23
  DCL.uulinz [fc2,fN],[0.0,0.0],1,23
  DCL.uulinz [fc1,fc1],[0.0,1.0],3,21
  DCL.uulinz [fc2,fc2],[0.0,1.0],3,21

  r = Lanczos::bp_filter( f, fc1, fc2, n )
  r_low = Lanczos::low_filter( f, fc2, n )
  r_high = Lanczos::high_filter( f, fc1, n )

  # Lanczos filter
  DCL.uulinz f,r,1,3

  DCL.uxsttl 't',"n=#{n}",-1
#}

  DCL.grcls
  
end

Duchon, C. E. (1979). Lanczos filtering in one and two dimensions. Journal of Applied Meteorology, 18, 1016-1022.

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