Pythonでフーリエ変換をRubyでやってみる

## はじめに

[Pythonでフーリエ変換](https://qiita.com/kaityo256/items/64a54bb2e2c477cc6fa1)という記事をみた。

数学は全くわからないのだが、Ruby向けのグラフ描出ツール[GR.rb](https://qiita.com/kojix2/items/456259f40d5563005cd1)を作っている人間としてはRubyで同じことができるか気になってしまう。

そこで最初のガウス分布のフーリエ変換を見様見真似で真似してみた。数学は全くわからないのだが。
グラフ描出には[GR.rb](https://qiita.com/kojix2/items/456259f40d5563005cd1)を使う。

## ガウス分布のフーリエ変換

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/07300a9d-0a5f-38ee-5437-6f4703aff7d2.png)


```ruby
require 'numo/narray'
require 'numo/pocketfft'
require 'gr/plot'

N = 4096     # サンプル数
s = N / 256  # 標準偏差
pi = Math::PI

y = Array.new(N) do |i|
  x = i - N / 2
  Math.exp(-x**2 / (2.0 * s**2)) / (Math.sqrt(2 * pi) * s)
end

freq = Numo::DFloat.linspace(0,1,N+1)
freq[N/2..-1].inplace - 1.0
freq = freq[0..-2]

theory = freq.map { |k| Math.exp(-(2 * pi * k)**2 * s**2 / 2.0) }

x = freq
y = Numo::Pocketfft.fft(Numo::NArray.cast(y)).abs

GR.plot([x, y, 'bx'], [x, theory, 'r'], xlim: [-0.05, 0.05])
GR.savefig('gaussian.png')
```

FFTを実行してくれそうなライブラリはRubyにも結構いろいろある。numo-fftwとか、numo-ffteとか。

今回は[Rumale](https://github.com/yoshoku/rumale)の作者である@yoshokuさんが作成している新しい[numo-pocketfft](https://github.com/yoshoku/numo-pocketfft)を利用した。

途中でつまずいたのは、`fftfreq` 関数がないことである。

最初は[PyCall](https://github.com/mrkn/pycall.rb)を使ってPythonを呼んでいたが

```ruby
require 'pycall/import'
include PyCall::Import
pyimport :numpy, as: :np
freq = np.fft.fftfreq(N) # Python配列
freq = freq.tolist.to_a # Rubyの配列に変換
```

結局上のサンプルでは、NArrayの機能を使って自分で作るようにした。
NArrayの `inplace`は特に深い意味はなくて `n.inplace + b` は `n+=b` と同じです。

数学はわからないですが、Rubyでも似たようなグラフが描けていると思う。

## おまけ: uplotの紹介

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/9b15ae10-21bc-971b-7d4f-386775408a71.png)

現在、標準入力から[UnicodePlot](https://github.com/red-data-tools/unicode_plot.rb)を使ってグラフを描出する[uplot](https://github.com/kojix2/uplot)というコマンドラインツールを作成中である。
これを使うとこんなこともできる。(元のガウス分布のグラフ)

```ruby:g.rb
N = 4096     # サンプル数
s = N / 256  # 標準偏差
pi = Math::PI

N.times do |i|
  x = i - N / 2
  puts Math.exp(-x**2 / (2.0 * s**2)) / (Math.sqrt(2 * pi) * s)
end
```

```sh
ruby g.rb | uplot line --xlim 1948..2148
```

完成したらQiitaに紹介記事を書きます。

## 追記

なぜかトレンドに上がってしまっているので慌てて他のものも追加。

#### 時系列（指数減衰）

```ruby
require 'numo/narray'
require 'numo/pocketfft'
require 'gr/plot'

N = 4096
h = 0.1
gamma = 0.2
v0 = 2.0

v = v0
f = []
x = []

N.times do |i|
  v -= gamma * v * h
  f << v
  x << i
end

y = Numo::Pocketfft.fft(Numo::NArray.cast(f))

freq = Numo::DFloat.linspace(0, 1, N + 1) * 10
freq[N / 2..-1].inplace - 10.0
freq = freq[0..-2]

GR.plot(freq, y.abs, 'g.', grid: false)
GR.savefig('decay.png')
```

freq検算用

```ruby
# require 'pycall/import'
# include PyCall::Import
# pyimport :numpy, as: :np
# freq2 = np.fft.fftfreq(N, d: 0.1) # Python配列
# freq2 = freq2.tolist.to_a         # Rubyの配列に変換
# p freq2 == freq                   # true
```

![decay.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/c8796c22-abd6-3603-3dd1-56937553adac.png)

```ruby
theory = freq.map { |w| v0 / Math.sqrt((w * 2 * Math::PI)**2 + gamma**2) / h }
GR.plot([freq, y.abs, 'g.'], [freq, theory, 'r.'],
        grid: false, xlog: true, ylog: true)
```

![decay.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/1d23a94c-907e-9bc1-6961-8a643f8a2f2f.png)

#### 時系列(ランジュバン方程式)

```ruby
require 'numo/narray'
require 'numo/pocketfft'
require 'gr/plot'

N = 4096
h = 0.1          # 時間刻み
gamma = 0.2      # 減衰係数
v0 = 0.0         # 初速

v = v0
x = []
f = []


N.times do |i|
  v = v -gamma * v * h + Numo::DFloat.new(1).rand_norm(0, Math.sqrt(h))[0]
  x << i
  f << v
end

GR.plot(x,f, xlim: 0..4096)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/6c373d00-3182-254b-2a01-d8dabe8c27a3.png)

```ruby
fw = Numo::Pocketfft.fft(Numo::NArray.cast(f))

freq = Numo::DFloat.linspace(0, 1, N + 1) * 10
freq[N / 2..-1].inplace - 10.0
freq = freq[0..-2]

GR.plot(freq, fw.abs, '.')
```

![fig.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/a0048aa1-50ca-420b-63a6-44e8335d0804.png)

```ruby
theory = freq.map{|w| Math.sqrt(N*h)/Math.sqrt((2*Math::PI*w)**2 + gamma**2)/h}

GR.plot([freq, fw.abs, '.'], [freq, theory, 'r'], 
         size: [100, 75], xlog: true, ylog: true)
```

![fig.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/0ee1d422-cc63-c7a7-88a6-362679b0b91a.png)

理論面は[Pythonの元記事](https://qiita.com/kaityo256/items/64a54bb2e2c477cc6fa1)を見てください。グラフはなんとなく正しそうだけど、理論面は私は全くわかってないので、一応確認してから使ってください。

この記事は以上です。