アルゴリズム (FFT)


目次

離散フーリエ変換(DFT)は、時間領域の信号を周波数領域に変換する処理です。(\[x_i\]) を長さNのデータとすると、そのDFTは次式で与えられるデータ (\[F_n\]) となります。

\[F_n=\sum_{i=0}^{N-1}x_ie^{-\frac{2\pi j}{N}ni}\]

Originは、FFTWライブラリを使ってフーリエ変換を実行します。Originは、変換されたデータを使って、振幅、マグニチュード、パワー密度を計算します。

FFTW

FFTWでは、変換されたデータの計算は、“codelets”と呼ばれるC言語のコードブロックで構成されるエグゼキュータにより実行されます。各codeletは、変換の一部に指定されています。これらのcodeletsを使って、エグゼキュータはCooley-TurkeyのFFTアルゴリズムを実行し、これの考えは入力信号のサイズをファクタリングするものです。再帰的ファクタリングによって、信号は短く分けられます。短く分けられた変換の結果は乗算され、最終的に元の信号の変換が計算されます。FFTWについての詳細は、http://fftw.org/をご覧下さい。

パワー密度

定義式により、パワー密度とスペクトルは次の式で計算できます。

\[P_{xx}(e^{j\omega })=\sum_{m=-\infty }^\infty r_{xx}(m)e^{-j\omega m}\]

ここで\(r_{xx}(m)\,\!\)は入力信号の自動相関関数です。

しかし、入力信号は有限であり、ある方法がパワースペクトルを推定するのに使えるだけなので、定義式でパワースペクトルを計算することはできません。Originで使われる方法は、Periodogramで、これはフーリエ変換した振幅からパワーを推定するものです。振幅の2乗がパワースペクトルの振幅に比例することは一般的ですが、各ドメインでパワースペクトルの規格化について議論するさまざまな規則が存在します。Originでは、二乗振幅(MSA)、二乗和振幅(SSA)、時間積分2乗振幅 (TISA)の3つを使っています。これらは次のように表されます。

\[Power Density(two-sided)=\begin{cases}\frac{{Re}^2+{Im}^2}{n^2},for MSA\\\frac{{Re}^2+{Im}^2}n,for SSA\\\frac{\Delta t({Re}^2+{Im}^2)}n,for TISA\end{cases}\]

ここで\(Re\,\!\)および\(Im\,\!\)は変換データの実数部と虚数部で、\(n\,\!\)は入力データの長さ、 \(\Delta t\,\!\)はサンプリング間隔です。

パワースペクトルは、スペクトルの種類 (st)で、両側(2)または片側(1)が選択されているかどうかにより、片側または両側にすることができます。片側のパワー密度を計算するには、最初に両側のパワーを計算する必要があります。そして、結果は次の式を使って片側のパワーに変換されます。

\[P_s(i)=P_d(i),i=0\,\!\]

\[P_s(i)=2P_d(i),i=1,2,\cdots \frac n2-1\]

ここで\(P_s(i)\,\!\)は片側のパワースペクトルで、\(P_d(i)\,\!\)は両側のパワースペクトルです。

窓関数が適用される場合、パワーの結果は、次式で定義される補償の係数で乗算されます。

\(N/{\sum_{n=0}^{N-1}w(n)^2}\) ここで\(w(n)\)は、下記で定義した窓関数です。

その他の結果

Originは変換したデータのマグニチュード、位相、振幅を計算できます。\(Re\,\!\)\(Im\,\!\)を変換したデータの実数部と虚数部にし、\(n\,\!\)を入力信号のサイズにします。\(\Delta t\,\!\)を使って、サンプリング間隔を表します。norma 変数を0にセットします。(正規化を使用しない)その他の出力は、以下の式で計算されます。

スペクトルの種類
両側
(i=1-n/2 ~ n/2)
スペクトルの種類
片側
(i=0 ~ n/2)

位相

\[\arctan (\frac{Im}{Re})\,\]

マグニチュード

\[\sqrt{Re^2+Im^2}\,\]

振幅

\[\sqrt{Re^2+Im^2}/n\,\]

\(\sqrt{Re^2+Im^2}/n, i=0\mbox{ or }i=n/2\,\)
\(2*\sqrt{Re^2+Im^2}/n, \mbox{ otherwise }\,\)

dB

\[20log(Amplitude)\,\]

dB単位で正規化された振幅

\[dB-max(dB)\,\]

RMS振幅

\[\frac{\sqrt{2}}2Amplitude\,\]

正規化

上記の計算は、実際には、norma 変数がfalseにセットされているという前提に基づいています。この変数がtrueにセットされていると、複素数、実数、虚数、マグニチュード、二乗マグニチュードが正規化されます。位相、パワー、振幅、正規化した振幅、dB、二乗振幅は、norma 変数の影響を受けません。

スペクトルの種類(st)で、両側 (2)が選択されていて、正規化 (norma)がTrueにセットされている場合、複素数、実数、虚数、マグニチュード、二乗マグニチュードの結果は、 \(n\,\)で除算されます。ここで、 \(n\,\)は入力信号の大きさです。

スペクトルの種類(st)で、片側 (1)が選択されていて、正規化 (norma)がTrueにセットされている場合、複素数、実数、虚数、マグニチュード、二乗マグニチュードの結果は、次のように正規化されます。\(res_s'\,\)を正規化した結果にします。

\[res_s'(i) = \begin{cases} res_s(i)/n, & \mbox{if } i=0 \\ 2*res_s(i)/n, & \mbox{otherwise} \end{cases}\]

サンプリング間隔の自動計算

自動的に計算されるサンプリング間隔は、時間データの増加の平均で、これは通常入力信号と結びついているXデータが使われます。結びついているX列が無ければ、行番号が使われます。Originが増加の平均を取得するのに失敗した場合、サンプリング間隔は1にセットされます。

周波数

周波数の列は、サンプリング間隔Fft1 help English files image014.gifおよび、入力データポイント Nの数から取得されます。n番目の周波数データは次式で与えられます。

\[f_n=\frac n{N\Delta t}\]

N個の入力データポイントがある場合、周波数領域も最大周波数 でN個のポイントを持ちます。最大周波数\(f_{\max }\,\!\)\(\frac 1{\Delta t}(1-\frac 1N)\)です。結果の移動オプションが選択されていない場合、変換は0から\(f_{\max }\,\!\)の間で表示されます。そうでない場合、移動した結果は \(-\frac{f_{\max }}2\)から\(\frac{f_{\max }}2\)の間で表示されます。

ウィンドウ法

ウィンドウ法は漏れを抑えるのに使用します。Originで利用できる異なるウィンドウの種類が次のように定義されます。

以下の式で、\(n\)はデータのインデックス、\(N\)はデータセットの総数です。

矩形ウィンドウ

\[w[n]=1\,\!\]

Welchウィンドウ

\[w[n]=1-\left( \frac{n-\frac 12(N-1)}{\frac 12(N+1)}\right) ^2\]

Triangularウィンドウ

奇数: \(w(n)=\frac 2{N+1}(\frac {N+1}2-|n+1-\frac {N+1}2|)\)
偶数: \(w(n)=\frac 2N(\frac N2-|n+1-\frac {N+1}2|)\)

Bartlettウィンドウ

\[w(n)=\frac 2{N-1}(\frac{N-1}2-|n-\frac{N-1}2|)\]

Hanningウィンドウ

\[w[n]=\frac 12[1-\cos (\frac{2\pi n}{N-1})]\]

Hammingウィンドウ

\[w[n]=0.54-0.46\cos (\frac{2\pi n}{N-1})\]

Blackmanウィンドウ

\[w[n]=0.42-0.5\cos (\frac{2\pi n}{N-1})+0.08\cos (\frac{4\pi n}{N-1})\]

Gaussianウィンドウ

\[w[n]=exp(-0.5(Alpha( \frac{2n}{N-1}-1 ))^2) \,\!\]

Kaiserウィンドウ

\[w[n]=I(beta*\sqrt{1-(\frac{2n}{N-1}-1)^2}) / I(beta) \,\!\]