Numo::NArrayをpocketfftでフーリエ変換するGemを作成した

はじめに

NumPyが、ver. 1.17から、フーリエ変換にFFTPACKではなくpocketfftを使う様になった。pocketfftはFFTPACKを修正したもので、速度と精度が改善されている。

ENH: Add pocketfft sources to numpy for testing, benchmarks, etc. by mreineck · Pull Request #11888 · numpy/numpy · GitHub

https://gitlab.mpcdf.mpg.de/mtr/pocketfft

pocketfftは、単一のC言語のソースファイルからなるので「これはソースを同梱する形でNumo::NArrayな配列をフーリエ変換するなにかが作れるぞ」と思い、早速作ってみた。Numo::NArrayでは、フーリエ変換のAPIは提供されておらず、別でNumo::FFTWやNumo::FFTEをインストールする形になっている。そこで、Numo::Pocketfftとして公開した。

numo-pocketfft | RubyGems.org | your community gem host

使い方

pocketfftを同梱しているので、普通にgemコマンドでインストールできる。Mac、Ubuntu、Windows 10（RubyInstallerな環境）でインストールできることを確認した。

$ gem install numo-pocketfft

メソッド名や動きはnumpy.fftに従った。普通の高速フーリエ変換をfftメソッド、ifftメソッドで、実数な配列を想定した高速フーリエ変換をrfftメソッド、irfftメソッドで行う。

Module: Numo::Pocketfft — Documentation for numo-pocketfft (0.1.0)

まずは、実数の1次元配列のフーリエ変換・逆変換は以下のようになる。

irb(main):001:0> require 'numo/pocketfft'
=> true
irb(main):002:0> a = Numo::DFloat.new(4).rand
=> Numo::DFloat#shape=[4]
[0.0617545, 0.373067, 0.794815, 0.201042]
irb(main):003:0> x = Numo::Pocketfft.rfft(a)
=> Numo::DComplex#shape=[3]
[1.43068+0i, -0.733061-0.172025i, 0.282461+0i]
irb(main):004:0> Numo::Pocketfft.irfft(x)
=> Numo::DFloat#shape=[4]
[0.0617545, 0.373067, 0.794815, 0.201042]
irb(main):005:0>

逆変換の結果が複素数で返ってきても良いなら、全てfft系メソッドでもよい。例えば実数の2次元配列のフーリエ変換・逆変換は以下のようになる。fft2と2をつける（これはnumpy.fftにあわせた）。

irb(main):001:0> require 'numo/pocketfft'
=> true
irb(main):002:0> a = Numo::DFloat.new(4, 4).rand
=> Numo::DFloat#shape=[4,4]
[[0.0617545, 0.373067, 0.794815, 0.201042],
 [0.116041, 0.344032, 0.539948, 0.737815],
 [0.165089, 0.0508827, 0.108065, 0.0687079],
 [0.904121, 0.478644, 0.342969, 0.164541]]
irb(main):003:0> Numo::Pocketfft.ifft2(Numo::Pocketfft.fft2(a))
=> Numo::DComplex#shape=[4,4]
[[0.0617545+0i, 0.373067+0i, 0.794815+0i, 0.201042+0i],
 [0.116041+0i, 0.344032+0i, 0.539948+0i, 0.737815+0i],
 [0.165089+0i, 0.0508827+0i, 0.108065+0i, 0.0687079+0i],
 [0.904121+0i, 0.478644+0i, 0.342969+0i, 0.164541+0i]]

任意の次元数を受け付けるfftnとrfftnメソッドもある。

irb(main):001:0> require 'numo/pocketfft'
=> true
irb(main):002:0> a = Numo::DFloat.new(2,2,2).rand + Complex::I * Numo::DFloat.new(2,2,2).rand
=> Numo::DComplex#shape=[2,2,2]
[[[0.0617545+0.165089i, 0.373067+0.0508827i],
  [0.794815+0.108065i, 0.201042+0.0687079i]],
 [[0.116041+0.904121i, 0.344032+0.478644i],
  [0.539948+0.342969i, 0.737815+0.164541i]]]
irb(main):003:0> Numo::Pocketfft.ifftn(Numo::Pocketfft.fftn(a))
=> Numo::DComplex#shape=[2,2,2]
[[[0.0617545+0.165089i, 0.373067+0.0508827i],
  [0.794815+0.108065i, 0.201042+0.0687079i]],
 [[0.116041+0.904121i, 0.344032+0.478644i],
  [0.539948+0.342969i, 0.737815+0.164541i]]]

といった感じで、フーリエ変換・逆変換できる。

git submodule による外部ライブラリの同梱

pocketfftの同梱は、Gitlab上のpocketfftのリポジトリを、submoduleとして取り込むことで実現した。submoduleの追加は、特別な工夫はなくgitコマンド通りになる。

$ cd ext/numo/pocketfft
$ git submodule add https://gitlab.mpcdf.mpg.de/mtr/pocketfft pocketfft

Ruby extensionでは、mkmfというライブラリでMakefileを作成する。mkmfでは、基本的には、全てのコードが同一のディレクトリ下（Numo::Pocketfftを例にすればext/numo/pocketfft下）にあることを想定している。submoduleで追加したpocketfftは、サブディレクトリにある。この階層構造を、mkmfに教えるようなことが必要になる。

# extconf.rb
# ...省略

# ext/numo/pocketfft下のコードを$srcsに追加する
$srcs = Dir.glob("#{$srcdir}/*.c").map { |path| File.basename(path) }

# pocketfftのコードを$srcsに追加する
$srcs << 'pocketfft.c'

# ext/numo/pocketfft/pocketfftをインクルードパスやmakeの探索パスに追加する
# ※ インクルードパスはいらないかも
Dir.glob("#{$srcdir}/*/") do |path|
  dir = File.basename(path)
  $INCFLAGS << " -I$(srcdir)/#{dir}"
  $VPATH << "$(srcdir)/#{dir}"
end

# ...省略

これで、rake compileでpocketfftも含めてコンパイルできる。そして、gemファイルを作る際に、submoduleなpocketfftを含めるために、gemspecファイルで、Gem::Specificationのfilesにsubmoduleへのパスを追加する。

# numo-pocketfft.gemspec
# ...省略

Gem::Specification.new do |spec|
# ...省略

  # submoduleのディレクトリ下のファイルをfilesに追加する
  # 参考: https://gist.github.com/mattconnolly/5875987
  gem_dir = __dir__ + '/'
  `git submodule --quiet foreach pwd`.split($OUTPUT_RECORD_SEPARATOR).each do |submodule_path|
    Dir.chdir(submodule_path) do
      submodule_relative_path = submodule_path.sub gem_dir, ''
      `git ls-files`.split($OUTPUT_RECORD_SEPARATOR).each do |filename|
        spec.files << "#{submodule_relative_path}/#{filename}"
      end
    end
  end

# ...省略

これでsubmoduleにより外部ライブリをgemに同梱できた。後で発見したが、同様のことがruby-eigenでも行われていた。ruby-eigenでは、C++の線形代数ライブラリEigenが、submoduleにより同梱されている。

おわりに

画像にフィルタを掛ける場合、画像とフィルタとの畳み込み演算なる。畳み込み演算は、フーリエ変換により複素数の世界に持っていくと、シンプルな乗算になる。これを利用すると、畳み込み演算をフーリエ変換・逆変換で実装できる。画像処理ライブラリのMagroの開発で、このフーリエ変換による畳み込みを必要としていたので、良い機会だと思いpocketfftによるフーリエ変換のgemを作った。NumPyのフーリエ変換では、正規化のオプションがあったりするが、Numo::Pocketfftでは、まずは、シンプルなフーリエ変換・逆変換を実装した。今後bugfixとかも含めて、バージョンアップしていくつもりでいる。

github.com