高速フーリエ変換(FFT)ツールの使い方

高速フーリエ変換(FFT)ツールのマニュアルです。

利用環境

Internet Explorerには対応していません。Google Chrome、Microsoft Edgeなどのブラウザをご使用ください。スマートフォンでの利用は推奨しません。パソコンでご利用ください。
入力された条件や計算結果などは、外部のサーバーには送信されません。計算はすべて、ご使用のパソコン上で行われます。

使用方法

1.データの入力

まず、「入力データ」欄でデータを入力します。

最初に、[データファイル読込]ボタンでファイルを指定し読み込みます。データファイルは、テキスト形式の時系列データで、時間と値の2列のデータ、

0.000976563	3.426735648
0.001953125	1.677066196
0.002929688	-0.785192212

もしくは、値のみの1列のデータ、

3.28471565
3.520201378
3.665739561

の形式です。2列の場合は、空白またはカンマで区切られているものとします。時系列データは一定の時間間隔で出力されている必要があります。

ファイルが読み込まれると次の各ボックスに値が表示されます。

データ点数: 読み込んだデータの数です。最大で131,072(=217)個のデータを読むことができます。

サンプリング数: FFTのサンプリング数です。読み込んだデータの数が入力されますが、変更することができます。サンプリング数がデータ点数より小さい場合、データの先頭からサンプリング数までで計算されます。サンプリング数がデータ点数より大きい場合、データはゼロ埋めされて計算されます。サンプリング数は2の累乗である必要はありません。サンプリング数の最大は、131,072(=217)です。

サンプリング周期: 時系列の時間間隔です。2列の時系列データを読み込んだ場合は、時間から自動的にサンプリング周期が計算されます。1列のデータの場合は、1が入力されます。いずれも変更可能です。また、時間の単位を選択します。

サンプリング周波数: サンプリング周期からサンプリング周波数が自動で計算され表示されます。

窓関数: 窓関数を選択します。「なし」は矩形窓に相当します。「ハニング」「ハミング」「バートレット」「ブラックマン」の各窓関数を選択できます。

フィルタ周波数: 「低」の周波数より低い周波数と、「高」より高い周波数をFFT結果から除去し逆FFTします。このフィルタ処理により、もとの入力波形を滑らかにすることができます。どちらか一方もしくは両方指定できます。周波数を指定しない場合はフィルタ処理されません。例えばローパスフィルタの場合は、「高」のみ指定します。

入力できたら、[FFT計算実行]ボタンでFFTを実行します。

2.入力データグラフの表示

計算が実行されると、「入力データグラフ」欄に入力データのグラフが表示されます。

軸x: x軸の「最小」「最大」「刻み」を指定できます。指定しない場合は自動で表示されます。

軸y: y軸の表示を指定します。項目は「入力データ」「入力データ(窓関数)」「フィルタ処理データ」から選択できます。軸の「最小」「最大」「刻み」を指定できます。指定しない場合は自動で表示されます。

入力グラフ表示: 軸設定や窓関数のチェックなどを変更した場合は、[入力グラフ表示]ボタンを押すとグラフが再描画されます。

画像保存: 表示されているグラフが input_graph.png というファイル名でPNG形式で保存されます。

データ保存: 表示されているグラフのデータが input_data.csv というファイル名でCSV形式で保存されます。

3.FFT結果グラフの表示

計算が実行されると、「FFT結果グラフ」欄にFFT結果のグラフが表示されます。

軸x: x軸の表示を指定します。項目は「周波数」「波数」から選択できます。周波数の場合は単位を選択します。軸の「最小」「最大」「刻み」を指定できます。「対数」にチェックを入れると、対数表示になります。

軸y: y軸の表示を指定します。項目は「振幅」「位相」「パワー」の各スペクトルから選択できます。各項目の詳細は次項で説明しています。軸の「最小」「最大」「刻み」を指定できます。「対数」にチェックを入れると、対数表示になります。

正規化: チェックを入れるとデータが正規化されて表示されます。正規化の詳細は次項で説明します。

窓関数補正: チェックを入れると窓関数の補正が考慮されたデータで表示されます。補正の詳細は次項で説明しています。

デシベル基準値: パワースペクトルのデシベル[dB]表示の際の基準値を指定します。デフォルトでは音圧レベルの基準値 20×10-6[Pa]が入っています。

結果グラフ表示: 軸設定などを変更した場合は、[結果グラフ表示]ボタンを押すとグラフが再描画されます。

画像保存: 表示されているグラフが fft_graph.png というファイル名でPNG形式で保存されます。

データ保存: 表示されているグラフのデータが fft_data.csv というファイル名でCSV形式で保存されます。

基礎式

離散フーリエ変換(DFT)

ある時間間隔 $\Delta$(サンプリング周期)で測定された $N$ 個のデータ $x$ を考えます。

$$ x(0), \ x(\Delta), \ x(2 \Delta), \ \cdots , \ x(n \Delta), \ \cdots , \ x((N-1) \Delta)$$

このとき、離散フーリエ変換(Discrete Fourier Transform: DFT)は、

$$ X(f_k) = \sum_{n=0}^{N-1} x(n \Delta) e^{-jkn2 \pi / N}$$

となります。$f_k$ は周波数、$k$ は波数で、

$$f_k = \frac{k}{N \Delta} \qquad (k=0, \ 1, \ 2, \ \cdots, \ N-1)$$

です。最も高い周波数は、$k=N$ の時で、サンプリング周波数といい $f_s = 1 / \Delta$ となります。離散フーリエ変換では、$f_0 = 0$ から $f_{N-1} = (N-1) f_\Delta$ まで $f_\Delta = 1 / (N \Delta)$ 刻みで、フーリエ係数 $X(f_k)$ が求まります。

また、離散逆フーリエ変換(Inverse Discrete Fourier Transform: IDFT)は、

$$ x(n \Delta) = \frac{1}{N} \sum_{k=0}^{N-1} X(f_k) e^{jkn2 \pi / N}$$

となります。

高速フーリエ変換(FFT)

離散フーリエ変換は測定データ数が大きくなると演算時間が急増するため、この演算を高速に計算するためのアルゴリズムとして高速フーリエ変換(Fast Fourier Transform: FFT)が考案されました。高速フーリエ変換はデータ数が2の累乗に限定されるものが一般的ですが、ここではBluesteinのFFTアルゴリズム(Chrip Z-変換)を使用しているため、データ数が2の累乗である必要はありません。

参考文献: Rabiner, L. R.; Schafer, R. W.; Rader, C. M. The chirp z-transform algorithm and its application. Bell Syst. Tech. J. 1969, vol. 48,  p. 1249–1292. 

結果項目

FFT結果で表示される項目について説明します。

振幅スペクトル
 振幅は次式で定義されます。
$$ |X(f_k)|$$

位相スペクトル
 位相は次式です。
$$\angle X(f_k) = \tan^{-1} (X_{Im} /X_{Re})$$

 $X_{Im}$ は$X(f_k)$ の虚部、$X_{Re}$ は実部を表します。

パワースペクトル
 パワーは次式です。
$$ |X(f_k)| ^2$$

パワースペクトル密度
 パワーを周波数分解能 $f_\Delta = 1 / (N \Delta)$ で割った値です。
$$ PSD = |X(f_k)| ^2 / f_\Delta$$

パワースペクトル[dB] 
 パワースペクトルのデシベル[dB]表示です。
$$10 \log_{10} \left(\frac{|X(f_k)| ^2}{p_0^2} \right)$$

 $p_0$ はデシベルの基準値です。音圧レベルの場合は 20×10-6[Pa]です。

A特性パワースペクトル[dB]
 A特性の音圧レベルを考慮した表示です。

C特性パワースペクトル[dB]
 C特性の音圧レベルを考慮した表示です。

正規化

正規化にチェックを入れると次の2つの正規化が行われたデータを表示します。

正規化1
離散フーリエ変換の定義から、求まる振幅 $ |X(f_k)|$ を入力データの振幅に合わせるため $1/N$ 倍して振幅を計算します。また、パワースペクトル $ |X(f_k)| ^2$ は $1/N^2$ 倍されます。

正規化2
フーリエ変換された $N$ 個のスペクトル(振幅やパワー) は、サンプリング周波数の 1/2 の周波数(ナイキスト周波数)を堺に左右対称となります。グラフ表示では、ナイキスト周波数までを表示しています。スペクトルの値は対になる対称成分を足し合わせたものが、入力データの実データと一致するため、正規化ではスペクトル値を2倍して表示します。ただし、DC成分($f_0$)やナイキスト周波数成分など対になる成分がない場合はそのまま表示されます。

窓関数補正

窓関数を使用すると入力信号が小さくなるため、振幅やパワーが小さくなります。そこで、もとの入力データと比較できるようにするため、振幅やパワーの値を補正して表示します。

波数 $k$ に対する窓関数を $w_k$ とすると振幅の補正値は、

$$ \sum_{k=0}^{N-1} w_k / N $$

パワーの補正値は、

$$ \sum_{k=0}^{N-1} w_k^2 / N $$

とします。この補正値は窓関数がある場合と窓関数が無い場合(矩形窓)との面積比に相当するものです。振幅やパワーの値をこれらの補正値で割ることで、出力値を補正しています。

スポンサーリンク
科学技術計算のご相談は「キャットテックラボ」へ

科学技術計算やCAEに関するご相談、計算用プログラムの開発などお困りのことは「株式会社キャットテックラボ」へお問い合わせください。

お問い合わせはこちら

フォローする