FFT

フーリエ変換は、FFT (Fast Fourier Transform) を用いて、文字どおり、高速に数値計算できます。

FFT は、サンプル点を 2の累乗に選ぶ必要があります。FFT の概念は、1965年に考案されましたが、当時は、まだ計算機の能力が低く、実用的に使用されるようになったのは、それから 20年程度経過してからのことでした。とはいっても、その頃は、高速に処理するには専用 LSIを用いる必要があり、筆者の経験では、ソフトウェアで処理するには、ワークステーションに高速演算用のハードウェア・オプションを追加しても、数分を要したと記憶しています。

エクセルによる FFT

今では、エクセルに組み込まれており、それほど高速ではない PC でも瞬時に計算できるようになりました。

時間関数を周波数関数に変換するのが FFT で、その逆は、iFFT (inverse FFT) といいます。エクセルのアドインは、バッチ処理となりますが、FFT のアルゴリズムは簡単に理解できるので、自分で組み立てると、比較的簡単にリアルタイム処理が可能となります。

ここでは、アドインを用いる方法について説明します。例題として、「ラプラス変換とフーリエ変換」に記載したパルス波形を用います。

Article header 128345 sc57 fig1  1
図1 定数の定義

まず、図1 のように定数を定義します。例では、周期 T=10ns、パルス幅 TW=3ns、立ち上がり時間 tr=1ns です。角周波数 ω のキザミ(ステップ)を B8 に求めておきます。キザミの周波数は 1/T なので、Δω=2π/T です。iFFT の処理をしたときの時間のキザミを B9 に入力します。キザミは、周期を FFT のサンプル数 N で割る、すなわち Δt=T/N です。今回の例では N=64 としました。

  • π::ギリシャ文字小文字のパイ
  • ω:ギリシャ文字小文字のオメガ
Article header 128345 sc57 fig2  1
図2 F(ω) の各要素

図2 は、F(ω) の各要素です。D列は ω で、0 から Δω キザミです。行は全部で 64 あります。文字色がグレーのセルは、その上のセルを下向きにコピーします。これ以降の列についても同様です。E列は、ωTW/2、F列は ωtr/2 です。G列は ωTW/2 の sinc 関数、すなわち、sin(x)/x です。同様に H列は ωtr/2 の sinc 関数です。sinc(0)=1 なので、ω=0 の行は、1 とします。I列は exp{-jω(TW+tr)/2} です。複素数の指数関数は、imexp です。D列から I列は、2~34 行まで 33 個計算します。N/2+1 です。I35 以下は折り返しの計算のためのインデックスです。

Article header 128345 sc57 fig3  1
図3 F(ω)

図3 は F(ω) です。全体の係数は、B4 の TW/T で、二つの sinc 関数(G列と H列)をかけて、complex 関数で複素数化し、improduct 関数で指数部(I列)をかけます。(J列)F(ω) は、正負の ω を求めます。J2 が直流で、J3 から J34 までが正の ω です。負の ω は、正の ω の F(ω) の共役複素数なので、J35 を中心に対称に共役複素数を求めます。そのためのインデックスが I35 から下です。J35 は、offset($J$35, I35) で、J35 の -1、すなわち J34 をみつけて、その共役複素数を imconjugate 関数で求めます。下に向かって、offset の負の数がひとつづつ増えていくので、J36 は J35 の -2、すなわち J33 をみつけて共役複素数を求めます。以下同様に、J65 まで求めます。

F(ω) が求まったので、iFFT(F(ω) から f(t)) を計算します。

エクセルの準備(アドイン)

FFT/iFFT を使うには、エクセルにフーリエ解析のアドインを追加する必要があります。その方法は、エクセルのバージョンによって異なります。

古いバージョン(2003 以前)
[準備]

  • ツール → アドイン → 分析ツール にチェック

[解析]

  • ツール → 分析ツール → フーリエ解析


新しいバージョン(2007 以降)
[準備]

  • ファイル → オプション → アドイン → 管理 Excel アドインの設定をクリック → 分析ツール にチェック

[解析]

  • データ → データ分析 → フーリエ解析

フーリエ解析

いずれのバージョンも上記「フーリエ解析」以降はほぼ同じです。
入力元 $J$2:$J$65
出力先 $K$2
逆変換にチェック

これで、K2:K65 に時間関数(虚数部が 0 の複素数)が出力されます。グラフ化するために K列の実数部を M列に取り出します。時間軸は、L列で、キザミは B9 の Δt です。

Article header 128345 sc57 fig4  1
図4 計算結果
Article header 128345 sc57 fig5  1
図5 計算結果(N=1024)

図4 は、時間を L2:L65、振幅を M2:M65 でグラフ化したものです。図1 の台形波が再現できていることがわかります。この例は、簡単のために N=64 の FFT で、図4 はわずかなリンギングが生じています。より細かく表示するには、N=1024 程度にしてみると、図5 のようにきれいな波形が再現できました。

今回は、エクセルの FFT について紹介しました。せっかくの機能なので、ぜひ活用してください。

碓井有三のスペシャリストコラムとは?

基礎の基礎といったレベルから入って、いまさら聞けないようなテーマや初心者向けのテーマ、さらには少し高級なレベルまでを含め、できる限り分かりやすく噛み砕いて述べている連載コラムです。

もしかしたら、他にも気になるテーマがあるかも知れませんよ!