Go言語(golang)で FFT してみたかった
高速フーリエ変換(FFT, Fast Fourier Transform)と逆高速フーリエ変換(IFFT, Inverse Fast Fourier Transform)のシンプルなサンプルです。
とりあえず、バイトのスライス(可変長配列)をフーリエ変換して、それを元に逆フーリエ変換するだけです。
既存のライブラリを使えば車輪の再発明とならずに済みます。
しかし、ツヨツヨなパイセンの言葉に「車輪を作ったことがなければ、車輪造りの難しさ(構造の理解)や、(ライブラリの)ありがたさを知ることはできないし、いざと言う時メンテナンスもできない」というのがあり、そもそも作ってみないと(触らないと)理解できないタイプなので、長年ペンディングにしていた(そっ閉じしていた)チャレンジをチートせずにトライしてみました。
それでも、やっぱり理解したとは言えないことを理解しました。とほほ。
- 参考文献
- Fast Fourier transform @ Wikipedia
- Cooley–Tukey FFT algorithm @ Wikipeida
- Fast Fourier Transformation FFT - Basics @ NTi Audio
- FFT Analysis (Fast Fourier Transform): The Ultimate Guide to Frequency Analysis by Søren Linnet Gjelstrup | Blog @ DWESoft
- fft 高速フーリエ変換 @ Mathworks
main.go
package main
import (
"fmt"
"math/cmplx"
)
func main() {
const sampleSpacing = 0.02 // 50 Hz = 1 / 50 = 0.02 seconds
originalData := []int{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}
complexData, paddedLength := Real2Compex(originalData)
y := FFT(complexData, paddedLength)
z := IFFT(y, paddedLength)
f := FFTFreq(paddedLength, sampleSpacing)
fmt.Println("The next power of 2 of the length of the original data is:", paddedLength)
fmt.Println(" K | DATA | FFT | IFFT | Freq (Hz) | AMPLITUDE")
for k := 0; k < paddedLength; k++ {
amplitude := cmplx.Abs(y[k])
fmt.Printf("%2d | %s |%8.3f %8.3f |%8.3f %8.3f | %9.3f | %9.3f\n",
k, val2Str(originalData[k]), real(y[k]), imag(y[k]), real(z[k]), imag(z[k]), f[k], amplitude)
}
restoredData := make([]int, len(z))
Complex2Real(z, restoredData)
fmt.Println("Original Data:", originalData)
fmt.Println("Restored Data:", restoredData)
//
// Output:
// The next power of 2 of the length of the original data is: 16
// K | DATA | FFT | IFFT | Freq (Hz) | AMPLITUDE
// 0 | 1 | 8.500 0.000 | 1.000 -0.000 | 0.000 | 8.500
// 1 | 2 | -0.500 2.514 | 2.000 0.000 | 3.125 | 2.563
// 2 | 3 | -0.500 1.207 | 3.000 -0.000 | 6.250 | 1.307
// 3 | 4 | -0.500 0.748 | 4.000 -0.000 | 9.375 | 0.900
// 4 | 5 | -0.500 0.500 | 5.000 -0.000 | 12.500 | 0.707
// 5 | 6 | -0.500 0.334 | 6.000 0.000 | 15.625 | 0.601
// 6 | 7 | -0.500 0.207 | 7.000 0.000 | 18.750 | 0.541
// 7 | 8 | -0.500 0.099 | 8.000 -0.000 | 21.875 | 0.510
// 8 | 9 | -0.500 0.000 | 9.000 -0.000 | -25.000 | 0.500
// 9 | 10 | -0.500 -0.099 | 10.000 0.000 | -21.875 | 0.510
// 10 | 11 | -0.500 -0.207 | 11.000 -0.000 | -18.750 | 0.541
// 11 | 12 | -0.500 -0.334 | 12.000 0.000 | -15.625 | 0.601
// 12 | 13 | -0.500 -0.500 | 13.000 0.000 | -12.500 | 0.707
// 13 | 14 | -0.500 -0.748 | 14.000 -0.000 | -9.375 | 0.900
// 14 | 15 | -0.500 -1.207 | 15.000 0.000 | -6.250 | 1.307
// 15 | 16 | -0.500 -2.514 | 16.000 -0.000 | -3.125 | 2.563
// Original Data: [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]
// Restored Data: [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]
}
fft.go
package main
import (
"fmt"
"math"
"math/bits"
"math/cmplx"
"reflect"
)
// FFT performs the Fast Fourier Transform on the input data.
func FFT(x []complex128, n int) []complex128 {
y := fft(x, n)
complex_n := complex(float64(n), 0.0)
for i := 0; i < n; i++ {
y[i] = y[i] / complex_n
}
return y
}
// IFFT performs the Inverse Fast Fourier Transform on the input data.
func IFFT(x []complex128, n int) []complex128 {
y := make([]complex128, n)
for i := 0; i < n; i++ {
y[i] = cmplx.Conj(x[i])
}
y = fft(y, n)
for i := 0; i < n; i++ {
y[i] = cmplx.Conj(y[i])
}
return y
}
// FFTFreq computes the Discrete Fourier Transform sample frequencies and returns
// a slice of length n containing the sample frequencies.
//
// Arguments:
// - n: The window length. This is the number of samples in the signal.
// - d: The sample spacing (inverse of the sampling rate). This is the time
// interval between samples.
func FFTFreq(n int, dt float64) []float64 {
ndt := float64(n) * dt
f := make([]float64, n)
for i := 0; i < n/2; i++ {
f[i] = float64(i) / ndt
f[n/2+i] = -float64(n/2-i) / ndt
}
return f
}
// RFFTFreq computes the frequency components for the real FFT.
func RFFTFreq(n int, dt float64) []float64 {
ndt := float64(n) * dt
f := make([]float64, n/2+1)
for i := 0; i < n/2+1; i++ {
f[i] = float64(i) / ndt
}
return f
}
// NextPowerOf2 returns the smallest power of 2 greater than or equal to n.
func NextPowerOf2(n int) int {
if n <= 1 {
return 1
}
return 1 << (bits.Len(uint(n - 1)))
}
// Number is a constraint that permits any numeric type.
type Number interface {
~int | ~int8 | ~int16 | ~int32 | ~int64 |
~uint | ~uint8 | ~uint16 | ~uint32 | ~uint64 |
~float32 | ~float64
}
// Real2Compex converts real data to complex data with zero imaginary parts,
// padding the result to the next power of 2.
func Real2Compex[T Number](data []T) ([]complex128, int) {
n := len(data)
nn := NextPowerOf2(n)
cmplxData := make([]complex128, nn)
for i := range n {
cmplxData[i] = complex(float64(data[i]), 0.0)
}
return cmplxData, nn
}
// Complex2Real converts complex data to real data by extracting the real parts.
func Complex2Real[T Number](from []complex128, to []T) {
for i, v := range from {
realValue := real(v)
switch any(to).(type) {
case []int:
to[i] = T(int(math.Round(realValue)))
case []int8:
to[i] = T(int8(math.Round(realValue)))
case []int16:
to[i] = T(int16(math.Round(realValue)))
case []int32:
to[i] = T(int32(math.Round(realValue)))
case []int64:
to[i] = T(int64(math.Round(realValue)))
case []uint:
to[i] = T(uint(math.Round(realValue)))
case []uint8:
to[i] = T(uint8(math.Round(realValue)))
case []uint16:
to[i] = T(uint16(math.Round(realValue)))
case []uint32:
to[i] = T(uint32(math.Round(realValue)))
case []uint64:
to[i] = T(uint64(math.Round(realValue)))
case []float32:
to[i] = T(float32(realValue))
case []float64:
to[i] = T(realValue)
default:
panic(fmt.Sprintf("unsupported type: %v", reflect.TypeOf(to).Elem()))
}
}
}
func val2Str[T any](v T) string {
dataValue := reflect.ValueOf(v)
var dataStr string
switch dataValue.Kind() {
case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64:
dataStr = fmt.Sprintf("%6d", dataValue.Int())
case reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64:
dataStr = fmt.Sprintf("%6d", dataValue.Uint())
case reflect.Float32, reflect.Float64:
dataStr = fmt.Sprintf("%6.1f", dataValue.Float())
default:
dataStr = fmt.Sprintf("%6v", v)
}
return dataStr
}
// fft is an internal function that performs the FFT.
func fft(a []complex128, n int) []complex128 {
x := make([]complex128, n)
copy(x, a)
j := 0
for i := 0; i < n; i++ {
if i < j {
x[i], x[j] = x[j], x[i]
}
m := n / 2
for m > 1 && j >= m {
j -= m
m /= 2
}
j += m
}
for kmax := 1; kmax < n; kmax *= 2 {
for k := 0; k < kmax; k++ {
theta := complex(0.0, -math.Pi*float64(k)/float64(kmax))
for i := k; i < n; i += 2 * kmax {
j := i + kmax
temp := x[j] * cmplx.Exp(theta)
x[j] = x[i] - temp
x[i] = x[i] + temp
}
}
}
return x
}
- オンラインで動作を見る @ GoPlayground