フーリエ変換@離散フーリエ変換

参考文献及びWEB

参考文献

(1)「ディジタル信号処理技術」玉井徳みち、長島厚、藤田泰弘、若井修造著日経BP社
(2)「ディジタル信号処理の基礎」三上直樹著 CQ出版
(3)「C言語によるディジタル信号処理入門」三上直樹著 CQ出版
(4)「アナログ&ディジタルフィルタ入門」小野浩司著日刊工業
(5)「フーリエの冒険」ヒッポファミリークラブ
(6)「C言語で始める医用情報処理」小高知宏著 Ohmsha
(7)「量子力学の冒険」ヒッポファミリークラブ
(8)「最新ウェーブレット実践講座」戸田浩、章忠、川畑洋昭著 SoftBankCreative社
(9)「ウェーブレットによる信号処理と画像処理」中野宏毅、山本鎭男、吉田靖夫著共立出版

あわせて読む

フーリエ級数展開
フーリエ級数展開@射影と直交
フーリエ級数展開@複素数表示
フーリエ変換@基本
フーリエ変換@離散フーリエ変換
フーリエ変換@高速フーリエ変換
サンプリング周波数と基本周波数
畳込みについて
時間領域と周波数領域

フーリエ変換

$G(f) = \int_{-\infty}^{\infty} f(t)e^{-i 2\pi f t}dt$ (1-1)
$f(t) = \int_{-\infty}^{\infty} G(f) e^{i 2\pi f t}df$ (1-2)

意味:どんな複雑な波でも(たとえ周期がなくとも)単純な波一つ一つに分解できる

離散フーリエ変換

$G(\frac{1}{\tau N}\cdot n) = \sum_{k=0}^{N} f(k\cdot\tau)e^{-i 2\pi (\frac{n}{N}) \cdot k}$ (2-2)
$f(k\cdot\tau) = \sum_{n=0}^{\infty} G(\frac{1}{\tau N}\cdot n) e^{i 2\pi (\frac{n}{N})\cdot k}$ (2-3)

N:値の数(何点取るか？)
$\tau$ : 値を読んだ間隔(サンプリング周期)
k:何番目の値か(時間波形の)？
n:基本周波数の何倍か？

標本化(sampling)

しかしながら、実際のデータは連続でなくデータとデータの間に間があります。
要するに離散的でしか表すことができません。

その場合ある一定の時間間隔(サンプリング周期)で標本値(サンプル値)を取らなければなりません。
これを標本化(sampling)と言います。

例えば、
いま250Hzのサンプリング周波数でデータをサンプリングしようとすると

1/250=0.004[s]
     =4 x 10-3 sec
     =4[ms]

の間隔でデータを取得していくということです。

図で表すと下図のようになります


                         *
     |         *         *   *
     |    *    *    *    *   *
     |____|____|____|____|___|
    0     4    8    12   16  20 [msec]

時間領域上でサンプリング周期を短くしていけばいくほど、
周波数領域上で見れる周波数帯が広がるということに注目してください。

□ナイキストのサンプリング定理

時間データを周波数変換し、周波数領域にしたとき、
サンプリング周波数が250Hzならば、1/2の125Hzまでが観測できる最大の周波数となります。
これをナイキストのサンプリング定理と呼びます。

量子化(quantitize:クオンティタイズ)

また信号の振幅はアナログ的な連続量のままでは無限の分解能をもつはずです。
これを必要な分解能に応じた値に振幅を変換することを量子化(quantitize:クオンティタイズ)と言います。

例えば
0から5Vまでを10bit=2^10=1024で量子化すると、

5/1024=0.0049

となり、
値が1のときは0.0049[V]の電圧ということがわかり、
値が1024のときは、5[V]の電圧ということが分かるのです。


0.0147-|                       *
0.0098-|         *         *   *
0.0049-|    *    *    *    *   *
       |____|____|____|____|___|
      0     4    8    12   16  20 [msec]

離散フーリエ変換を導く

では離散化されたデータをフーリエ変換するにはどうすれば良いのでしょう。
実際のフーリエ変換は離散値を扱うため少々式を変更する必要がありそうです。

$G(f) = \int_{-\infty}^{\infty} f(t)e^{-i 2\pi f t}dt$ (1-1)

今、離散間隔(サンプリング周期と呼びます)を $\tau$ とします。

またとびとびの点で表される波を読み込んでフーリエ変換をかけようとしたとき
それが飛び飛びの値を取るなら飛び飛びの値の最大個数が分かるはずです。
つまり0～N個までの飛び飛びの値を読んでそれを一周期としてフーリエ変換を
かけようと考えることができるはずです。

そうすると波の1周期を $T=\tau*N$ と考えることができます。

また、
関数 f(t) は　 $f(k*\tau)$ (ただし $k=0,1,2,3,\cdots,N$ )
と考えることができます。

積分は離散間隔を限りなく0に近づけたときの面積です。
ところが今は離散間隔は $\tau$ という値を持っていますから
G(f) の右の式は

$= \sum_{k=0}^{N} f(k\cdot\tau)e^{-i 2\pi f \cdot (k\tau)}dt$ (2-1)

になることが分かります。
面積ですから幅 $\tau$ と高さ $f(k\tau)$ の掛け算をして0からN-1まで足し合わせることになるのです。

ところで全体の周期が $T=\tau*N$ ですから
基本周波数は $\frac{1}{(\tau*N)}$ になります。

ある基本周波数の整数倍の振幅を取り出せるのがフーリエ変換の定義でした。
それがまた離散的にしか取り出すことができなくなってしまった。
要するに今までは基本周波数が0に限りなく近かった。
今度は基本周波数が分かっていて、しかも離散的だから取り出すときは
$\frac{1}{(\tau*N)} \times n$
$(n=0,1,2,3,\cdots,\infty)$
の周波数成分が取り出せるということになる。
つまり

$G(\frac{1}{\tau N}\cdot n) = \sum_{k=0}^{N} f(k\cdot\tau)e^{-i 2\pi (\frac{1}{\tau N}\cdot n) \cdot (k\tau)}dt$
$G(\frac{1}{\tau N}\cdot n) = \sum_{k=0}^{N} f(k\cdot\tau)e^{-i 2\pi (\frac{n}{N}) \cdot (k)}dt$ (2-2)

として各周波数成分が取り出せるはずです。

ここで

N:値の数(何点取るか？)
$\tau$ : 値を読んだ間隔(サンプリング周期)
k:何番目の値か(時間波形の)？
n:基本周波数の何倍か？

となります。
これを離散フーリエ変換と言います。

離散フーリエ変換をプログラミングする

DFTとFFTのテスト
 DFTとFFTのソースコード

$G(\frac{1}{\tau N}\cdot n) = \sum_{k=0}^{N} f(k\cdot\tau)e^{-i 2\pi (\frac{1}{\tau N}\cdot n) \cdot (k\tau)}dt$
を実装すればよい。

/**
  * 離散周波数変換
  */
public void dfft_time(double dt,double df){
  double yrr=0,yii=0;
  //窓関数の適用
  window(yr);
   
  //周波数を整数倍を表しているG(n) (sampling/N)*n=df*nに相当するものは表示時にする
  for(int i=0; i<xr.length; i++){
     yrr=0;
     yii=0;
     //k=0..Nに相当する部分k\tauに相当する部分はj*dt 1/(\tau N) x nに相当する部分はi*dfとなる       
     for(int j=0; j < xr.length ; j++){
        yrr = yrr+xr[j]*Math.cos(2*Math.PI*(i*df)*(j*dt));
        yii = yii-xr[j]*Math.sin(2*Math.PI*(i*df)*(j*dt));
     }
     yr[i]=yrr;
     yi[i]=yii;
   }
 }
/**
 表示時
*/ 
for(int i = 0; i < drawBuffer.length - 1 ; i++){
 //表示時に基本周数数dfを掛ければ良い  
 mc.drawLine(i*df, drawBuffer[i],(i + 1)*df, drawBuffer[i + 1],2);
}