lynxeyedの電音鍵盤

MBDとFPGAと車載で使うデバイスの備忘録

画像認識におけるフーリエ変換とフーリエ記述子(その1)

はじめに

一連のFPGA関連ネタとして、画像認識を扱おうと思います。
ここでは周波数スペクトルを推定する離散フーリエ変換からはじまり、フーリエ記述子、ウェーブレット変換、ニューラルネットワークに至るまで取り扱うつもりです。

画像の特徴を抽出する

 入力画像から特定の形状を抽出して、基準画像と比較することは画像認証ではよく行う方法です。そのためには、入力画像データから図形を抽出し、数式化する必要があります。
 例えば一辺の長さが2の正方形を画像の中から抽出し数式化しなければならないとします。もし図1の様に、画像の縦横に対して図形が平行なら、数式化は難しくありません。
f:id:Lynx-EyED:20130612191453p:plain
図1:一辺が2の正方形

この正方形は以下の式で表せます。
f:id:Lynx-EyED:20130612235544p:plain

しかし、図2の様に正方形に傾きがある場合、数式化は難しいかもしれません。不可能ではないですが、その数式から正方形であるかどうか判別するのに時間がかかるはずです。

f:id:Lynx-EyED:20130612191505p:plain
図2:傾きのある正方形

このように傾きがある場合でも、画像中に正方形が存在する事を数式として表現できるでしょうか。
一つの答えは、図3のように正方形上を移動する点Pを用意する事です。
f:id:Lynx-EyED:20130612191514p:plain
図3:正方形上を移動する点P

言葉で表現すると以下の様になります。

  1. 距離2だけ進み、開始時のベクトル方向より正の方向(反時計回り)に90度回転
  2. 距離2だけ進み、開始時のベクトル方向より正の方向(反時計回り)に180度回転
  3. 距離2だけ進み、開始時のベクトル方向より正の方向(反時計回り)に270度回転
  4. 距離2だけ進み、開始時のベクトル方向より正の方向(反時計回り)に360度回転


このように表現すると正方形の傾きは関係なくなります。
グラフ化してみましょう。点Pが移動する距離をs、開始時のベクトルと現時点でのベクトルのなす角度をθとします。先ほどdegreeで角度を表現しましたが、便宜上radian表示に直しました。
f:id:Lynx-EyED:20130612191530p:plain
図4:点Pの軌跡をグラフ化

少し違和感のあるグラフになりました。点Pはそもそも正方形上を移動しているので2x4=8だけ進むと周期的に元の場所に戻って来ているはずです。

そこで、点Pの軌跡を「正規化」します。*1
先ほどの様に言葉で表現するとこうなります。

  1. 距離2だけ進み、正の方向(反時計回り)に90度回転
  2. (1)を4回繰り返す


グラフで表現しましょう。「正規化」しているので少し複雑になります。

θの正規化偏角関数をθN(s)、
正方形の全長をL、
開始時のベクトル方向をθ(0)、
現在の角度とθ(0)とのなす角をθ(s)

とすると、
f:id:Lynx-EyED:20130612192029p:plain
となり、グラフ化すると図5の様になります。
f:id:Lynx-EyED:20130612191540p:plain
図5:正方形に対する正規化偏角関数

図5から、正方形を傾きに依存せず周期的な関数で表せる事がわかります。*2

このように、図形の持つ角度(または周波数成分)から、どのような形状をしているかを知る(推定する)ことが可能であり、情報圧縮や特徴抽出に利用されます。

 取得したデータにどの周波数成分がどの程度含まれているかを知る上でフーリエ変換は有効と言えます。

フーリエ変換

周期T[s]の周期関数の基本角周波数ω[rad/s]は
f:id:Lynx-EyED:20130612234001p:plain…(1)

となります。
フーリエ級数における基本的な概念は、
「全ての周期信号はバイアス成分とω[rad/s]およびその整数倍の角周波数をもつ三角関数の和として表現できる」
というものです。
ある時間変化する周期信号 s(t)に対して、
f:id:Lynx-EyED:20130612234014p:plain…(2)

が成り立ちます(実フーリエ級数展開)。
ただし三角関数の直交性から
f:id:Lynx-EyED:20130612234023p:plain…(3)

となります。
式(2)を複素指数関数で表現すると、以下の様になります。
f:id:Lynx-EyED:20130612234031p:plain…(4)


式(4)は、分解して考えると
f:id:Lynx-EyED:20130612234040p:plain…(5)


式(5)と式(2)の係数をオイラーの公式を考慮しつつ対応させて考えると、以下の様になります。
f:id:Lynx-EyED:20130612234048p:plain…(6)


式(3)でも示しましたが、フーリエ級数展開は周期関数に対して有効な手がかりです。しかし非周期信号や孤立波を対象にする場合はこの方法は有効ではありません。どうすれば良いでしょうか?

非周期信号を周期無限大の周期関数とみなせばよいのです。*3

具体的には積分を±無限大の範囲で行う事になります。
しかし数学的にはその手法をとる事は可能ですが、コンピュータでは離散化された信号が対象であり、有限回で計算を終える必要があります。そのために

時間Tのあいだにサンプリングした全データNが周期的に繰り返されるものとみなす

という工夫をして、この問題を解決*しようと*しています。(なぜ「解決しています」では無いのかは後述) これが離散フーリエ変換の考え方になります。

離散フーリエ変換(DFT)

前述のフーリエ級数展開では、信号s(t)は連続時間tの間に変化する信号であったので時間軸に対して積分を行っていました。ここではN[個]のサンプリングデータに対して演算を行います。従って複素フーリエ係数を求める際には時間tではなく、i=0番目から(N-1)番目のサンプリングデータに対して積分(正確に言うと総和)を求めます。

アナログ入力信号: s(t)
サンプリング周期: τ[s/回]
総サンプリングデータ: N[回]またはN[個]
サンプリング時間: T[s]

とする。s(t)を標本化したデータのうちi番目のものs[i]は、

f:id:Lynx-EyED:20130614223033p:plain …(7)

時間t[s]は標本化されたiτ[s]に置き換えられます。
また、
f:id:Lynx-EyED:20130614223108p:plain …(8)

f:id:Lynx-EyED:20130614223119p:plain …(9)


前述した通りN[個]のサンプリングデータに対して演算を行うため、式(5)で求めた複素フーリエ係数は離散化したデータにおいて、
f:id:Lynx-EyED:20130614223132p:plain …(10)

式(10)が離散フーリエ変換と呼ばれるものになります。また、式(4)を直接離散化すると、 逆離散フーリエ変換(IDFT)が得られ、

f:id:Lynx-EyED:20130614223144p:plain …(11)


となります。しかしT=Nτはサンプリング時間であり、信号の周期と必ずしも一致しません。例えば図6のようにサンプリング周期Tと実際の波形の周期T’が異なっている場合を考えましょう。
f:id:Lynx-EyED:20130614223205p:plain
図6:サンプリング周期と実際の周期が異なる場合

この波形に対してDFTを行い、その後IDFTを行うと、図7の様に不連続な波形となってしまい、スペクトルには多くの高次成分が現れてしまいます。
f:id:Lynx-EyED:20130614223231p:plain
図7:不連続な波形として再現されてしまい、多くの高次成分が現れる

この様な現象を抑えるためハミングや、ハニング関数が用いられます。

ちょっと長いので位相回転因子やFFTの話題、ハニング窓関数のお話はまた次回に。


◆参考文献
酒井幸市 (2003);ディジタル画像処理入門 OpenDesignBooks, pp.151-156. CQ出版社
酒井幸市 (2007);改訂版 ディジタル画像処理の基礎と応用 ディジタル信号処理シリーズ, pp.119-125. CQ出版社
後藤尚久 (2001);なっとくする電気数学, pp.99-119. 講談社
佐野理 (1992);キーポイント 微分方程式, pp.108-138. 岩波書店
昌達慶仁 (2010);圧縮処理プログラミング, pp.318-319. ソフトバンククリエイティブ

*1:詳細は後述しますが、不連続点が多いと収束が遅くなるのでなるべく単純な周期関数になるようにします

*2:詳細は後述する、有限個のフーリエ係数で図形の特徴を記述するフーリエ記述子の項で扱います

*3:厳密に言うと周期無限大=非周期関数