高速フーリエ変換（FFT)

　　（プログラムを参考にしてください。）

   Excel Macro   FFT   

最初にExcel Macroがあります。

（Excelのマクロのみです。）

Sheet1に

N,1,1,DT,

Time(i),X(i),Y(i)

のデータを貼り付け、実行すると、

Sheet2に計算結果が出力されます。

Sheet2には、周波数分析結果の

Fxr(i),Fxi(i),Power X(i),Fyr(i),Fyi(i),Power Y(i),TF(i),Phase(i)

が出力されます。

これをグラフにしたもの伝達関数の図をつけました。

macroのコピーのやり方。

新しいExcelのファイルを用意する。(***.xlsmで）

表示→マクロの表示→マクロ作成のフォームが

開くので、sub名をFFTとして入れる。

マクロ作成画面になるので、End subを消して、

1.Excel macro   FFT  MACROと書かれた、

以下のマクロをコピーして（Sub FFT()を除く）

貼り付ける。

Sheet1にデータを貼り付け実行すれば、

Sheet 2に計算結果が表示されます。

Sheet1のデータ形式は、マクロを読んで

考えてね。

貼り付けるデータは、impuls_response.txtを

開いてSheet1の[A1]セルに貼り付ける。

カンマで区切られたデータです。

このデータの作成方法は、

FORTRANのimplus_response.forをコンパイル

して実行すると、impuls_response.txtができます。

試してみてね。固有振動数は1Hzです。

このマクロ（1．Excel macro  FFT  MACRO)は、

実験データを貼り付けて、実測データの周波数分析ができます。応用してね。（データ数は、2^nです。、データ長により周波数分析した振幅が替ります。伝達関数は、Y/Xで相殺され、振幅は、応答倍率となります）

貼り付けるデータの1行目は、

N：データ数、1，1，dtです。

この、1，1、は意味はありません。

ここでは、4－DFTに掲載されているインパルス応答のプログラムの出力結果を入力データとして用い、駆動点伝達関数を求めました。伝達関数の図は、それです。

謝辞

日野幹雄先生の本をだいぶ参考にさせて頂きました。良く書かれた本で、奥深いところまで書いてあります。

FFTを参考にさせて頂きました。FORTRANのCOMPLEXが使用されておりましたが、私が式を展開して実数だけで計算できるように変更しました。そしてそれをExcelのマクロに書き替えました。

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

日野幹雄著　スペクトル解析　朝倉書店

1．Excel macro        FFT  MACRO

       Sub FFT()

      Dim xx(16384) As Single

      Dim yy(16384) As Single

      Dim xC(16384) As Single

      Dim xS(16384) As Single

      Dim yC(16384) As Single

      Dim yS(16384) As Single

      Dim fxr(8192) As Single

      Dim fxi(8192) As Single

      Dim fyr(8192) As Single

      Dim fyi(8192) As Single

      Dim px(8192) As Single       ' x   power spectrum

      Dim py(8192) As Single       ' y   poer spectrum

      Dim tfr(8192) As Single      ' real part of transfer function   TF=y/x  x : input   y : output

      Dim tfi(8192) As Single      ' imaginary part of transfer function   TF=y/x

      Dim tfa(8192) As Single      ' amplitude of transfer function

      Dim phi(8192) As Single      ' phase of transfer function

      Worksheets("Sheet1").Select

      [A1].Select

      N = ActiveCell.Offset(0, 0)

      dt = ActiveCell.Offset(0, 3)

      For i = 1 To N

      xx(i) = ActiveCell.Offset(i, 1)

      yy(i) = ActiveCell.Offset(i, 2)

      Next i

      ND = 16384

    '  N2 = N

    '  NC2 = N2 / 2

      NB = N

      NB2 = NB / 2

      Dim A(8192) As Single

      Dim B(8192) As Single

      Dim freq(8192) As Single

      Tend = (N - 1) * dt

      f1 = 1# / Tend

      aka = Tend / N

      bka = -Tend / N

      For i = 1 To N

      xC(i) = xx(i)

      yC(i) = yy(i)

      xS(i) = 0#

      yS(i) = 0#

      Next i

      For i = 1 To N2

      A(i) = 0#

      B(i) = 0#

      Next i

      IND = 1

        '                  8192,4096,8192,2048

      Call subFFT(xC, A, B, ND, N2, NB, NB2, IND)

      For i = 1 To NB2

      fxr(i) = A(i) * aka

      fxi(i) = B(i) * bka

      freq(i) = f1 * i

      Next i

      Call subFFT(yC, A, B, ND, N2, NB, NB2, IND)

      For i = 1 To NB2

      fyr(i) = A(i) * aka

      fyi(i) = B(i) * bka

      Next i

      pai = 3.14159265358979

      For m = 1 To NB2

      px(m) = fxr(m) * fxr(m) + fxi(m) * fxi(m)

      py(m) = fyr(m) * fyr(m) + fyi(m) * fyi(m)

      If (px(m) = 0) Then

      tfr(m) = 1#

      tfi(m) = 1#

      Else

      tfr(m) = (fxr(m) * fyr(m) + fxi(m) * fyi(m)) / px(m)   ' real part of transfer function

      tfi(m) = (fxr(m) * fyi(m) - fyr(m) * fxi(m)) / px(m)   ' imaginary part of transfer function

      End If

      tfa(m) = Sqr(tfr(m) * tfr(m) + tfi(m) * tfi(m))        ' amplitude of transfer function

      atfr = tfr(m)

      atfi = tfi(m)

      phi(m) = arctan(atfr, atfi)                            ' phase of transfer function

      Next m

      Worksheets("Sheet2").Select

      [A1].Select

      ActiveCell.Offset(0, 0) = " freq"

      ActiveCell.Offset(0, 1) = "fxr"

      ActiveCell.Offset(0, 2) = "fxi"

      ActiveCell.Offset(0, 3) = "Power X"

      ActiveCell.Offset(0, 4) = "fyr"

      ActiveCell.Offset(0, 5) = "fyi"

      ActiveCell.Offset(0, 6) = "Power Y"

      ActiveCell.Offset(0, 7) = "TFr"

      ActiveCell.Offset(0, 8) = "TFi"

      ActiveCell.Offset(0, 9) = "TF"

      ActiveCell.Offset(0, 10) = "PHASE"

      For i = 1 To NB2

      ActiveCell.Offset(i, 0) = freq(i)

      ActiveCell.Offset(i, 1) = fxr(i)

      ActiveCell.Offset(i, 2) = fxi(i)

      ActiveCell.Offset(i, 3) = fxr(i) * fxr(i) + fxi(i) * fxi(i)

      ActiveCell.Offset(i, 4) = fyr(i)

      ActiveCell.Offset(i, 5) = fyi(i)

      ActiveCell.Offset(i, 6) = fyr(i) * fyr(i) + fyi(i) * fyi(i)

      ActiveCell.Offset(i, 7) = tfr(i)

      ActiveCell.Offset(i, 8) = tfi(i)

      ActiveCell.Offset(i, 9) = tfa(i)

      ActiveCell.Offset(i, 10) = phi(i) * 180 / pai

      Next i

      End Sub

      Sub subFFT(xC, A, B, ND, N2, NB, NB2, IND)

      'temp , THETA

      Dim xS(16384) As Single

      j = 1

      For i = 1 To NB

      If (i >= j) Then GoTo B110:

      tempc = xC(j)

      temps = xS(i)

      xC(j) = xC(i)

      xS(j) = xS(i)

      xC(i) = tempc

      xS(i) = temps

B110:   m = NB / 2

B120:   If (j <= m) Then GoTo B130:

      j = j - m

      m = m / 2

      If (m >= 2) Then GoTo B120:

B130:   j = j + m

      Next i

      KMAX = 1

B150:   If (KMAX >= NB) Then GoTo B200:

      ISTEP = KMAX * 2

      For K = 1 To KMAX

      thetac = 0#

      thetas = 3.141592654 * (IND * (K - 1)) / (KMAX)

      For i = K To NB Step ISTEP

      j = i + KMAX

      acexpr = Cos(thetas)

      acexpi = Sin(thetas)

      tempc = xC(j) * acexpr - xS(j) * acexpi

      temps = xS(j) * acexpr + xC(j) * acexpi

      xC(j) = xC(i) - tempc

      xS(j) = xS(i) - temps

      xC(i) = xC(i) + tempc

      xS(i) = xS(i) + temps

     ' temp = x(j) * CEXP(THETA)

     ' x(j) = x(i) - temp

     ' x(i) = x(i) + temp

      Next i

      Next K

      KMAX = ISTEP

      GoTo B150:

B200: For i = 1 To NB2

    '  A(i) = REAL(x(i))

    '  B(i) = IMAG(x(i))

      A(i) = xC(i)

      B(i) = xS(i)

      Next i

      End Sub

      Function arctan(x, y)

      pai = 3.14159265358979

      If x = 0 And y = 0 Then

      arctan = 0

      ElseIf x > 0 And y = 0 Then

      arctan = ((pai / 2) * 0)

      ElseIf x = 0 And y > 0 Then

      arctan = ((pai / 2) * 1)

      ElseIf x < 0 And y = 0 Then

      arctan = ((pai / 2) * 2)

      ElseIf x = 0 And y < 0 Then

      arctan = (-(pai / 2))

      ElseIf x > 0 And y > 0 Then

      arctan = Atn(y / x)

      ElseIf x < 0 And y > 0 Then

      arctan = Atn(y / x) + pai

      ElseIf x < 0 And y < 0 Then

      arctan = -pai + Atn(Abs(y) / Abs(x))

      ElseIf x > 0 And y < 0 Then

      arctan = Atn(y / x)

      End If

      End Function

マクロは、ここまで。