AlxisDataBlog-プログラミング


MathJax = {
  chtml: {
    matchFontHeight: false
  },
  tex: {
    inlineMath: [['$', '$']]
  }
};






以前フーリエ変換をやりましたね。今回は説明し損ねていたフーリエ変換の高速化の話です。離散フーリエ変換(DFT)というものはご存知でしょうか。以下のような配列 $f(j)$ を配列 $F(k)$ へ変換するアルゴリズムのことです。

\begin{align}
F(k)=\sum_{j=0}^{N-1}f(j)e^{-2\pi ijk/N},\ k=0,1,\dots,N-1
\end{align}

フーリエ級数の係数を取り出すための積分演算が和で置き換わっていますね。DFT はデジタル信号処理にはなくてはならない非常に重要な技術です。しかし計算量に着目すると、1 つの $k$ に対して $F(k)$ を計算するためにはループ文を走らせなければなりません。したがってすべての $k$ に対して $F(k)$ を計算して変換を完了させるためには、2 重のループ文を実行し、ループ文の中身を $N^2$ 回計算する必要があります。


このループ回数を劇的に減らしたアルゴリズムが、これから説明する高速フーリエ変換(FFT)です。FFT の登場により信号処理分野は大きく発展しました。


FFT は $N$ の素因数分解を使います。説明を簡単にするために 2 つの素数 $p,q$ により $N=pq$ となっている前提で話を進めましょう。


配列 $f,F$ のレコード数を、$p\times q$ の 2 次元配列と考えてみましょう。$p$ と $q$ のどちらを行と考えるかにより 2 通りの表現方法があります。

\begin{align}
j=qx_1+x_2,\ k=y_1+py_2
\end{align}

このような 2 次元配列として考えると $jk/N$ は以下のようになります。

\begin{align}
\dfrac{jk}{N}&=\dfrac{pqx_1y_2+qx_1y_1+px_2y_2+x_2y_1}{pq}\\
&=x_1y_2+\dfrac{x_1y_1}{p}+\dfrac{x_2y_2}{q}+\dfrac{x_2y_1}{pq}
\end{align}

これを指数関数の引数に代入すると

\begin{align}
e^{-2\pi ijk/N}&=e^{-2\pi ix_1y_2}e^{-2\pi i x_1y_1/p}e^{-2\pi ix_2(y_2/q+y_1/pq)}\\
&=e^{-2\pi ix_1y_1/p}e^{-2\pi ix_2(y_2/q+y_1/pq)}
\end{align}

$x_1y_2$ は整数であることと変数 $x_1,x_2$ の分離という点に注意して式を追ってみてください。変数分離によって DFT の総和計算も分離されます。

\begin{align}
F(y_1+py_2)&=\sum_{x_1=0}^{p-1}\sum_{x_2=0}^{q-1}f(qx_1+x_2)e^{-2\pi ix_1y_1/p}e^{-2\pi ix_2(y_2/q+y_1/pq)}\\
&=\sum_{x_2=0}^{q-1}\Big(\sum_{x_1=0}^{p-1}f(qx_1+x_2)e^{-2\pi ix_1y_1/p}\Big)e^{-2\pi ix_2(y_2/q+y_1/pq)}
\end{align}

変数分離によって計算量が減っていることに気が付くでしょうか。実際に数えてみましょう。内側のループ文はすべての $x_2,y_1$ に対して実行しなければならないため、計算回数は $pN$ です。外側のループ文に関しても同様にすべての $y_1,y_2$ に対して実行しなけれなならず、計算回数は $qN$ です。したがって合計で $(p+q)N$ 回の計算が必要になります。


DFT と FFT を比較すると $N$ の部分が $p+q$ になって減少していることが分ります。このことから更に $N$ の素因数分解が細かいほど高速化の効果は高い、すなわち $N$ が 2 の冪乗となる場合に最も高い高速化効果が得られることがわかります。


ご興味あれば $N$ が 2 の冪乗となる場合の FFT を勉強してみると良いでしょう。次回は FFT を用いた畳込み演算の高速化の話をします。(written by J.Nishigata)


Copy Right 2020-現在 ALXIS DATA, Inc. All Right Reserved.

高速フーリエ変換、高速化の秘密を知ろう

AlxisDataBlog-画像処理


MathJax = {
  chtml: {
    matchFontHeight: false
  },
  tex: {
    inlineMath: [['$', '$']]
  }
};






前回の宣言通りニッチな話をします。広い文脈で捉えるならばデコンボリューションの応用事例として考えることもできますが、あえてデコンボリューションを使わないアイデアを紹介します。


最小エレメント幅以上の径でぼやけたバーコードは、2 値化による読み取り困難です。バーコードの利用現場では、カメラの性能やオペレータの操作の不安定性から、しばしば最小エレメント幅以上のぼやけが発生します。このような状況でも安定してバーコードを読み取るためには、ぼやけに強い符号化が必要になります。


バーコードのラインセンシング画像を関数と見立てたとき、関数の表現を変えると聞いてピンとくるものは何でしょうか。おそらくフーリエ変換が思い浮かぶでしょう。


しかし、フーリエ変換は Blur に対して劣化が著しく、不安定です。Blur といえば、単純に積分してしまえは情報量は少ないですが安定しています。

\begin{align}
\lambda_1=\dfrac{1}{L}\int_0^Lf_0(x)\,dx
\end{align}

同様にして

\begin{align}
&f_n(x)=\dfrac{1}{L}\int_0^xf_{n-1}(y)\,dy-\dfrac{\lambda_nx}{L},\\
&\lambda_{n+1}=\dfrac{1}{L}\int_0^Lf_n(x)\,dx
\end{align}

で計算される数列 $\lambda$ も安定です。したがって数列 $\lambda$ が、バーコードの示す各キャラクタに対して良く分離しているのであれば、安定した符号化として使えます。実際に使う場合は以下のステップを実行します。


1. 理論的な各キャラクタから数列 $\lambda$ を算出(これを基準列と呼ぶ)
2. 実測されたキャラクタから数列 $\lambda$ を計算して各基準列と比較
3. 最も「近い」基準列のキャラクタを返す


最後の「近い」ことの定義には若干の癖があります。基準列との差分の $L_2$ をとったとしても、なかなか上手くいきません。その理由は、比較する基準列に依存して距離の縮尺が変わるからです。この縮尺も含めて最初のステップで計算しておけば良いです。


なお、実際に実装する場合には数列 $\lambda$ の計算を

\begin{align}
(-1)^{n-1}\lambda_n=\dfrac{1}{(n-1)!L^n}\int_0^Lx^{n-1}f_0(x)\,dx-\dfrac{\lambda_1}{n!}+\dfrac{\lambda_2}{(n-1)!}+\dots+(-1)^{n-1}\dfrac{\lambda_{n-1}}{2!}
\end{align}

のように実行すると計算量が少なくなります。PassMark 700 前後の環境で単位キャラクタ当たり 11μs 強で実行できます。この符号化を使うと、なんと最小エレメント幅の 2 倍の Blur に対しても安定してバーコードを読み取ることができます。次回もニッチな話です。(written by J.Nishigata)


Copy Right 2020-現在 ALXIS DATA, Inc. All Right Reserved.

FIR フィルタリングを FFT で高速化してみよう