サルでも分かるFFT（1）

「サルでも分かる」と書きましたが、始めに断っておきます。サルは私です。

大学時代、信号処理の授業は「優」で通っていたはずなのに自分の無恥さを知りました。

今日は、フーリエ変換の概念とDFTのプログラムサンプルを書きます。（間違いがありましたら、ご指摘ください。）

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

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

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

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

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

Xkが入力データ（ある時刻における振幅）、fjがDFTの結果で周波数に直したときの特性を示します。

具体的には、

ABS(fj) がその周波数のスペクトル、ArcTan (Img(fi)/Re(fj))がその周波数成分の位相を表します。（fjは複素数であることに注意）

と、ここまでは知ってました。今日改めて気付いたのは、サンプリング点がN個の場合は、N個の周波数成分に分解できるということです。連続信号を離散で取っているので、周波数成分への分割の方法は複数パターン可能（特にNが小さければ無数にある）だと思いますが、DFTの結果って、どういう分割になっているのだろう。。。

あとfjのjって何？？っていう疑問がわきました。これは、サンプリング間隔周波数のj倍の周波数成分と予想してるのですが、どうでしょう・・・。（うわ、全然分かってないな。。あとで調べよう。 OR　誰か教えてください笑）

以下、今日書いてみたDFTのコード。

複素数の扱いが簡単なPythonで書きました。（例題サンプルをいくつか見つけて検算したので合っているはず。）

view plainprint?

1. import sys  
2. import math  
3.   
4. # input data  
5. x = []  
6. for line in sys.stdin:  
7.   x.append(float(line))  
8.   
9. # construct matrix  
10. n = len(x)  
11. W = [[0 for j in range(n)] for i in range(n)]  
12. for i in range(n):  
13.   for j in range(n):  
14.       W[i][j] = math.cos(2*i*j*math.pi/n) - math.sin(2*i*j*math.pi/n)*1J  
15.   
16. # calc  
17. f = [0] * n  
18. for i in range(n):  
19.   for j in range(n):  
20.       f[i] += W[i][j] * x[j]  
21.   f[i] /= n  
22.   
23. # display result  
24. for i in range(n):  
25.   print f[i]  

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]