« 能登半島地震 記録と雑感 | トップページ | ホウ素とフッ素の排水基準、温泉旅館への適用は延期へ »

2007年3月29日 (木)

高速フーリエ変換 Excel VBA用FFTプログラム

ちまたには高速フーリエ変換(FFT)のプログラムが出回っているが、これをExcel VBA(いわゆるマクロ)で使えるようにしたところ、机上作業の効率が飛躍的にアップした。下記にプログラム等をノートしておく。ちまたに出回っているFFTのプログラムは利用者の知らないうちに窓関数がかけられていたり、周波数を自動的に計算したりするので、あまり教育的ではない。一番シンプルなFFTのプログラムで、生の処理を体験すべきである。

ちなみに、フーリエ変換とは、理工系の大学で習う数学的処理の方法の一つで端的に言えば、横軸が時間となっているグラフをフーリエ変換すると、横軸が周波数のグラフとなるもの。

Sin_t 例えば左図は横軸が時間で、周波数の異なるサインカーブを3つ足し合わせたグラフである。0.001秒毎に(サンプリング周波数1 kHz)、1024個のデータを取得したもの。横軸が時間だと、この波の周波数を読み取るのに一苦労である。

Sin_f これをフーリエ変換すると左のようになる。横軸が周波数となって、何Hzの波が含まれているかがクリアにわかる。時間領域と周波数領域の間を互いに変換するのがフーリエ変換である。

このフーリエ変換のプログラムを、CooleyとTukeyの発見による感動的なアルゴリズムで高速化したのが高速フーリエ変換(Fast Fourier Transform, FFT)である。「高速」のつかないフーリエ変換のプログラムもあるにはあるが普通は用いられない(データ数が2の整数乗個でなくてもよいという点だけがメリット)。

注意事項

  • 高速フーリエ変換を行う元のデータも、変換後のデータもすべて複素数。実際には元のデータが複素数であることはまずないので、虚数部分を0とする。変換後のデータは、複素数の絶対値を求めて使用する。
  • データの数は2の整数乗個でなければならない。例 1024, 2048, 4096
  • 変換後の周波数領域のデータは左右対称となるので左半分を使用する。右半分は削除する。横軸の右端はサンプリング周波数の半分となる。
  • サンプリング周波数の半分以上の周波数成分が信号に含まれていたとしても、観察されるはずがないことを直観で理解して欲しい。同じようにサンプリングした時間よりも周期が大きい波、つまり時間領域のグラフに波1個さえ入れない低周波成分も出てくるはずがない。特に後者は要注意で上の周波数領域のグラフで原点に向かって立ち上がりがあるが、別に低周波成分が含まれていることを示すわけではない。サンプリングしている間に、波が複数入ってないと周期的変動として検知できないのである。

使用法

  1. プログラム中にデータ数を指定する(5行目のnの部分と7行目の括弧の中)。下のプログラム例ではデータ数は1024個になっている。データの個数は2の整数乗でなければならない。
  2. ExcelのA1~A1024に元データの実数部分、B1~B1024に虚数部分を入力する。ただし、虚数部分があるデータはまれだと思われるので、B1~B1024にはゼロを入力しておけばよい。実数部分で、データ数が2の整数乗に足りない場合は余ったところにゼロを入力しておく人もいれば、同じデータを繰り返しておく人もいる(厳密にはよろしくないと思われる)。
  3. Excelの[表示]→[ツールバー]→[Visual Basic]で、VBAのメニューバーを表示、Visual Basic Editorのアイコンをクリックする。Visual Basic Editor上で[挿入]→[標準モジュール]、さらに[挿入]→[プロシージャ]でSub、Publicを選択し、名前を適当に入力して、下記プログラムをコピー&ペースト(最初と最後の1行は重複するので削除)
  4. メニューバーからプログラムを実行する(▲が横向きのアイコン)。
  5. 4列目(Dの列)にフーリエ変換結果の実数部分、5列目に虚数部分、6列目に(複素数の)絶対値が表示される。
  6. 普通は6列目を使用する。グラフを描くと、左右対称になっていることがわかるが、左半分のみを使用する。横軸の左端は0 Hzとなるが、左半分の右端(左右対称のグラフだと中央、512個目のデータ)はサンプリング周波数の半分に対応する。例えば、元のデータが0.001秒毎に取得したデータであれば、サンプリング周波数は1000 Hzなので、グラフの右端は500Hzとなる。この500Hz、つまりサンプリング周波数の半分のことをナイキスト周波数と呼ぶこともある。

Public Sub fft()
Dim g, h, i, j, k, l, m, n, o, p, q As Integer
i = 0: j = 0: k = 0: l = 0: p = 0: h = 0: g = 0: q = 0
'データ数nを指定(2の整数乗である必要あり)
n = 1024: m = Log(n) / Log(2)
'xr,xiはデータ数以上、s,cはデータ数の半分以上を指定しておく
Dim xr(1024), xi(1024), xd, s(512), c(512), a, b As Single 
'データの読み込み 1列目を実数部、2列目を虚数部とする
For i = 1 To n
xr(i - 1) = Cells(i, 1): xi(i - 1) = Cells(i, 2)
Next i
'FFTの計算
a = 0: b = 3.14159265359 * 2 / n
For i = 0 To n / 2
s(i) = Sin(a): c(i) = Cos(a): a = a + b
Next i
l = n: h = 1
For g = 1 To m: l = l / 2: k = 0
For q = 1 To h: p = 0
  For i = k To l + k - 1: j = i + l
   a = xr(i) - xr(j): b = xi(i) - xi(j)
   xr(i) = xr(i) + xr(j): xi(i) = xi(i) + xi(j)
   If p = 0 Then
    xr(j) = a: xi(j) = b
   Else
    xr(j) = a * c(p) + b * s(p): xi(j) = b * c(p) - a * s(p)
   End If
   p = p + h: Next i
k = k + l + l: Next q
h = h + h: Next g
j = n / 2
For i = 1 To n - 1: k = n
If j < i Then
xd = xr(i): xr(i) = xr(j): xr(j) = xd: xd = xi(i): xi(i) = xi(j): xi(j) = xd
End If
k = k / 2
Do While j >= k
j = j - k: k = k / 2
Loop
j = j + k
Next i
'データの出力
For i = 1 To n   '4列目に実数部、5列目に虚数部、6列目に絶対値を表示
Cells(i, 4) = xr(i - 1): Cells(i, 5) = xi(i - 1)
Cells(i, 6) = (xr(i - 1) ^ 2 + xi(i - 1) ^ 2) ^ 0.5
Next i
End Sub

先人の業績をExcel VBAで使えるようにしたものなので、ご自由にご利用いただいて構いません。

参考文献

科学計測のための波形データ処理(南茂夫編著、CQ出版) amazonへリンク 古書としてのみ入手可能・・・本書は名著中の名著。FFTに関しても面白くてためになる。しかし、PCの方が進歩しすぎたためか、残念ながら絶版。

ニューメリカルレシピ・イン・シー C言語による数値計算のレシピ(技術評論社) amazonへリンク・・・これも名著。内容が充実している。FFTでは実際にはデータが実数で、虚数部分のゼロが計算の無駄になるので、計算を工夫して、2つ分のフーリエ変換を同時に行うプログラムが同書に提案されている。

補足 逆フーリエ変換(周波数領域→時間領域)について

逆フーリエ変換は、上と同じプログラムで実行できる。ただし、注意点として、元となる周波数領域のデータ(複素数)は、左右対称のデータを用いる必要がある(フーリエ変換直後のデータと同様、右半分と左半分が対称なものを作成する)。また、縦軸の倍率が変わるので、縦軸の数値を全データ数(1024とか4096とか)で割った上、正負を逆にする必要がある。

時間領域のデータから出発して、フーリエ変換→逆フーリエ変換を行って、元のデータと比較すれば理解できる。フーリエ変換→逆フーリエ変換を行った後、時間領域のデータの虚数成分は完全なゼロにはならず、ゼロに近い数値になるが、これは誤差のためである。

補足 C言語のFFTについて

上述の技術評論社のニューメリカルレシピ・イン・シー(C言語による数値計算のレシピ)のp.390に記載されています。

|

« 能登半島地震 記録と雑感 | トップページ | ホウ素とフッ素の排水基準、温泉旅館への適用は延期へ »

コメント

FFTのプログラム、とても参考になりました。
ありがとうございます。

ところで、一点教えていただきたいことがあります。
ハードウェアがまだできていないので、波形はまだ
見ていないのですが、FFTで周波数解析を行い、
低周波と高周波のノイズ除去を行いたいのです。

このような場合、ディジタルフィルタを用いるのが
最適でしょうか?別の方法をご存知であれば
アドバイスをいただきたく、よろしくお願いします。

投稿: あんちゃん | 2007年8月 9日 (木) 16時38分

コメントありがとうございます。
他にはローパスフィルタやハイパスフィルタを用いて、ハード的にノイズ除去を行うこともできますが、どちらが最適かはケースバイケースと思います。

投稿: tsuyumoto | 2007年8月 9日 (木) 21時06分

アドバイス、ありがとうございました。

甘えついでに、もう一点お願いがあります。
逆FFTプログラムの掲載をしていただけると、
非常にありがたいです。

ご検討、よろしくお願い申し上げます。

投稿: あんちゃん | 2007年8月10日 (金) 09時47分

補足として加筆しました。同じプログラムで可能です。ただし、周波数領域のデータは左右対称なものを用い、縦軸の倍率に注意してください。

投稿: tsuyumoto | 2007年8月10日 (金) 12時52分

返事おそくなってもうしわけございません。
補足文書、ありがとうございました。
とても参考になります。

今後ともよろしくお願い申し上げます。

投稿: あんちゃん | 2007年8月20日 (月) 14時09分

始めまして。
56歳の男性です。
FFTのプログラムで、サンプリングの周波数は、どこに入力するのでしょうか、教えてください。

投稿: jiisan | 2007年8月25日 (土) 19時31分

サンプリングの周波数は入力する必要はありません。処理後の横軸の最大値が、サンプリング周波数の半分になりますので。

投稿: tsuyumoto | 2007年8月25日 (土) 20時29分

早速の回答ありがとうございます。
質問の内容が不十分でした。
改めて、質問します。
紹介されているFFTのプログラムを応用する場合、
(例えば、サンプリングタイム0.02sec,サンプリング数1024の場合)
どのように変更すれば、よいのでしょうか。

宜しく、お願いします。

投稿: jiisan | 2007年8月25日 (土) 23時01分

プログラムを変更する必要はありません。
サンプリングタイム0.02 secならサンプリング周波数が50 Hzだと思いますので、周波数領域(変換後データ)の最大値が25 Hzとなります。

投稿: tsuyumoto | 2007年8月26日 (日) 11時06分

ありがとうございまあした。

投稿: jiisan | 2007年8月26日 (日) 14時31分

うちまちがいをしました。
改めて、FFTプログラムを使用させていただきます。
ありがとうございました。

投稿: jiisan | 2007年8月26日 (日) 14時34分

これから仕事で使おうかと思っていたところで, この記事を見つけ大変参考になりました.

瑣末なことで恐縮ですが1点気になったことがあります.

dim i,j as Integer

のような記述がありますが, Excel VBA (2003) では j は Integer になるものの, i は Variant になってしまうのではないかと思います.
(最新版の Excel は持っていないので確認できません. 申し訳ありません.)

プログラムの動作には影響ないかと思いますが,
私自身過去にこれではまったことがありましたので, コメントさせていただきました.

今後の記事も楽しみにしています.

投稿: ぼやっと | 2007年9月 5日 (水) 12時54分

ご指摘ありがとうございます。これは知りませんでした。

確かに

dim i as Integer

で、i を整数型に定義しても、 i はVariantのままになってます。
本当に盲点でした。

今後のプログラミングには気をつけます。

投稿: tsuyumoto | 2007年9月 6日 (木) 00時01分

VBのプロシージャ作成にとても参考になりました。
ひとつ発見がありましたの報告いたします。
先の投稿で dim i,j as Integer に関してコメントがあり、
その部分については、動作に影響ありませんでしたが、
Dim g, h, i, j, k, l, m, n, o, p, q As Integer 
の箇所については、一部影響する箇所がございます。
このまま使用していますと問題ないのですが、
すべて Integer 宣言しますと、kが Integer の場合
 For i = 1 To n - 1: k = n
 If j < i Then
 xd = xr(i): xr(i) = xr(j): xr(j) = xd: xd = xi(i): xi(i) = xi(j): xi(j) = xd
 End If
 k = k / 2
 Do While j >= k
 j = j - k: k = k / 2
 Loop
 j = j + k
 Next i
の For ~ Next で 例えば、n=4 の場合は i=3, n=8 の場合は i=7 の時
  Do While j >= k
  j = j - k: k = k / 2
  Loop
の箇所で、 j=0, k=0 となり 無限ループに陥ってしまいます。
k = k / 2 のところで k / 2 の k が 1 となり、1/2 の Integer
で 0 となってしまい、Do While 0>= 0 で抜け出せなくなって
しまいます。
気がつきましたので参考まで報告させていただきました。
逆フーリエのセルもあるとうれしいのですが。

 

投稿: kimura | 2007年10月 2日 (火) 23時27分

一瞬ドキッとしましたが、このままで問題ないようです。
実際の使用時には、n,j,kが無限ループを起こす数をとることはありません。

投稿: tsuyumoto | 2007年10月 4日 (木) 16時27分

このプログラムを実行した結果。
横軸は周波数になると思うのですが、縦軸は何になるのですか?
教えてください

投稿: aki | 2007年11月15日 (木) 11時12分

掲載されていたプログラム、とても参考になりました。
少しの改造で自分の使いたいように変えることができ、非常に助かります。

実際の観測データを細かく解析したいときに使わせていただきます。

投稿: Yohei | 2008年3月 2日 (日) 00時39分

このプログラムを用いてローパスフィルタ処理をしようと思っています。
データ数は1024、サンプリング間隔は0.0001mです。
このとき、400cycles/m以上をカットしたいのですが、どうすればよろしいですか?
自分で行った方法です。
波形データから周波数領域には変換しました。
このあと、右側半分を削除しました。
512個目が5000cycles/mということで、400cycles/m以上のデータカットとして41番目以降を削除しました。(5000/400で計算しまし
た)
もとのデータが1024のために、不足部分には0(実数部、虚数部)を入力して、左右対称になるようにコピー貼り付けにて、データを作成
しました。このあとマクロを実行しても、思ったような波形になりませんでした。

投稿: takeshi | 2008年4月24日 (木) 00時59分

実際のデータを見ていないので、あくまでも予測です。

フルスケールが5000Hzで、400Hz以上をカットするのはバランスが悪い気がします(ほとんどがカットされている)。たとえば、フルスケール1000Hzすなわち、サンプリング間隔0.5msくらいで、サンプリングされる方がよいのではないでしょうか。

投稿: tsuyumoto | 2008年4月24日 (木) 23時18分

掲載されたプログラム、ものすごく便利になります。
変換と逆変換を実行しましたが、確かに言われた通りで計算誤差が発生しました。ところで、この計算誤差をどうやって小さくすることができますか。
是非ご教授ください。

投稿: れい 其俊 | 2008年7月21日 (月) 16時40分

時間データのサンプリングが2.56倍なのですが、周波数はどう見れば良いですか。周波数を1.28倍と見ないといけないですか?

投稿: たなか | 2009年8月18日 (火) 18時32分

大変参考になり、さっそく使用させていただきました。
このFFTを利用して、ハイパス、ローパス、バンドパスフィルタの作り方を教えていただけませんか。
EXCEL VBAで実現する方法がよくわかりません。
よろしくお願いします。

投稿: miguel | 2010年8月28日 (土) 22時42分

4096個のデータについてFFTを実行する必要があり、こちらのプログラムを参考に、配列サイズを変更して、g, h, i, j, k, l, m, n, o, p, qをすべてIntegerになるように書き換えたところ、n=4096のとき無限ループに なってしまいました。

投稿: さいと | 2011年9月12日 (月) 19時10分

65536個のデータについてFFTの実行をしたところ、オーバーフローしてしまいました。これを回避するにはどうすればよいのでしょうか。

投稿: takasaki | 2011年12月 5日 (月) 18時10分

3.14159265359 は システム定数がありますので
計算精度として定数が望ましいと思います。

投稿: システム定数 | 2013年5月 7日 (火) 11時53分

周波数20Hzの正弦波を作成し、こちらのプログラムでフーリエ変換したところピーク周波数が20.5Hzとなってしまいました。この誤差の原因はなんでしょうか。

投稿: kawasaki | 2015年5月 6日 (水) 21時04分

今更感もあるかもしれませんが、各変数をintegerにしたところ、さいとさんのおっしゃるとおり、無限ループに陥りました。

Do While j >= k
j = j - k: k = k / 2
Loop

のところで、整数で丸められるため、j = k = 0 となり、無限ループするようです。
ですので、j と k のみは single ないし double で宣言する必要があります。
こうすると最後に、 j = 0 , k = 0.5 となり、ループを抜けることを確認しました。

最後になりましたが、本当に便利なVBAをありがとうございます。
これですぐに周波数特性を確認できる環境を整えられました。

投稿: tanuki | 2016年8月 3日 (水) 21時44分

コメントを書く



(ウェブ上には掲載しません)




トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/187788/14449114

この記事へのトラックバック一覧です: 高速フーリエ変換 Excel VBA用FFTプログラム:

« 能登半島地震 記録と雑感 | トップページ | ホウ素とフッ素の排水基準、温泉旅館への適用は延期へ »