フーリエ変換@高速フーリエ変換

参考文献及びWEB

参考文献

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

あわせて読む

ソースコード

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

離散フーリエ変換

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

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

高速フーリエ変換

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

いま、簡単化のために $\Delta{t}=1$ とすると
$G(\frac{n}{N}) = \sum_{k=0}^{N-1} f(k)e^{-i 2\pi (\frac{n}{N}) \cdot (k)}$

Nは定数となりますので、

$W = e^{-i2\pi/N}$

と表します。を回転子と言います。

すると式は

$G(\frac{n}{N}) = \sum_{k=0}^{N-1} f(k)W^{nk}$

となります。ではN=8のときを計算してみましょう。

$G(\frac{n}{8}) = \sum_{k=0}^{7} f(k)W^{nk}$

これを n=0～7まで計算しましょう。

$G(0/8) = f(0)W^{0} + f(1)W^{0} + f(2)W^{0} + f(3)W^{0} + f(4)W^{0} + f(5)W^{0} + f(6)W^{0} + f(7)W^{0}$
$G(1/8) = f(0)W^{1*0} + f(1)W^{1*1} + f(2)W^{1*2} + f(3)W^{1*3} + f(4)W^{1*4} + f(5)W^{1*5} + f(6)W^{1*6} + f(7)W^{1*7}$
$G(2/8) = f(0)W^{2*0} + f(1)W^{2*1} + f(2)W^{2*2} + f(3)W^{2*3} + f(4)W^{2*4} + f(5)W^{2*5} + f(6)W^{2*6} + f(7)W^{2*7}$
$G(3/8) = f(0)W^{3*0} + f(1)W^{3*1} + f(2)W^{3*2} + f(3)W^{3*3} + f(4)W^{3*4} + f(5)W^{3*5} + f(6)W^{3*6} + f(7)W^{3*7}$
$G(4/8) = f(0)W^{4*0} + f(1)W^{4*1} + f(2)W^{4*2} + f(3)W^{4*3} + f(4)W^{4*4} + f(5)W^{4*5} + f(6)W^{4*6} + f(7)W^{4*7}$
$G(5/8) = f(0)W^{5*0} + f(1)W^{5*1} + f(2)W^{5*2} + f(3)W^{5*3} + f(4)W^{5*4} + f(5)W^{5*5} + f(6)W^{5*6} + f(7)W^{5*7}$
$G(6/8) = f(0)W^{6*0} + f(1)W^{6*1} + f(2)W^{6*2} + f(3)W^{6*3} + f(4)W^{6*4} + f(5)W^{6*5} + f(6)W^{6*6} + f(7)W^{6*7}$
$G(7/8) = f(0)W^{7*0} + f(1)W^{7*1} + f(2)W^{7*2} + f(3)W^{7*3} + f(4)W^{7*4} + f(5)W^{7*5} + f(6)W^{7*6} + f(7)W^{7*7}$

これは、

$G(0/8) = f(0)W^{0} + f(1)W^{0} + f(2)W^{0} + f(3)W^{0} + f(4)W^{0} + f(5)W^{0} + f(6)W^{0} + f(7)W^{0}$
$G(1/8) = f(0)W^{0} + f(1)W^{1} + f(2)W^{2} + f(3)W^{3} + f(4)W^{4} + f(5)W^{5} + f(6)W^{6} + f(7)W^{7}$
$G(2/8) = f(0)W^{0} + f(1)W^{2} + f(2)W^{4} + f(3)W^{6} + f(4)W^{8} + f(5)W^{10} + f(6)W^{12} + f(7)W^{14}$
$G(3/8) = f(0)W^{0} + f(1)W^{3} + f(2)W^{6} + f(3)W^{9} + f(4)W^{12} + f(5)W^{15} + f(6)W^{18} + f(7)W^{21}$
$G(4/8) = f(0)W^{0} + f(1)W^{4} + f(2)W^{8} + f(3)W^{12} + f(4)W^{16} + f(5)W^{20} + f(6)W^{24} + f(7)W^{28}$
$G(5/8) = f(0)W^{0} + f(1)W^{5} + f(2)W^{10} + f(3)W^{15} + f(4)W^{20} + f(5)W^{25} + f(6)W^{30} + f(7)W^{35}$
$G(6/8) = f(0)W^{0} + f(1)W^{6} + f(2)W^{12} + f(3)W^{18} + f(4)W^{24} + f(5)W^{30} + f(6)W^{36} + f(7)W^{42}$
$G(7/8) = f(0)W^{0} + f(1)W^{7} + f(2)W^{14} + f(3)W^{21} + f(4)W^{28} + f(5)W^{35} + f(6)W^{42} + f(7)W^{49}$

こうなります。
64回の積和演算をしなくはいけません。これを減らすことができるのでしょうか？
表にしてみましょう。

図1-1

先ほどWは回転子と言うと言いました。
ではなぜ回転子と言うのでしょうか？
ここでWは何だったのかを思い出してみましょう。

$W = e^{-i2\pi/N} = \cos(2\pi/N) - i\sin(2\pi/N)$

は円周状の点yを正弦波,xを余弦波で表されるものだから結局周期 $2\pi$
でぐるぐる回っている現象を表している。

円と正弦波のアプリも参考

そして今は、 N=8 なのでこの回転は図のように8分割されて回転していることになるのです。
ということは

$W^0 = W^8, W^1 = W^9, \cdots$ と書きかえられるので、

図1-2

となります。

偶数項と奇数項に分ける

図1-3

すると偶数組みの方は上半分と下半分が同じです。
これを使えば計算量を半分にできそうです。
でも奇数組がいまいちです。
ひとまず、偶数組と奇数組に分けてみましょう。

$G(\frac{n}{N}) = \sum_{k=0}^{N-1} f(k)W^{nk}$
$= \sum_{k=0}^{N/2-1} f(2k)W^{n\cdot2k} + \sum_{k=0}^{N/2-1} f(2k+1)W^{n\cdot(2k+1)}$

ここで、
f(2k) = p(k)
f(2k+1) = q(k)

として書き直します。

$G(\frac{n}{N}) = \sum_{k=0}^{N/2-1} p(k) W^{n\cdot2k} + \sum_{k=0}^{N/2-1} q(k) W^{n\cdot(2k+1)}$
$= \sum_{k=0}^{N/2-1} p(k) W^{n\cdot2k} + \sum_{k=0}^{N/2-1} q(k) W^{2nk}W^{n}$
$= \sum_{k=0}^{N/2-1} p(k) W^{n\cdot2k} + W^{n}\sum_{k=0}^{N/2-1} q(k) W^{2nk}$

ここで
p(k),q(k) は、 $0,\cdots,N/2-1$ までで
f(2k) = p(k), f(2k+1)=q(k) ですから
$p(0) = f(0),p(1) = f(2),\cdots$ という関係を忘れないようにすること。

ところで奇数項は W^n が飛び出たので W^n を引いてもう一度図を書き直しましょう。

図1-4

これで
奇数　偶数とも上半分だけで計算することができるようになりました。
計算回数は

偶数16回+奇数16回 + W^nをかける回数4回=36回

で計算することができるようになりました。

$W^{2nk} = W^{'nk}$ とおく

$W = e^{-i2\pi/N}$
$W^{nk} = e^{(-i2\pi/N)nk}$

$W^{2nk} = e^{(-i2\pi/N)2nk} = e^{(-i2\pi/(N/2))nk}$

となって $W^{2nk}$ になると円は4等分になります。そこで

$W^{2nk} = W^{'nk}$

とすることで $p=0,\cdots,3$ までの式で表すことができます。つまり

$G(\frac{n}{N}) = \sum_{k=0}^{N/2-1} p(k) W^{n\cdot2k} + W^{n}\sum_{k=0}^{N/2-1} q(k) W^{2nk}$
$G(\frac{n}{N}) = \sum_{k=0}^{N/2-1} p(k) W^{'nk} + W^{n}\sum_{k=0}^{N/2-1} q(k) W^{'nk}$

と表すことができます。
ではそれに習いをにしてしまいましょう。
の肩を2で割れば良いのです。

図1-5

今度はそう上半分と下半分が同じなので
上半分だけ考えてこれをまた2等分できるかどうかやってみる。

偶数と奇数に分ける

図1-6

式は
$G(\frac{n}{N}) = \sum_{k=0}^{N/2-1} p(k) W^{'nk} + W^{n}\sum_{k=0}^{N/2-1} q(k) W^{'2nk}$

$G(\frac{n}{N}) =$ $\sum_{k=0}^{N/4-1} p(2k) W^{'2nk} + \sum_{k=0}^{N/4-1} p(2k+1) W^{'n(2k+1)}$
$+ W^{n} \Big{ \sum_{k=0}^{N/4-1} q(2k) W^{'2nk} + \sum_{k=0}^{N/4-1} q(2k+1) W^{'n(2k+1)} \Big}$

ここで

$W^{'2nk+n} = W^{'2nk} \cdot W^{'n}$

なので先ほどと同じように $W^{'n}$ を前に出す。

$G(\frac{n}{N}) =$
$\sum_{k=0}^{N/4-1} p(2k) W^{'2nk} + W^{'n} \sum_{k=0}^{N/4-1} p(2k+1) W^{'2nk}$
$+ W^{n} \Big{ \sum_{k=0}^{N/4-1} q(2k) W^{'2nk} + W^{'n} \sum_{k=0}^{N/4-1} q(2k+1) W^{'2nk} \Big}$

となる。

とおく

$G(\frac{n}{N}) =$ $\sum_{k=0}^{N/4-1} a(k) W^{'2nk} + W^{'n} \sum_{k=0}^{N/4-1} b(k) W^{'2nk}$
$+ W^{n} \Big{ \sum_{k=0}^{N/4-1} c(k) W^{'2nk} + W^{'n} \sum_{k=0}^{N/4-1} d(k) W^{'2nk} \Big}$

また同じように $W^{'n}$ でくくって引くと

図1-7

図のようになります。
ここで a(0) とか d(1) はなんだったのでしょう。
思い出してみます。

図1-8

ここで
W^n は8等分のWだから45度位相がずれる波となる。

$W^{'n}$ は4等分のWだから90度位相がずれる波となります。

だから d(1) なんかは、まず45度ずれてさらに90度ずれる波を求めていることになります。

結局

: 4回
$W^{'n}(b(0),b(1))$ 4回＋2回
$W^{n}(c(0),c(1))$ 4回＋2回
$W^{n}(W^{'n}(d(0),d(1)))$ 4回＋2回＋2回

計24回の計算時間の短縮になるというわけです。

実際にプログラミングする

シグナルフローグラフを作る

演算の流れをもう一度よく見てシグナルフローグラフを作ってみましょう。

をで表す

図1-5

図は,2個の4点DFTになっていました。

偶数の4点DFTを P(n) ,奇数の4点DFTを Q(n) として表しましょう。
ためしに最初の G(n/8) の偶数部分を P(n/4) で表すと
n=0 の場合は、

$P(0/4) = p(0)\cdot W^{'0} + p(1)\cdot W^{'0} + p(2)\cdot W^{'0} + p(3)\cdot W^{'0}$

となります。

$W^{'}$ は4等分された回転子ですのでとしましょう。

はとしましょう。

書き直すと

$P(0/4) = p(0)\cdot W_4^{0} + p(1)\cdot W_4^{0} + p(2)\cdot W_4^{0} + p(3)\cdot W_4^{0}$

です。全部表して見ましょう。

〇偶数

$P(0/4) = p(0)\cdot W_4^{0} + p(1)\cdot W_4^{0} + p(2)\cdot W_4^{0} + p(3)\cdot W_4^{0}$
$P(1/4) = p(0)\cdot W_4^{0} + p(1)\cdot W_4^{1} + p(2)\cdot W_4^{2} + p(3)\cdot W_4^{3}$
$P(2/4) = p(0)\cdot W_4^{0} + p(1)\cdot W_4^{2} + p(2)\cdot W_4^{0} + p(3)\cdot W_4^{2}$
$P(3/4) = p(0)\cdot W_4^{0} + p(1)\cdot W_4^{3} + p(2)\cdot W_4^{2} + p(3)\cdot W_4^{2}$

○奇数

$Q(0/4) = q(0)\cdot W_4^{0} + q(1)\cdot W_4^{0} + q(2)\cdot W_4^{0} + q(3)\cdot W_4^{0}$
$Q(1/4) = q(0)\cdot W_4^{0} + q(1)\cdot W_4^{1} + q(2)\cdot W_4^{2} + q(3)\cdot W_4^{3}$
$Q(2/4) = q(0)\cdot W_4^{0} + q(1)\cdot W_4^{2} + q(2)\cdot W_4^{0} + q(3)\cdot W_4^{2}$
$Q(3/4) = q(0)\cdot W_4^{0} + q(1)\cdot W_4^{3} + q(2)\cdot W_4^{2} + q(3)\cdot W_4^{2}$

奇数はpがqになっただけです。

○は

$G(0/8) = P(0) + W_8^{0}Q(0)$
$G(1/8) = P(1) + W_8^{1}Q(1)$
$G(2/8) = P(2) + W_8^{2}Q(2)$
$G(3/8) = P(3) + W_8^{3}Q(3)$

$G(4/8) = P(0) + W_8^{4}Q(0)$
$G(5/8) = P(1) + W_8^{5}Q(1)$
$G(6/8) = P(2) + W_8^{6}Q(2)$
$G(7/8) = P(3) + W_8^{7}Q(3)$

は、180度回転することは j*j=-1 をかけることになるから

$G(0/8) = P(0) + W_8^{0}Q(0)$
$G(1/8) = P(1) + W_8^{1}Q(1)$
$G(2/8) = P(2) + W_8^{2}Q(2)$
$G(3/8) = P(3) + W_8^{3}Q(3)$

$G(4/8) = P(0) - W_8^{0}Q(0)$
$G(5/8) = P(1) - W_8^{1}Q(1)$
$G(6/8) = P(2) - W_8^{2}Q(2)$
$G(7/8) = P(3) - W_8^{3}Q(3)$

と実際は計算されます。ここで G(4) 以降の第2項に-が付いているのは
$W_N^{\frac{N}{2}+k} = -W_N^{k}$

という規則があるからです。
回転子が90度回転するのと虚数jを掛けるのは同じです。
つまり半周期(180度)回転したものはjxj=-1を掛けたものに等しいという条件です。
これらすべてを加味してグラフにすると以下のようになります。

ピンクはWが掛かった線,青線は-の線です。
これを見ると4回の掛け算+8回の加算で計算することができます。
このようなグラフをシグナルフローグラフ、
FFTのこのような演算をバタフライ演算といいます。

この図はP,Qが求まったとして話が始まっています。
今度はP,Qについて同じことをしなければなりません。

をで表す

図1-7

$A(0) = a(0)\cdot W_4^{0} + a(1)\cdot W_4^{0}$
$A(1) = a(0)\cdot W_4^{0} + a(1)\cdot W_4^{2}$
$B(0) = b(0)\cdot W_4^{0} + b(1)\cdot W_4^{0}$
$B(1) = b(0)\cdot W_4^{0} + b(1)\cdot W_4^{2}$
$C(0) = c(0)\cdot W_4^{0} + c(1)\cdot W_4^{0}$
$C(1) = c(0)\cdot W_4^{0} + c(1)\cdot W_4^{2}$
$D(0) = d(0)\cdot W_4^{0} + d(1)\cdot W_4^{0}$
$D(1) = d(0)\cdot W_4^{0} + d(1)\cdot W_4^{2}$

$P(0) = A(0) + W_4^{0}B(0)$
$P(1) = A(1) + W_4^{1}B(1)$

$P(2) = A(0) - W_4^{0}B(0)$
$P(3) = A(1) - W_4^{1}B(1)$

$Q(0) = C(0) + W_4^{0}D(0)$
$Q(1) = C(1) + W_4^{1}D(1)$

$Q(2) = C(0) - W_4^{0}D(0)$
$Q(3) = C(1) - W_4^{1}D(1)$

4点DFT(P,Q)も2点DFT で表すことができました。

をで表す

最後に1点DFTの形を作ってa(0),a(1),b(0),b(1),c(0),c(1),d(0),d(1)までで表します。

$A(0) = a(0) + W_2^{0}a(1)$
$A(1) = a(0) - W_2^{0}a(1)$

$B(0) = b(0) + W_2^{0}b(1)$
$B(1) = b(0) - W_2^{0}b(1)$

$C(0) = c(0) + W_2^{0}c(1)$
$C(1) = c(0) - W_2^{0}c(1)$

$D(0) = d(0) + W_2^{0}d(1)$
$D(1) = d(0) - W_2^{0}d(1)$

をに戻す

a(0)=f(0)
a(1)=f(4)
b(0)=f(2)
b(1)=f(6)
c(0)=f(1)
c(1)=f(5)
d(0)=f(3)
d(1)=f(7)

で戻す。

すべてを組み合わせます

これをプログラミングすることで周波数領域を求めることができます。

ビットリバース

このf(0),f(4),f(2),f(6),f(1),f(5),f(3),f(7)という並びですがこれどうやってこの順番割り出すのでしょう。

これはビットで見ると分かります。
最初私達は２つのグループを偶数と奇数に分けました。

そこで
(1) 偶数を0　奇数を1として並べてみます。

次にそのそれぞれのグループをまた2つに並べました。そこで
(2) また0,1で並べます。

(3) 最後にまたまた分割したものを並べます。

0 0 0 ->0
1 0 0 ->4
0 1 0 ->2
1 1 0 ->6
0 0 1 ->1
1 0 1 ->5
0 1 1 ->3
1 1 1 ->7

■逆に表示してみる

0 0 0->0
0 0 1->1
0 1 0->2
0 1 1->3
1 0 0->4
1 0 1->5
1 1 0->6
1 1 1->7

逆に表示してみると今度はきれいな数字の列になりました。

コンピュータでビットリバースを求めるには？

コンピュータで計算するときは

与えられたデータの最上位Bit(MSB)と最下位Bit(LSB)を交換すればよい
そして得られた番号に格納されたデータと最初のデータを交換すればよい

これはFFTをやる前とした後のどちらにも適用できますが本図ではFFTをする前の場合です。

このFFTをCooley Tukeyの時間間引きFFTと言います。

プログラムする

回転因子、ビットリバース配列

#include "common.h"
/***********************************************************************
 * NAME    : fft_table
 * FUNCTION: FFTで使用するテーブルを作成する
 * PROCESS : 省略
 *
 * INPUT   : @param double wn_FFT[] -回転因子Wnのための配列
 *         : @param short  br_FFT[] -ビット逆順のための配列
 *         : @param int    N_FFT    -FFTする全データ数
 *
 * OUTPUT  : @param double wn_FFT[] -生成された回転因子Wnのための配列
 *         : @param short  br_FFT[] -生成されたビット逆順のための配列
 *
 * RETURN  : none
 * CALL    : none
/**********************************************************************/
void fft_table(double wn_FFT[],short br_FFT[],int N_FFT)
{
   int i,n_half,ne,jp;
   double arg;
   int fft34=0;

   //calculation of twiddle factor
   arg = 2*PI / N_FFT;
    
   fft34 = (N_FFT*3)/4;
    
   for(i=0;i<fft34;i++){
      wn_FFT[i] = (double)cos(arg*i);
 
#if DEBUG
    printf("cos=%lf\tsin=%lf\n",cos(arg*i),sin(arg*i));
#endif
   } 

   //calculation of bit reversal table
   n_half = N_FFT/2;
   br_FFT[0] = 0; 

#if DEBUG
    printf("br_FFT[%d]=0x%x\n",0,br_FFT[0]);
#endif
 
 for(ne=1;ne<N_FFT;ne=ne<<1){
 
#if DEBUG
    printf("ne=%d\n",ne);
#endif

   for(jp=0;jp<ne;jp++){
     br_FFT[jp+ne] = br_FFT[jp] + n_half; 

#if DEBUG
     printf("br_FFT[%d]=0x%x\n",jp+ne,br_FFT[jp+ne]);
#endif
   } 

   n_half = n_half >> 1;
 }
}//end of fft_table

Swap部位

/***************************************************************
 * NAME    : swap
 * FUNCTION: データをひっくり返す
 * PROCESS : 省略
 * INPUT   : @param double *a -データ
 *         : @param double *b -データ
 *         
 * RETURN  : none
 * CALL    : none
**************************************************************/
void swap(double *a,double *b){
   double tmp;
   tmp = *a;
   *a = *b;
   *b = tmp;
}

FFT部位

/***************************************************************
 * NAME    : fft5
 * FUNCTION: 時間間引きFFT
 * PROCESS : 直接sin cosを計算する 三上直樹氏のプログラムを改良
 *
 * INPUT   : @param double real[]       - FFTしたいデータ
 *         : @param double imag[]       - FFTしたいデータ(全データ0を入れる)
 *         : @param short  bitreverse[] - ビット逆順配列
 *         : @param int    fftNum       - FFTしたい全データ数(2の倍数)
 *
 * OUTPUT  : @param double real[]       - FFTされたデータ
 *         : @param double imag[]       - FFTされたデータ
 *
 * RETURN  : none
 * CALL    : none
 **************************************************************/
void time_fft(double real[],double imag[],short bitreverse[],int fftNum)
{
  double xtmpr,xtmpi;

  int j,jnh,k,jxC,jxS,n_half,n_half2;
  int step;
  double arg;

  for(j=0;j<fftNum;j++){
      //条件判定ない場合　ひっくり返したあとまたひっくり返すので注意!!
      if(j<bitreverse[j]){
        swap(&real[j],&real[bitreverse[j]]);
        swap(&imag[j],&imag[bitreverse[j]]);
      }
  }

  n_half = 1;
  n_half2 = 2;

  for(step=fftNum/2;step >= 1 ;step = step/2){
       for(k=0;k<fftNum;k=k+n_half2){

#if DEBUG
            printf("k=%d\n",k);
#endif
			
          jxC = 0;
          jxS = fftNum / 4;
	    
          for(j=k;j< (k+n_half);j++){

#if DEBUG
             printf("j=%d\t",j);
#endif
            //jに対してn_half分移動したところと両者を演算する
            jnh = j + n_half;
	
            xtmpr = real[jnh];
            xtmpi = imag[jnh];
	
            arg = 2*PI / fftNum;

            real[jnh] = xtmpr*cos(arg*jxC) + xtmpi * sin(arg*jxC);
            imag[jnh] = xtmpi*cos(arg*jxC) - xtmpr * sin(arg*jxC);
		 
#if DEBUG
            printf("arg*jxC=%lf\n",arg*jxC);
#endif
            xtmpr = real[j];
            xtmpi = imag[j];

            real[j] = xtmpr + real[jnh];
            imag[j] = xtmpi + imag[jnh];
            real[jnh] = xtmpr - real[jnh];
            imag[jnh] = xtmpi - imag[jnh];
	
            //回転子は次の演算領域に入ると1/2倍ずつ増える
            jxC = jxC + step;
            //jxS = jxS + step;
         }
      }
    n_half = n_half * 2;
    n_half2 = n_half2 * 2;
  }
}

Matlabでフーリエ変換

サンプリング時間を $\Delta{t_s}$
サンプリング周波数を $\Delta{f_s} = \frac{1}{\Delta{t_s}}$
データ総数を
基本周波数 $\Delta{f} = \frac{1}{\Delta{t_s}\cdot N}$

時間信号を $f(n \Delta{t_s})$ とすると

それを離散フーリエ変換したものは、 $F(k\Delta_f)$ となる。

したがって、
$F(k\Delta_f) = F_r (k\Delta_f) + j F_i (k\Delta_f)$

振幅スペクトル:
$| F(k\Delta_f) |= \sqrt{F^2_r (k\Delta_f) + F^2_i (k\Delta_f)}$

位相スペクトル
$\angle F(k\Delta_f) = \arctan\frac{F_i (k\Delta_f)}{ F_r (k\Delta_f) }$

Matlabでは

fft

という命令だけでFFTしてくれるが使い方に注意！
以下にサンプルで示す。

用意されている窓関数

Matlab ソース

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab FFT sample program
%                 programming by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

%サンプリング周期 50[ms]
dt=50e-3; 

%サンプリング周波数20Hz 10Hz以下まで表すことができる
df=1/dt;

%離散時間（横軸）を作る
t=[0:dt:50];

%50msごとの離散時間が何個で構成されるか調べる
[gyo retu]=size(t)

%gyo行、retu列の0データを作る。
zerodata = zeros(gyo,retu);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1Hz+2Hz+5Hzの正弦波を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,1);

f1=1; %1[Hz]
f2=2; %2[Hz]
f3=5; %5[Hz]

if 1
 y1=sin(2*pi*f1*t);
else 
 y1=zerodata; %0でもうまくいく
end

if 1
 y2=sin(2*pi*f2*t);
else
 y2=0;
end

if 1
 y3=sin(2*pi*f3*t);
else
 y3=0;
end

y=y1+y2+y3;

% 1+2+5Hzの正弦波を表示する
plot(t,y,'ro-','linewidth',2);
%stem(t,y,'ro-');

hold on;

axis([0 5 -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Time[s]','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 周波数スペクトルを表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,2);
data=[];

%基本周波数の導出 基本周波数=サンプリング周波数/データ数
basicfre = df/retu

%FFTの結果は複素数である。複素数はテキストデータで保存できない
fftdata = fft(y(1:retu))';

%分解して、テキストとして保存する
fftdata_re = real(fftdata);
fftdata_im = imag(fftdata);

%離散周波数を作る
data = [data,[0:retu-1]'*df/retu];
%FFTの実部を格納する
data = [data,fftdata_re];
%FFTの虚部を格納する
data = [data,fftdata_im];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 周波数スペクトルの絶対値を求める
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%FFTの絶対値を求める
power=abs(fftdata).^2;
maxpower = max(power);

%Normalize(正規化する)
Nmpower = power/maxpower;

%Normalize + [dB/dec]
NmpowerDB = 20*log10(power/maxpower);

%正規化データを格納する
data = [data,Nmpower];
%正規化データ[dB/dec]を格納する
data = [data,NmpowerDB];

plot(data(:,1),data(:,4),'b-','linewidth',2);

axis([0 df/2  0 1]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Frequency[Hz]','Fontsize',15);
ylabel('Normalize','Fontsize',15);

save fftdata.txt data -ascii
print -djpeg fft_sample.jpg

■グラフ

■saveで保存したデータ列

 周波数          FFTの実部       FFTの虚部       FFTの絶対値     FFTの絶対値[dB/dec]
 0.0000000e+000 -8.5997875e-013 -0.0000000e+000  2.9701028e-030 -5.9054457e+002
 1.9980020e-002  1.0238257e-004  3.2621866e-002  4.2738365e-009 -1.6738364e+002
 3.9960040e-002  4.0986708e-004  6.5296747e-002  1.7123646e-008 -1.5532808e+002
 5.9940060e-002  9.2346767e-004  9.8078048e-002  3.8634755e-008 -1.4826044e+002
 7.9920080e-002  1.6448870e-003  1.3101996e-001  6.8950813e-008 -1.4322921e+002
 9.9900100e-002  2.5765350e-003  1.6417786e-001  1.0827599e-007 -1.3930936e+002
 1.1988012e-001  3.7215561e-003  1.9760876e-001  1.5687814e-007 -1.3608875e+002
 1.3986014e-001  5.0838652e-003  2.3137174e-001  2.1509300e-007 -1.3334747e+002
 1.5984016e-001  6.6681924e-003  2.6552845e-001  2.8332961e-007 -1.3095416e+002

データはこのような形で保存しておくとあとで利用しやすい。

パルス波形を生成し、それをフーリエ変換する

Matlabソース

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% パルス波形を表示
%                     programming by embedded samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

%サンプリング周期 1[ms]
dx=1e-3;
% xの範囲は -pi/2から0.001単位で2piまで
x=[-pi/2:dx:2*pi]; 
xx = x/pi;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 注意！ 
% cos(n*x) n=omega
% n=omega=1,2,..,n である。
% omega = 2*pi*f = 1,2,...,n の意味なので
% 周波数的には、
% f= 1/2pi,2/2pi,3/2pi,...,n
%  = 0.1592[Hz],0.3183[Hz],0.4775[Hz],0.6366[Hz],
%  = 0.7958[Hz],0.9549[Hz],1.1141[Hz],....,n
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1[Hz]を表示
% ここはテスト
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);

ty1=sin(2*pi*1*x);

%1[Hz]に近い
ty2=sin(6.2832*x);

plot(xx,ty1,'r-','linewidth',2);
hold on;
plot(xx,ty2,'b--','linewidth',3);

% xラベルとその文字の大きさ、線の太さの設定
str={'-p/2','0','p/2','p','3p/2','2p'}
set(gca,'FontName','symbol','xtick',[-0.5:.5:2],'xticklabel',str)
set(gca,'LineWidth',2,'FontSize',18)

% x-y範囲
axis([-Inf Inf -Inf Inf]);
grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% n=1のとき
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);

y1=0;
for m=0:1:0
 y=(4/pi)*(1/(2*m+1))*(-1)^m*cos((2*m+1).*x);
 y1=y1+y;
end

plot(xx,y1,'r-','linewidth',3);
hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% n=10のとき
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y2=0;
for m=0:1:10
 y=(4/pi)*(1/(2*m+1))*(-1)^m*cos((2*m+1).*x);
 y2=y2+y;
end

plot(xx,y2,'g-','linewidth',3);
hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% n=100のとき
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y3=0;
for m=0:1:100
 y=(4/pi)*(1/(2*m+1))*(-1)^m*cos((2*m+1).*x);
 y3=y3+y;
end

plot(xx,y3,'m-','linewidth',3);
hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 理想パルス
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data=[];

for n=x
 if (n < pi/2)
   y=1;
   data=[data,y];
 elseif (n >= pi/2) && (n <= 3*pi/2) 
   y=-1;
   data=[data,y];
 elseif ( n > 3*pi/2)
  y=1;
  data=[data,y];
 end
end

plot(xx,data,'b-','linewidth',3);
hold on;

% xラベルとその文字の大きさ、線の太さの設定
str={'-p/2','0','p/2','p','3p/2','2p'}
set(gca,'FontName','symbol','xtick',[-0.5:.5:2],'xticklabel',str)
set(gca,'LineWidth',2,'FontSize',18)

% x-y範囲
axis([-Inf Inf -Inf Inf]);
grid on

% xラベル、yラベル
xlabel('Positon t','Fontsize',20,'FontName','Century');
ylabel('pulsu wave y(t)','Fontsize',20,'FontName','Century')

% キャプション
h=legend('n=1',...,
'n=10',...,
'n=100',...,
'ideal pulsu',...,
1);
set(h,'FontSize',20,'FontName','Century');

print -djpeg pulsuwave_t.jpeg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3);

%サンプリング周波数 1kHz
df=1/dx
%xの数 7854
[gyo retu]=size(x)

%基本周波数 0.1273[Hz]
bf=df/retu

%FFTの結果は複素数である
fftdata1 = fft(y1)';
fftdata2 = fft(y2)';
fftdata3 = fft(y3)';
fftdata4 = fft(data)';

%FFTの絶対値を求める
power1=abs(fftdata1).^2;
power2=abs(fftdata2).^2;
power3=abs(fftdata3).^2;
power4=abs(fftdata4).^2;

%離散周波数を作る
freq = [0:retu-1]'*bf

plot(freq,power1,'r-','linewidth',3);
hold on
plot(freq,power2,'g-','linewidth',3);
hold on
plot(freq,power3,'m+--','linewidth',2);
hold on
plot(freq,power4,'bo-','linewidth',1);
hold on

axis([0 1  0 Inf]);

set(gca,'LineWidth',2,...,
'FontSize',15);

xlabel('Frequency[Hz]','Fontsize',20,'FontName','Century');
ylabel('Normalize','Fontsize',20,'FontName','Century');

% キャプション
h=legend('n=1',...,
'n=10',...,
'n=100',...,
'ideal pulsu',...,
1);
set(h,'FontSize',20,'FontName','Century');

print -djpeg pulsuwave_f.jpg

/***********************************************************/
/* FFT Test Program by embedded.samurai                    */
/*                                                         */
/***********************************************************/

#include <stdio.h>
#include <math.h>
#include "common.h"

#define INPUT_RECTANGLE 0

int main(void)
{
   short i;
   double wn_FFT[nFFT3_4]; //for twiddle factor
   short  br_FFT[nFFT];    //for bit reversal
   double xr[nFFT],xi[nFFT];
   double yr[nFFT],yi[nFFT];
   double arg;
	
   FILE *fin,*fout;

#if INPUT_RECTANGLE
    //波形を作る
    rectangle(xr,nFFT);
#else
     
    fin = fopen("sin20hz.txt","r");
    //入力データを配列に格納する
    i=0;
    while(1){
      float buf;
      fscanf(fin,"%f",&buf);

      if(ferror(fin)){
          printf("read error\n");
          return(-1);
      }
 		
      if(feof(fin)){
           break;
       }
		
       xr[i++]= buf;
       if(i>=nFFT){
          fclose(fin);
          break;
       }
    }//end of while(1)
#endif

   //出力
   fout = fopen("output_vc.txt","w");

   //テーブルを作る generate tables for FFT
   fft_table(wn_FFT,br_FFT,nFFT);
    
   //出力データの配列を作る
   for(i=0;i<nFFT;i++){
      yr[i] = xr[i];
      yi[i] = xi[i] = 0;
   }
   //FFTする
   fft2(yr,yi,br_FFT,nFFT);
    
   //FFT4に必要
   //arg = 2*PI / nFFT; 
   //fft4(nFFT, arg, yr, yi);

   //IFFTする
   i_fft1(yr,yi,br_FFT,nFFT);

   //結果を出力ファイルに格納する
   for(i=0;i<nFFT;i++){
     fprintf(fout,"%lf\t%lf\n",yr[i],yi[i]);
   }
  
   fclose(fout);
   return 0;
}

FFT

/*============================================================*/
/*        FFT programming by embedded.samurai                 */
/*                                                            */
/*============================================================*/
#include "common.h"
/***********************************************************************
 * NAME    : fft_table
 * FUNCTION: FFTで使用するテーブルを作成する
 * PROCESS : 省略
 *
 * INPUT   : @param double wn_FFT[] -回転因子Wnのための配列
 *         : @param short  br_FFT[] -ビット逆順のための配列
 *         : @param int    N_FFT    -FFTする全データ数
 *
 * OUTPUT  : @param double wn_FFT[] -生成された回転因子Wnのための配列
 *         : @param short  br_FFT[] -生成されたビット逆順のための配列
 *
 * RETURN  : none
 * CALL    : none
/**********************************************************************/
void fft_table(double wn_FFT[],short br_FFT[],int N_FFT)
{
   int i,n_half,ne,jp;
   double arg;
   int fft34=0;

   //calculation of twiddle factor
   arg = 2*PI / N_FFT;
    
   fft34 = (N_FFT*3)/4;
    
   for(i=0;i<fft34;i++){
      wn_FFT[i] = (double)cos(arg*i);
 
#if DEBUG
    printf("cos=%lf\tsin=%lf\n",cos(arg*i),sin(arg*i));
#endif
   }

   //calculation of bit reversal table
   n_half = N_FFT/2;
   br_FFT[0] = 0;

#if DEBUG
    printf("br_FFT[%d]=0x%x\n",0,br_FFT[0]);
#endif
 
 for(ne=1;ne<N_FFT;ne=ne<<1){
 
#if DEBUG
    printf("ne=%d\n",ne);
#endif

   for(jp=0;jp<ne;jp++){
     br_FFT[jp+ne] = br_FFT[jp] + n_half;

#if DEBUG
     printf("br_FFT[%d]=0x%x\n",jp+ne,br_FFT[jp+ne]);
#endif
   }

   n_half = n_half >> 1;
 }
}//end of fft_table

/***************************************************************
 * NAME    : swap
 * FUNCTION: データをひっくり返す
 * PROCESS : 省略
 * INPUT   : @param double *a -データ
 *         : @param double *b -データ
 *         
 * RETURN  : none
 * CALL    : none
**************************************************************/
void swap(double *a,double *b){
   double tmp;
   tmp = *a;
   *a = *b;
   *b = tmp;
}

/***************************************************************
 * NAME    : fft1
 * FUNCTION: テーブルを使ったFFT
 * PROCESS : テーブルを使ったFFT 三上直樹氏のプログラム
 *
 * INPUT   : @param double xR[]     - FFTしたいデータ
 *         : @param short  xI[]     - FFTしたいデータ(全データ0を入れる)
 *         : @param double wn_FFT[] - 回転因子Wn配列
 *         : @param short  br_FFT[] - ビット逆順配列
 *         : @param int    N_FFT    - FFTしたい全データ数(2の倍数)
 *
 * OUTPUT  : @param double xR[]     - FFTされたデータ
 *         : @param short  xI[]     - FFTされたデータ
 *
 * RETURN  : none
 * CALL    : none
 **************************************************************/
void fft1(double xR[],double xI[],double wn_FFT[],short br_FFT[],int N_FFT)
{
   double xtmpR,xtmpI;
   int j,jnh,k,jxC,jxS,ne,n_half,n_half2;
   
   n_half = N_FFT / 2;

   for(ne=1;ne<N_FFT;ne = ne<<1){
       n_half2 = n_half *2;
#if DEBUG
       printf("n_half2=%d\n",n_half2); 
#endif
       for(k=0;k<N_FFT;k=k+n_half2){
#if DEBUG
           printf("k=%d\n",k);
#endif
           jxC = 0;
           jxS = N_FFT / 4;

           for(j=k;j< (k+n_half);j++){

#if DEBUG
               printf("j=%d\t",j);
#endif
               jnh = j + n_half;
               //beginning of butterfly operations
               xtmpR = xR[j];
               xtmpI = xI[j];
               xR[j] = xtmpR + xR[jnh];
               xI[j] = xtmpI + xI[jnh];
               xtmpR = xtmpR - xR[jnh];
               xtmpI = xtmpI - xI[jnh];
  
#if DEBUG
               printf("jnh=%d,jxC=%d,jxS=%d\n",jnh,jxC,jxS);
#endif
               xR[jnh] = xtmpR*wn_FFT[jxC] - xtmpI * wn_FFT[jxS];
               xI[jnh] = xtmpR*wn_FFT[jxS] + xtmpI * wn_FFT[jxC];
      
               //end of butterfly operations
               jxC = jxC + ne;
               jxS = jxS + ne;
           }//end of for(j=k;j< (k+n_half);j++)
       }//end of for(k=0;k<N_FFT;k=k+n_half2){
	
    n_half = n_half >> 1;
    
  }//end of for(ne=1;ne<N_FFT;ne = ne<<1)
    
  //Bit reverse
  //Bit reverse時のデータ取替えに注意
  for(j=0;j<N_FFT;j++){
   if(j<br_FFT[j]){
     swap(&xR[j],&xR[br_FFT[j]]);
     swap(&xI[j],&xI[br_FFT[j]]);
   }
  }
}


/***************************************************************
 * NAME    : fft2
 * FUNCTION: テーブルを使用しないFFT
 * PROCESS : 直接sin cosを計算する 三上直樹氏のプログラムを改良
 *
 * INPUT   : @param double real[]       - FFTしたいデータ
 *         : @param double imag[]       - FFTしたいデータ(全データ0を入れる)
 *         : @param short  bitreverse[] - ビット逆順配列
 *         : @param int    fftNum       - FFTしたい全データ数(2の倍数)
 *
 * OUTPUT  : @param double real[]       - FFTされたデータ
 *         : @param double imag[]       - FFTされたデータ
 *
 * RETURN  : none
 * CALL    : none
 **************************************************************/
void fft2(double real[],double imag[],short bitreverse[],int fftNum)
{
   double xtmpr,xtmpi;

   int j,jnh,k,jxC,jxS,n_half,n_half2;
   int step;
   double arg;

   n_half = fftNum / 2;

#if DEBUG
   for(j=0;j<fftNum;j++){
     printf("real[%d]=%lf\n",j,real[j]);
   }
#endif

   for(step=1;step < fftNum;step = step*2){
      n_half2 = n_half *2;
#if DEBUG
      printf("n_half2=%d\n",n_half2);
#endif
	
      for(k=0;k<fftNum;k=k+n_half2){
#if DEBUG
          printf("k=%d\n",k);
#endif
          jxC = 0;
          jxS = fftNum / 4;

          for(j=k;j< (k+n_half);j++){
#if DEBUG
             printf("j=%d\t",j);
#endif
             //jに対してn_half分移動したところと両者を演算する
             jnh = j + n_half;
 			
             xtmpr = real[j];
             xtmpi = imag[j];
		
             //バタフライ演算の上側には回転子の項がないため楽
             real[j] = xtmpr + real[jnh];
             imag[j] = xtmpi + imag[jnh];
		
             //バタフライ演算の下の項には回転子の項があるため注意!!
             xtmpr = xtmpr - real[jnh];
             xtmpi = xtmpi - imag[jnh];
		
             /* 
                (A+jB)*(cosx+jsinx)
                 Acosx + jAsinx +jBcosx-Bsinx
                 realpart-> Acosx - Bsinx
                 imagpart-> Asinx + Bcosx
             */

             arg = 2*PI / fftNum;
             //real[jnh] = xtmpr*cos(arg*jxC) - xtmpi * sin(arg*jxC);
             //imag[jnh] = xtmpr*sin(arg*jxC) + xtmpi * cos(arg*jxC);

             real[jnh] = xtmpr*cos(arg*jxC) + xtmpi * sin(arg*jxC);
             imag[jnh] = xtmpi*cos(arg*jxC) - xtmpr * sin(arg*jxC);

             printf("jnh=%d,jxC=%d,jxS=%d\n",jnh,jxC,jxS);
             printf("jnh=%d,jxC=%d,jxS=%d\n",jnh,jxC,jxS);
             
             //回転子は次の演算領域に入ると２倍ずつ増える
             jxC = jxC + step;
             jxS = jxS + step;
         }
       }
       n_half = n_half / 2;
     }

  //Bit reverse 
  for(j=0;j<fftNum;j++){
    //条件判定ない場合　ひっくり返したあとまたひっくり返すので注意!!
    if(j<bitreverse[j]){
       swap((double *)&real[j],(double *)&real[bitreverse[j]]);
       swap((double *)&imag[j],(double *)&imag[bitreverse[j]]);
    }
  }
}

/************************************************************************
 * NAME    : fft4
 * FUNCTION: 周波数間引きFFT　programming by 大浦拓哉
 * PROCESS : このプログラムは，"FFT (高速フーリエ・コサイン・サイン変換)
 *           の概略と設計法"のページ :
 *           http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/
 *           のCソースです。
 *
 * INPUT   : @param int n             - FFTしたい全データ数(2の倍数)
 *         : @param double theta      - 
 *         : @param double ar[]       - FFTしたいデータ
 *         : @param double ai[]       - FFTしたいデータ(全データ0を入れる)
 *
 * OUTPUT  : @param double ar[]       - FFTされたデータ
 *         : @param double ai[]       - FFTされたデータ
 *
 * RETURN  : none
 * CALL    : none
*************************************************************************/
void fft4(int n, double theta, double ar[], double ai[])
{
   int m, mh, i, j, k;
   double wr, wi, xr, xi;

   /* ---- scrambler ---- */
   i = 0;
   for (j = 1; j < n - 1; j++) {
       for (k = n >> 1; k > (i ^= k); k >>= 1);
       if (j < i) {
           xr = ar[j];
           xi = ai[j];
           ar[j] = ar[i];
           ai[j] = ai[i];
           ar[i] = xr;
           ai[i] = xi;
       }
   }
   theta *= n;

#if DEBUG
    printf("theta=%lf\n",theta);
#endif
 
   for (mh = 1; (m = mh << 1) <= n; mh = m) {
 
#if DEBUG
   printf("mh=%d\n",mh);
   printf("m=%d\n",m);
#endif
 
       theta *= 0.5;
 
#if DEBUG
   printf("theta=%lf\n",theta);
#endif

       for (i = 0; i < mh; i++) {

           wr = cos(theta * i);
           wi = sin(theta * i);

#if DEBUG
           printf("theta*i=%lf\n",theta*i);
#endif

           for (j = i; j < n; j += m) {
               k = j + mh;

#if DEBUG
               printf("j=%d\tk=%d\n",j,k);
#endif
                
               xr = wr * ar[k] - wi * ai[k];
               xi = wr * ai[k] + wi * ar[k];
		
               ar[k] = ar[j] - xr;
               ai[k] = ai[j] - xi;

               ar[j] += xr;
               ai[j] += xi;
           }
       }
   }
}
/***************************************************************
 * NAME    : fft5
 * FUNCTION: 時間間引きFFT
 * PROCESS : 直接sin cosを計算する 三上直樹氏のプログラムを改良
 *
 * INPUT   : @param double real[]       - FFTしたいデータ
 *         : @param double imag[]       - FFTしたいデータ(全データ0を入れる)
 *         : @param short  bitreverse[] - ビット逆順配列
 *         : @param int    fftNum       - FFTしたい全データ数(2の倍数)
 *
 * OUTPUT  : @param double real[]       - FFTされたデータ
 *         : @param double imag[]       - FFTされたデータ
 *
 * RETURN  : none
 * CALL    : none
 **************************************************************/
void fft5(double real[],double imag[],short bitreverse[],int fftNum)
{
   double xtmpr,xtmpi;

   int j,jnh,k,jxC,jxS,n_half,n_half2;
   int step;
   double arg;

   for(j=0;j<fftNum;j++){
       //条件判定ない場合　ひっくり返したあとまたひっくり返すので注意!!
       if(j<bitreverse[j]){
         swap(&real[j],&real[bitreverse[j]]);
         swap(&imag[j],&imag[bitreverse[j]]);
       }
   }

   n_half = 1;
   n_half2 = 2;

   for(step=fftNum/2;step >= 1 ;step = step/2){
        for(k=0;k<fftNum;k=k+n_half2){

#if DEBUG
            printf("k=%d\n",k);
#endif
			
           jxC = 0;
           jxS = fftNum / 4;
	    
           for(j=k;j< (k+n_half);j++){

#if DEBUG
             printf("j=%d\t",j);
#endif
             //jに対してn_half分移動したところと両者を演算する
             jnh = j + n_half;
	
             xtmpr = real[jnh];
             xtmpi = imag[jnh];
	
             arg = 2*PI / fftNum;

             real[jnh] = xtmpr*cos(arg*jxC) + xtmpi * sin(arg*jxC);
             imag[jnh] = xtmpi*cos(arg*jxC) - xtmpr * sin(arg*jxC);
		 
#if DEBUG
             printf("arg*jxC=%lf\n",arg*jxC);
#endif
             xtmpr = real[j];
             xtmpi = imag[j];

             real[j] = xtmpr + real[jnh];
             imag[j] = xtmpi + imag[jnh];
             real[jnh] = xtmpr - real[jnh];
             imag[jnh] = xtmpi - imag[jnh];
	
             //回転子は次の演算領域に入ると1/2倍ずつ増える
             jxC = jxC + step;
             //jxS = jxS + step;
          }
       }
     n_half = n_half * 2;
     n_half2 = n_half2 * 2;
   }
}

IFFT フーリエ逆変換

#include "common.h"

/***************************************************************
 * NAME    : i_fft1
 * FUNCTION: 逆FFT
 * PROCESS : 直接sin cosを計算する
 *
 * INPUT   : @param double real[]       - FFTされたデータ
 *         : @param double imag[]       - FFTされたデータ
 *         : @param short  bitreverse[] - ビット逆順配列
 *         : @param int    fftNum       - 逆FFTしたい全データ数(2の倍数)
 *
 * OUTPUT  : @param double real[]       - 逆FFTされたデータ
 *         : @param double imag[]       - 逆FFTされたデータ
 *
 * RETURN  : none
 * CALL    : none
 **************************************************************/
void i_fft1(double real[],double imag[],short bitreverse[],int fftNum)
{
   double xtmpr,xtmpi;

   int j,jnh,k,jxC,jxS,n_half,n_half2;
   int step;
   double arg;

   n_half = fftNum / 2;
    
#if DEBUG
   for(j=0;j<fftNum;j++){
     printf("real[%d]=%lf\n",j,real[j]);
   }
#endif

   for(step=1;step < fftNum;step = step*2){

       n_half2 = n_half *2;

#if DEBUG
       printf("n_half2=%d\n",n_half2);
#endif

       for(k=0;k<fftNum;k=k+n_half2){

#if DEBUG
           printf("k=%d\n",k);
#endif
		
           jxC = 0;
           jxS = fftNum / 4;

           for(j=k;j< (k+n_half);j++){

#if DEBUG
		printf("j=%d\t",j);
#endif

              //jに対してn_half分移動したところと両者を演算する
              jnh = j + n_half;
			
              xtmpr = real[j];
              xtmpi = imag[j];
		
              //バタフライ演算の上側には回転子の項がないため楽
              real[j] = xtmpr + real[jnh];
              imag[j] = xtmpi + imag[jnh];

              //バタフライ演算の下の項には回転子の項があるため注意!!
              xtmpr = xtmpr - real[jnh];
              xtmpi = xtmpi - imag[jnh];
		
              /* (A+jB)*(cosx+jsinx)
                 Acosx + jAsinx +jBcosx-Bsinx
                 realpart-> Acosx - Bsinx
                 imagpart-> Asinx + Bcosx
              */

              arg = 2*PI / fftNum;
              real[jnh] = xtmpr*cos(arg*jxC) - xtmpi * sin(arg*jxC);
              imag[jnh] = xtmpr*sin(arg*jxC) + xtmpi * cos(arg*jxC);
              //real[jnh] = xtmpr*cos(arg*jxC) + xtmpi * sin(arg*jxC);
              //imag[jnh] = xtmpi*cos(arg*jxC) - xtmpr * sin(arg*jxC);

#if DEBUG
               printf("jnh=%d,jxC=%d,jxS=%d\n",jnh,jxC,jxS);
               printf("jnh=%d,jxC=%d,jxS=%d\n",jnh,jxC,jxS);
#endif

              //end of butterfly operations
              //回転子は次の演算領域に入ると２倍ずつ増える
              jxC = jxC + step;
              jxS = jxS + step;
           }
        }
        n_half = n_half / 2;
   }

  //Bit reverse 
  for(j=0;j<fftNum;j++){
    //条件判定ない場合　ひっくり返したあとまたひっくり返すので注意!!
    if(j<bitreverse[j]){
        swap(&real[j],&real[bitreverse[j]]);
        swap(&imag[j],&imag[bitreverse[j]]);
    }
  }

  for(j=0;j<fftNum;j++){
      real[j] = real[j] / fftNum;
      imag[j] = imag[j] / fftNum;
  }
}