Houdiniでフーリエ変換を可視化してみる

はじめに

はじめまして。
最近、フーリエ変換についてHoudini Indieを活用して学習する機会がありました。
……のですが、数学オンチの自分としてはなかなかに大変でした。。。

本記事では、私同様に「数学は苦手だけど業務の関係でフーリエ変換を学びたい！」という人がいるかもしれない！という思いから、私がどのような過程で学習し、成果物を出せるに至ったかをまとめていきたいと思います。

最終的な成果物

2次元離散フーリエ変換を用い画像をスペクトル化した後、cos波の重ね合わせで元画像を復元しています。一般的な画像処理におけるフーリエ変換同様、高周波カットによって画像をぼかしたりなども可能になっています。

おことわり

• フーリエ変換がどういったものかという説明は記載しません。
  学習過程とHoudiniでの実装方法に焦点を当てています。
• ここではフーリエ変換を自分のスキルの範疇で実装することで、自身の「引き出し」を増やすことを目標としています。故に数学的に正確であることを必ずしも重視せず、「とりあえずそれらしいものができる」状態に達することを第一として学習しました。
  • 未解決の疑問があっても先に進むことを優先している場合があります。
• 上記に加え、筆者の数学知識は決して深くないので、おかしな部分があるかもしれません。
  • 用語用法がおかしいなど。よろしければご指摘いただけると嬉しいです。
  • 高校数学が怪しいくらいのレベルを想定してください。。。

0. 基本戦略

目標としたのは先にお見せした、画像処理的なデモをHoudiniを用いて作ることですが、いきなり2次元バージョンをつくるのは無理だろうとわかっていたので、まずは1次元バージョンから挑戦してみます。

有名だと思うのですが、3Blue1BrownJapanによるこちらの動画の考え方をベースにしました。

詳しくは実際に視聴していただきたいのですが、動画ではフーリエ変換は「波を原点回りにぐるぐる巻きつけることで周波数を分解するマシン」と紹介されています。

Houdiniで可視化するにあたり、ポイントを繋ぐことで波をカーブとして描くのは問題ないでしょう。
そしてそれを原点に巻き付けるように変形させることも、やればできそうな気がします。

なんとなくはじめの一歩が見えた気がしますね。
それでは、次の項から実際に「巻き付けマシン」をHoudiniで作っていってみましょう。

1. 実践編 1次元フーリエ変換

数式を解釈する

先に紹介した3Blue1BrownJapanの動画中では、フーリエ変換の式は以下のような形だと述べられています。

\hat{g}(f) = \int_{-\infty}^{+\infty}g(t)e^{-2{\pi}ift}dt

動画の内容に沿って式中の各要素の意味を見ていくと……

• $e^{-2{\pi}ift}$
  • $e^{-2{\pi}i}$ ・・・1倍で時計回りに1周になるように（オイラーの公式）
  • $^f$ ・・・ 巻き付け周波数
  • $^t$ ・・・ 時間（元の波のxの値）によって回転する
• $g(t)$・・・ 元の波(信号の強さと時間の関数)。 $e^{-2{\pi}ift}$　に掛け算することで原点回りに巻き付けることになる
• $\int_{-\infty}^{+\infty}$および$dt$ ・・・積分でグラフの重心をとる。ただし時間の区間で割っていないので原点から重心位置はスケールアップされる

……という風になっています。

この操作を元の波を構成する各点に適用してやれば、1回分の巻き付けがうまくいくというわけですね。
虚数については無視してしまいます。重心が求まるならば複素平面上でなくてもよいと判断しました。

そして動画によると、最終的に求めたいのは、巻き付け周波数を段々と大きくしていった時の、重心の位置の変化をプロットしたグラフであるようです。
となると、forEachLoopを使って、周波数を変えながら巻き付けを繰り返していってあげればよさそうです。

どのようにネットワークを組めばいいか、徐々に想像できてきました。

SOPネットワークを組む

実際に組んでみます。
ひとまず図のような感じになりました。

最初に十分なポイント数を持ったLineを作り、AttributeWrangleSOPで波形にします。

その後できた波形をforEachLoop内に投入します。
1回のループの中では、

1. 巻き付け
2. 重心を求めポイントを作る
3. 重心ポイントに現在の巻き付け周波数を持たせる

ということをしています。
巻き付けの部分は実に単純で、sinとcosで入力波形のポイントごとのx位置（=数式の「時間t」に相当）に対応する向きを求め（オイラーの公式の部分）、波の高さをベクトルの長さとしてスケールしてあげるだけです。

重心を求めるところは、そのままズバリなExtractCentroidSOPがあるのでそれで代用してしまいました。
その代わりループの外で、巻き付ける波形の区間の長さ（今回の場合は最初のLineの長さ）を出力された重心ポイント一つ一つに対し掛け算することで、重心を求める過程で行われる割り算を打ち消しています。

あとはループ内で付与しておいた巻き付け周波数をx軸、重心位置をY軸としたグラフになるよう再プロットしてあげます。

この時、重心位置はx成分のみ考慮します。y成分は本来複素数の虚部だから要らないんじゃないか、という判断です。

ここまでの内容を試してみると、無事最初に入力した周波数の部分だけがスパイクとなるようなスペクトルとして現れてくれました。（ただし、振幅の大きさは再現できていません）

次はここで得られたスペクトルから、元の波を復元できるか試してみます。

逆変換を試してみる

残念ながら3Blue1BrownJapanの動画では逆フーリエ変換（以後逆変換と呼びます）までは触れられていませんでした。
仕方がないのでWikipediaに乗っている逆変換の数式を見て考えていくことにします。

まず確認ですが、Wikipediaだとフーリエ変換の方の記号が動画と異なっていますが、違う記号を使っているだけで数式としては同じのようです。
そして逆変換のほうをみると……

f(x):=\int_{-\infty}^{\infty}\hat{f}(\xi)e^{2{\pi}ix\xi}d\xi

なんか普通のフーリエ変換の式とすごい似てますね。。。
見た感じオイラーの公式も重心を求める積分も同じように存在しています。
違うのは元の波形を掛け算していたところが先ほど求めたスペクトルのグラフに置き換わっている点ですね。

ということは、スペクトルをフーリエ変換同様原点に巻き付けていけば、元の波形が得られるのでは…？

試してみましょう。

やることはフーリエ変換の方とほとんど同じにします。

スペクトルを入力として、ループの中で巻き付け、重心計算をして、重心の変化をグラフとしてプロットしなおします。
ただしスペクトルのY軸が大きすぎるとうまくいかなそうだったので、「multiplyTimeRange」の掛け算を打ち消してみています。
また、2回目の「multiplyTimeRange」で掛ける値はスペクトルの長さに変更しておきます。

結果は……

おお、それっぽい！

……あれ、でも徐々に振幅が弱くなっていってしまいますね。。。
これは元の波形には見られなかった現象です。

この現象についてはいろいろと考えたのですが、結局未解決になっています。
ただ、そもそも コンピューター上で計算しているのだから離散フーリエ変換として考えないとおかしいのでは？ と思いました。
ですので、まずは離散フーリエ変換に直して、どうなるかを見ることにしました。

2. 実践編 1次元離散フーリエ変換

離散フーリエ変換については、主に以下の記事を参考にさせていただきました。

まずは数式を確認していきます。

\hat{g}(f) = \sum_{t=0}^{N-1}g(t)e^{-\frac{2{\pi}ift}{N}}

といっても、積分が総和になっただけで大きくは変わりません。
これならネットワークも大体同じに組めばよさそうです。
ただしeの指数に、時間≒標本点番号を正規化するようなNが追加されている点には注意が必要そうです。

ではネットワークを作っていきます。

まず、入力波形をforEachLoopで作るようにしました。この方が離散的な値であることが分かりやすかったからです。また、各種パラメータをコントロールするためのNullSOPを追加しました。

入力波形の周波数は、
離散フーリエ変換(DFT)の仕組みを完全に理解する によると

また，周波数は 1 Hz, 2 Hz の様に，離散的であり，1.5 Hz のような中途半端な周波数は使いません．これは単位時間内で周期的である必要があるので，単位時間内にピッタリ収まる波でないといけないからです．

とのことなので、
１Hz（1標本区間で波が一往復）の整数倍になるようにします。

ループの中では離散でないフーリエ変換同様、まず巻き付け処理をAttributeWrangleSOPで行います。
変更点は２つあります。

1つ目は入力波形のポイントごとのx位置の代わりにポイント番号を時間tとして使っている点です。
これはx位置をそのまま使うと小数が計算に入ってしまい、離散でなくなってしまうであろう……と考えたためです（あまり自信はない）。
2つ目は数式にあったように、標本点番号の正規化を追加した点です。

重心を求める処理は、ExtractCentroidSOPを使っていたのを念のためAttributeWrangleSOPによる処理に変更しています。

また、スペクトルを正しい振幅にするための処理を追加しています。

これは、以下のページの記載を参考に追加したものです。
（ただし、逆変換にかけるときはバイパスする必要があるようです）

ループを抜けたら重心の変化をプロットします。
ここまでの内容で、スペクトルは正しく得られているようです！

離散でないフーリエ変換の時と違い振幅の大きさもぴったり一致しています。

次に逆変換です。

スペクトルを入力として、離散フーリエ変換と同じ手順をもう一度行います。
前述の通り、「calcCorrectAmplitudeOnGraph」は一旦バイパスしておきます。
すると、完全にもとの波形と一致しました！

ただし注意点があって、巻き付け周波数の最大値>入力波形の長さ/2 となるようなケースで、スペクトルの内容が鏡映したようになってしまうようなケースがありました。

これは エイリアシング と呼ばれる現象のようで、 ナイキスト周波数 と呼ばれる周波数を境に起きてしまうようです。
この現象に対する完璧な対策は調べても見つけられなかったので、巻き付け周波数の最大値がナイキスト周波数を超えないようパラメータを調整して済ませています。

3. 実践編 2次元離散フーリエ変換で簡単な画像処理

さて、元の波形の完全再現ができたので、いよいよ当初の目標である、2次元離散フーリエ変換を使った画像の復元へとステップアップしていきます。

このあたりはこちらのスライドを参考にしました。

まずは画像を読み込んでポイント情報にするところから始めます。
本題から外れるので詳しい解説はしませんが、以下のような流れで行います。
解像度は128x128相当としておきましょう。

使用する画像はいらすとやで見つけたフーディニ（人名）のイラストにしました。

そしてこの段階で、 画像を4枚に複製 する処理を挟みます。

この理由については後程説明します。

次にスペクトルを求めるまでの処理を組んでいきます。

二次元離散フーリエ変換においては、forEachLoopを二重にしたうえで、巻き付け計算部分を画像のように修正する必要があります。

これまで同様、二次元離散フーリエ変換の数式をコードに起こしているだけです。

X_{u,v} = \sum_{h=0}^{H-1}\sum_{w=0}^{W-1}x_{h,w}e^{-2{\pi}j(\frac{uh}{H}+\frac{vx}{W})}

同様に総和を計算し……

スペクトルのプロット用のアトリビュートを付与し……

ループの外で実際にスペクトルを構築します。
2次元のスペクトル（3次元のグラフ）なので、二重ループの2つの巻き付け周波数をx,zとし、重心の変化をyとしてプロットする形になります。

そして肝となるスペクトルから元画像の復元は、
以下のようになっています。


二重ループで、2次元cos波の縦横の周波数を少しずつ変えながら足し合わせていきます。
足し合わせの際の係数として、先ほど得られたスペクトルの対応する座標（ここではポイント番号で代替）の振幅を使います。
注意点として、縦横それぞれの周波数が0になっているときはなぜか出力される画素値が二倍になってしまっていたので、例外的に2で割る処理を入れています。

さてこれでどうなるか見てみると……

元画像がしっかり再現できています！！ヤッター！！

そしてループの初期値と最大値を変更できるようにしているので、結果的に高周波カット、低周波カットもできるようになっています！
高周波をカットすると、一般に言われるように画像がしっかりボケるのを確認できました。

最後に先ほど説明をスキップした、画像を4枚に複製してから二重ループに投入している理由ですが、これは例のナイキスト周波数が絡んでいる部分なのかな、と考えています。
最初は画像1枚分で処理を行っていたのですが、下の図のようにぼんやりした画像になってしまっていました。

何かの情報が欠落しているように思えたので、試しに画像を4枚分に増やしてみると、

どうやら元画像を上下左右反転して重ねたような見た目に…？
さらに4枚の画像をそれぞれ反転してあげると、

無事元画像が元の画素値で復元できました！
正確な理屈はわかりませんが、上記のような結果をもとにこのような処理を入れています。

最終的なネットワークの全体像はこちらです。

まとめ

今回の検証を通じて、フーリエ変換がどのように成り立っているのか深く知ることができ、また数学に対する苦手意識をいくらか払拭できたように思います。
そして汎用性はさておきますが、自分なりに実装する方法も見出すことができました。

本記事を読んでいただいた皆さんにも、何かしら得るところがありましたら幸いです。

参考