Search on the blog

2011年5月16日月曜日

サルでも分かるFFT(1)

「サルでも分かる」と書きましたが、始めに断っておきます。サルは私です。
大学時代、信号処理の授業は「優」で通っていたはずなのに自分の無恥さを知りました。
今日は、フーリエ変換の概念とDFTのプログラムサンプルを書きます。(間違いがありましたら、ご指摘ください。)

 まず、フーリエ級数展開について。一言でいうと、「すべての周期関数は、複数のサイン波の重ね合わせで表現できる。」です。ここでポイントなのは、対象が”周期関数”であるということです。

 そして、フーリエ級数展開を非周期関数に適用可能にしたものがフーリエ変換です。実際の信号は周期性がないものがほとんどですので、フーリエ変換が必要になるわけです。

 次に、離散フーリエ変換(DFT)。実際の信号(音や光など)は連続関数です。しかし、コンピュータ上で連続な信号を扱うことは難しいです。なので、サンプリングという操作を行います。実際には連続な信号を0.1msに1個とかいう間隔で信号の値を取得して離散信号として扱います。この離散化された信号(デジタル信号)に対するフーリエ変換がDFTです。

 そして、いよいよFFT。FFTはDFTを高速に計算するためのアルゴリズムです。別にFFTなくてもデジタル信号のフーリエ変換はできます。どれくらい速くなるかというと、O(N^2)がO(Nlog N)まで減ります。基本ソートとクイックソートくらい違います。

 で、DFTを実装してみます。Wikipediaにある下の式を使います。

 f_j = \sum_{k=0}^{N-1} x_k \exp\left( - \frac{2 \pi i j k}{N} \right) \qquad j = 0,\dots,N-1

Xkが入力データ(ある時刻における振幅)、fjがDFTの結果で周波数に直したときの特性を示します。
具体的には、
ABS(fj) がその周波数のスペクトル、ArcTan (Img(fi)/Re(fj))がその周波数成分の位相を表します。(fjは複素数であることに注意)

と、ここまでは知ってました。今日改めて気付いたのは、サンプリング点がN個の場合は、N個の周波数成分に分解できるということです。連続信号を離散で取っているので、周波数成分への分割の方法は複数パターン可能(特にNが小さければ無数にある)だと思いますが、DFTの結果って、どういう分割になっているのだろう。。。
あとfjのjって何??っていう疑問がわきました。これは、サンプリング間隔周波数のj倍の周波数成分と予想してるのですが、どうでしょう・・・。(うわ、全然分かってないな。。あとで調べよう。 OR 誰か教えてください笑)

以下、今日書いてみたDFTのコード。
複素数の扱いが簡単なPythonで書きました。(例題サンプルをいくつか見つけて検算したので合っているはず。)
import sys
import math

# input data
x = []
for line in sys.stdin:
x.append(float(line))

# construct matrix
n = len(x)
W = [[0 for j in range(n)] for i in range(n)]
for i in range(n):
for j in range(n):
W[i][j] = math.cos(2*i*j*math.pi/n) - math.sin(2*i*j*math.pi/n)*1J

# calc
f = [0] * n
for i in range(n):
for j in range(n):
f[i] += W[i][j] * x[j]
f[i] /= n

# display result
for i in range(n):
print f[i]

0 件のコメント:

コメントを投稿