Excelにはフーリエ変換をする拡張機能がついていますが、とても使いにくいです。シートにデータがあって、取りあえずFFTをかけてみたい、どんな周期を持っているのか知りたい、という時にパッとは使えません。
そこで、取りあえずFFTをかけてグラフを作る、というマクロを作りました。周波数解析（周期解析）をしたい時に便利です。
機能やマクロの内容は次の通りです。
１．Excelのシートにある縦のデータ（列データ）をFFTに必要な個数で取得する。
２．そのデータをFFTにかけます。
３．窓関数（Hann窓、Hamming窓）の処理ができます。
４．FFT後の周期波形、位相波形をグラフにします。
５．周期波形の強度が強いものを１００％ととして指定のレベル以上（カットオフレベル）の周期（周波数）だけを用いて逆FFTをかけて時間波形に戻します。
６．元の波形とレベル選別した波形をグラフにします。
７．この二つの波形の相関係数を計算します。
８．グラフはシートにデータを書かずに作成します。
９．ポイント数は２の乗数に制限されます。
2,4,8,16,32,64,128,256,512,1024,2048,4096,.....です。
使い方は、
１．FFTをかけたい列データのどこかのセル１つを選択します。
２．マクロ「FFTグラフ作成」を起動します。
３．「データの開始行」、「データ数」、「カットオフレベル」、「窓関数」を聞いてきますので順に指定してください。
４．三個のグラフが作成されます。
取得するデータの部分を変更していただければ都合のいいように使いまわせると思います。VBAはいつものことながら、あまりきれいなコードではありません。すいません。
結果の数値データがほしい時は、グラフにしているデータをシートに出すようにしていただければよいと思います。
グラフはExcelのチャートですので、ポイント数が多いと重くなります。ポイント数が多い時（５１２ポイントぐらい以上）に、作成後のグラフを修正しようとするとエラーになることもあります。場合によっては初めからエラーになったりするかもしれません。
Sub FFTグラフ作成()
Dim AC As Range, CLS As Range, CLE As Range
Dim iyS As Long, iyE As Long, ix As Long
Dim inFFT As Long, iwind As Long
Dim fftdataR() As Double, fftdataI() As Double, fftdataAbs() As Double
Dim fftPeriod() As Double, fftFreq() As Double, theta() As Double
Dim fftdataOrg() As Double, DataNum() As Long
Dim PeakLevel As Double, CorrelationCoef As Double, CutOffPercent As Double
Dim strtemp As String, i As Long
Dim fftdataABSperinFFT() As Double, CutOffLevel() As Double
Dim objChart As Object, chartTtl As String, ItemName As String
Const pai = 3.14159265
Application.ScreenUpdating = False
'データの範囲と項目名
Set AC = ActiveCell
Set CLS = AC.End(xlUp) 'データの一番上
Set CLE = AC.End(xlDown) 'データの一番下
If IsNumeric(CLS.Value) = False Then 'データの一番上が数値でない時
ItemName = CLS.Value 'そこを項目名とみなす
iyS = CLS.Offset(1, 0).Row 'その１行下からデータとみなす
Else 'データの一番上が数値の時
ItemName = CLS.End(xlUp).Value 'さらに離れた上を項目名とみなす
iyS = CLS.Row 'その行からがデータとみなす
End If
ix = AC.Column 'データの列
iyE = CLE.Row 'データの一番下の行
Set CLS = Nothing
Set CLE = Nothing
'開始行
strtemp = InputBox("データの開始行を入力してください。" & vbCrLf & "対象の列、開始行、終了行 = " & ix & ", " & iyS & ", " & iyE, "FFTグラフ作成 " & ItemName, iyS)
If strtemp = "" Then Exit Sub
iyS = Val(strtemp)
'ポイント数の決定。２のｎ乗個
inFFT = 2 ^ Int(Log(iyE - iyS + 1) / Log(2))
iyE = iyS + inFFT - 1
strtemp = InputBox("データ数を入力してください。" & vbCrLf & "対象の列、開始行、終了行 = " & ix & ", " & iyS & ", " & iyE, "FFTグラフ作成 " & ItemName, inFFT)
If strtemp = "" Then Exit Sub
inFFT = Val(strtemp)
'カットオフレベルの決定
CutOffPercent = 50
strtemp = InputBox("元データにフィットするときのカットオフレベルを入力してください。", "FFTグラフ作成", CutOffPercent)
CutOffPercent = Val(strtemp)
'データのウィンドウ関数の選択
strtemp = InputBox("元データに窓関数をかけますか？ 0:そのまま（矩形窓） 1:Hann窓 2:Hamming窓。", "FFTグラフ作成", 0)
iwind = Val(strtemp)
If iwind < 0 Or iwind > 2 Then Exit Sub
'配列の設定
ReDim fftdataR(inFFT), fftdataI(inFFT), fftdataAbs(inFFT), fftdataOrg(inFFT), fftFreq(inFFT), DataNum(inFFT)
ReDim fftPeriod(1 To inFFT / 2), theta(1 To inFFT / 2)
ReDim fftdataABSperinFFT(1 To inFFT / 2), CutOffLevel(1 To inFFT / 2)
For i = 0 To inFFT - 1
DataNum(i) = i
fftdataR(i) = Cells(iyS + i, ix) '実時間データ
fftdataI(i) = 0#
Next i
fftdataOrg = fftdataR '元データを残す
'FFT。時間領域データを周波数領域に変換
FFT inFFT, fftdataR, fftdataI, iwind
'周波数領域データの出力。周期を横軸。
PeakLevel = 0
For i = 1 To inFFT / 2 'i=0の直流データは使わない
fftdataAbs(i) = Sqr(fftdataR(i) * fftdataR(i) + fftdataI(i) * fftdataI(i))
fftFreq(i) = i / (inFFT / 2)
fftPeriod(i) = inFFT / i '(timeunit * nnn) / i
If PeakLevel < fftdataAbs(i) Then PeakLevel = fftdataAbs(i)
'-1*atan()/(2pai)*T 過去が+、将来が-なので * -1
theta(i) = WorksheetFunction.Atan2(fftdataR(i), fftdataI(i)) * 180# / pai
Next i
'グラフ用データ
For i = 1 To inFFT / 2
fftdataABSperinFFT(i) = fftdataAbs(i) / inFFT
CutOffLevel(i) = PeakLevel * CutOffPercent / 100 / inFFT 'ピークレベルを使う
Next i
'周期グラフの作成
chartTtl = "周期 " & ItemName & " CutLevel[%]=" & CutOffPercent
Set objChart = plotChart(xlXYScatterLinesNoMarkers, fftPeriod, fftdataABSperinFFT, , , chartTtl) ', , ActiveCell)
Set objChart = plotAddSeries(objChart, fftPeriod, CutOffLevel)
objChart.Chart.Axes(xlCategory).ScaleType = xlLogarithmic
objChart.Chart.Axes(xlCategory).LogBase = 4
objChart.Chart.FullSeriesCollection(1).ChartType = xlXYScatterLines
Call moveChartPosVisibleRange(objChart, 3, 3)
'位相グラフの作成
chartTtl = "位相 " & ItemName
Set objChart = plotChart(xlXYScatterLinesNoMarkers, fftPeriod, theta, , , chartTtl) ', , ActiveCell.Offset(1, 1))
objChart.Chart.Axes(xlCategory).ScaleType = xlLogarithmic
objChart.Chart.Axes(xlCategory).LogBase = 4
Call moveChartPosVisibleRange(objChart, 2, 2)
'周波数領域データの出力。データのカットレベル。
'周波数領域生データの出力
For i = 0 To inFFT - 1
fftdataAbs(i) = Sqr(fftdataR(i) * fftdataR(i) + fftdataI(i) * fftdataI(i))
fftFreq(i) = i / (inFFT / 2)
Next i
'カットレベル以下を０にする
For i = 0 To inFFT - 1
If fftdataAbs(i) <= PeakLevel * CutOffPercent / 100 Then
fftdataR(i) = 0
fftdataI(i) = 0
End If
Next i
'RFT。レベルカットした周波数領域のデータを時間領域に変換
RFT inFFT / 2, fftdataR, fftdataI
'レベルカット前後のデータの相関
CorrelationCoef = CorrCoef(fftdataOrg, fftdataR, 0, inFFT - 1)
CorrelationCoef = 100 * Round(CorrelationCoef, 3)
'時間領域グラフの作成(元データとレベルカット後)
ReDim Preserve fftdataR(inFFT - 1), fftdataOrg(inFFT - 1), DataNum(inFFT - 1)
chartTtl = "FFT波形 " & ItemName & " CutOff " & CutOffPercent & "%" & " Corr " & CorrelationCoef & "%"
Set objChart = plotChart(xlLine, DataNum, fftdataR, , "Filterd", chartTtl, True) ', ActiveCell.Offset(2, 2))
Set objChart = plotAddSeries(objChart, DataNum, fftdataOrg, "Raw")
objChart.Chart.Axes(xlCategory).ReversePlotOrder = True
objChart.Chart.FullSeriesCollection(2).Format.Line.Weight = 0.5
Call moveChartPosVisibleRange(objChart, 1, 1)
AC.Select
Set AC = Nothing
Set objChart = Nothing
Application.ScreenUpdating = True
End Sub
'逆フーリエ変換
'元の周波数領域のデータは、nnn/2を中心に左右対称(real)、左右対称かつ上下反転(imag)（nnn/2を中心点の点対称）なので
'前半部分は正負を逆転させる。縦軸は数値をデータ数nnnで割っておく。
Public Function RFT(nnn As Long, xRorg() As Double, xIorg() As Double)
Dim xr() As Double, xi() As Double
Dim i As Long
ReDim xr(2 * nnn), xi(2 * nnn)
'nnn/2を中心点の点対称データを作る。nnn=4なら -x(0)=単独データ, -x(1)=x(7), -x(2)=x(6), -x(3)=x(5), x(4)=中心データ
For i = 1 To nnn
xr(i - 1) = xRorg(i - 1) / (2 * nnn)
xi(i - 1) = -xIorg(i - 1) / (2 * nnn)
xr(2 * nnn - i) = xRorg(i) / (2 * nnn)
xi(2 * nnn - i) = xIorg(i) / (2 * nnn)
Next i
'FFT実行
FFT 2 * nnn, xr, xi
'結果を元の配列に戻す
For i = 0 To 2 * nnn - 1
xRorg(i) = xr(i)
xIorg(i) = xi(i)
Next i
End Function
'参考
'http://tsuyu.cocolog-nifty.com/blog/2007/03/publi.html
'Numerical Recipes in C
'フーリエ変換。結果のf=i/(n*deltaT), 振幅=abs(x)/(n/2)
'xR実数部、xI虚数部。結果で元データを上書き。0 - nnn-1。
Public Function FFT(nnn As Long, xr() As Double, xi() As Double, Optional iwind As Long = 0)
Dim g As Long, h As Long, i As Long
Dim j As Double, k As Double
Dim l As Long, m As Long, p As Long, q As Long ', n As Long, o As Long
Dim s() As Double, c() As Double
Dim xd As Double, a As Double, b As Double
'窓関数
If iwind <> 0 Then '1:Hann, 2:Hamming
Call Windowing(nnn, xr, xi, iwind)
End If
's,cはポイント数の半分
ReDim s(nnn / 2), c(nnn / 2)
i = 0
j = 0# '最後のdo loop を抜けるためにdoubleとする
k = 0# '最後のdo loop を抜けるためにdoubleとする
l = 0
p = 0
h = 0
g = 0
q = 0
'nnn^m
m = Log(nnn) / Log(2)
'FFTの計算
a = 0
b = 3.14159265359 * 2 / nnn
For i = 0 To nnn / 2
s(i) = sin(a)
c(i) = Cos(a)
a = a + b
Next i
l = nnn
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 = nnn / 2
For i = 1 To nnn - 1
k = nnn
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
End Function
'xR実数部、xI虚数部。結果で元データを上書き。0 - nnn-1。
'WindowFunction w(i)=a-(1-a)*cos2πi/N, i=0 to N-1
'iwind:=1 Hannning Window a=0.5, 0.5-0.5*cos2πx,0<x<1
'iwind:=2 Hammning Window a=0.54, 0.54-0.46*cos2πx,0<x<1
Function Windowing(nnn As Long, xr() As Double, xi() As Double, iwind As Long)
Dim i As Long, alpha As Double
Const PI = 3.14159265358979
If iwind = 1 Then
alpha = 0.5
ElseIf iwind = 2 Then
alpha = 0.54
Else 'if iwind=0 then
alpha = 1#
End If
For i = 0 To nnn - 1
xr(i) = xr(i) * (alpha - (1 - alpha) * Cos(2 * PI * i / nnn))
xi(i) = xi(i) * (alpha - (1 - alpha) * Cos(2 * PI * i / nnn))
Next i
End Function
'チャートを見える範囲内に移動する。addOffsetは内側へ入れるセル数
Function moveChartPosVisibleRange(chartObj As Object, _
Optional addOffsetRow As Long = 1, Optional addOffsetCol As Long = 1)
Dim VR As Range
Dim cl As String, TL As Range, BR As Range
Dim mvRow As Long, mvCol As Long
Dim AC As Range
Set AC = ActiveCell
chartObj.Activate
With chartObj
'画面範囲と現在のチャートの左上、右下
Set VR = ActiveWindow.VisibleRange
Set TL = .TopLeftCell
Set BR = .BottomRightCell
'現在のチャートの左上、右下が画面範囲外にある時の移動量を決める。両方の場合右下を優先する（後で計算する）。
If TL.Row < VR(1).Row + addOffsetRow Then
mvRow = VR(1).Row - TL.Row + addOffsetRow
End If
If TL.Column < VR(1).Column + addOffsetCol Then
mvCol = VR(1).Column - TL.Column + addOffsetCol
End If
If BR.Row > VR(VR.Count).Row - addOffsetRow Then
mvRow = VR(VR.Count).Row - BR.Row - addOffsetRow
End If
If BR.Column > VR(VR.Count).Column - addOffsetCol Then
mvCol = VR(VR.Count).Column - BR.Column - addOffsetCol
End If
'ありえない位置への移動は修正
If TL.Row + mvRow < 1 Then mvRow = 1 - TL.Row
If BR.Row + mvRow > Rows.Count Then mvRow = Rows.Count - BR.Row
If TL.Column + mvCol < 1 Then mvCol = 1 - TL.Column
If BR.Column + mvCol > Columns.Count Then mvCol = Columns.Count - BR.Column
'チャートの左上の位置を変更
cl = TL.Address(False, False)
cl = Range(cl).Offset(mvRow, mvCol).Address(False, False)
.Top = Range(cl).Top
.Left = Range(cl).Left
End With
AC.Select
Set AC = Nothing
Set VR = Nothing
Set TL = Nothing
Set BR = Nothing
End Function
'配列からグラフを作る。参考にしたサイト。
'https://blog.goo.ne.jp/nish0707/e/8d2d2b7ec4bc52f6b23ad42d8c81af50
'http://foundknownanddone.dots2pattern.com/2016/02/excel-vba-plot-column-chart.html
'配列からグラフを作る。xlLine,xlXYScatter系
Function plotChart(ctype, xdata As Variant, ydata As Variant _
, Optional xttl As String, Optional yttl As String, Optional chartTtl As String _
, Optional bLegend As Boolean = False) ', Optional posChart As Range)
Dim oChart As Object
Dim AC As Range
Set AC = ActiveCell
'チャートオブジェクトを適当な位置で作成する
ActiveCell.Resize(5, 5).Select
Set oChart = ActiveSheet.ChartObjects.Add(10, 10, 300, 200)
oChart.Activate
With oChart
With .Chart 'タイトル
.ChartType = ctype 'グラフ形式を設定
With .SeriesCollection.NewSeries
.Values = ydata '１軸目のy値
If yttl <> "" Then
.Name = yttl
End If
If IsEmpty(xdata) = False Then '１軸目のx値
.XValues = xdata
End If
End With
If bLegend = True Then
.HasLegend = True '凡例表示
.Legend.Position = xlLegendPositionBottom 'グラフの下
Else
.HasLegend = False '凡例表示しない
End If
If chartTtl <> "" Then
.HasTitle = True
With .ChartTitle
.Text = chartTtl
.Format.TextFrame2.TextRange.Font.Size = 11
.Format.TextFrame2.TextRange.Font.Bold = msoFalse
End With
End If
If xttl <> "" Then
.Axes(xlCategory, xlPrimary).HasTitle = True
.Axes(xlCategory, xlPrimary).AxisTitle.Characters.Text = xttl
End If
If yttl <> "" Then
.Axes(xlValue, xlPrimary).HasTitle = True
.Axes(xlValue, xlPrimary).AxisTitle.Characters.Text = yttl
End If
.PlotArea.Format.Line.ForeColor.RGB = RGB(64, 64, 64) 'グラフ領域の枠線
End With
End With
AC.Select
Set plotChart = oChart
Set AC = Nothing
Set oChart = Nothing
End Function
'配列からグラフの線を追加する。xlLine,xlXYScatter系
Function plotAddSeries(objChart As Variant, xdata As Variant, ydata As Variant, Optional yttl As String)
objChart.Activate
With objChart
'これを繰り返すと複数のラインが描ける
With .Chart.SeriesCollection.NewSeries
.Values = ydata '追加軸の値
.XValues = xdata
If yttl <> "" Then
.Name = yttl
End If
End With
End With
Set plotAddSeries = objChart
End Function
'x(明示なし)が同じ２つのy1(xxx),y2(yyy)の配列の部分の相関を返す
Function CorrCoef(xxx() As Double, yyy() As Double, iStart As Long, iEnd As Long) As Double
Dim SX As Double, SX2 As Double, SXY As Double, Sy As Double, SY2 As Double
Dim nnn As Long
Dim i As Long
nnn = iEnd - iStart + 1
For i = iStart To iEnd
SX = SX + xxx(i)
SX2 = SX2 + xxx(i) * xxx(i)
Sy = Sy + yyy(i)
SY2 = SY2 + yyy(i) * yyy(i)
SXY = SXY + yyy(i) * xxx(i)
Next i
If ((nnn * SY2 - Sy * Sy) * (nnn * SX2 - SX * SX)) <> 0 Then
CorrCoef = (nnn * SXY - SX * Sy) * (nnn * SXY - SX * Sy) / ((nnn * SY2 - Sy * Sy) * (nnn * SX2 - SX * SX))
CorrCoef = Sqr(CorrCoef)
Else
CorrCoef = 0
End If
'Debug.Print CorrCoef
End Function