Go言語（golang）で FFT してみたかった

高速フーリエ変換（FFT, Fast Fourier Transform）と逆高速フーリエ変換（IFFT, Inverse Fast Fourier Transform）のシンプルなサンプルです。

とりあえず、バイトのスライス（可変長配列）をフーリエ変換して、それを元に逆フーリエ変換するだけです。

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
}

// CeilPowerOf2 returns the smallest power of 2 greater than or equal to n.
func CeilPowerOf2(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 := CeilPowerOf2(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.
// Based on Cooley–Tukey algorithm.
func fft(a []complex128, n int) []complex128 {
    // copy 'a' to avoid disruptive changes
	x := make([]complex128, n)
	copy(x, a)

    // bit-reversal
	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
	}

    // butterfly operation. Basic FFT calculations (sum-difference and rotation
    // of complex numbers).
	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) // create twiddle factor
                
				x[j] = x[i] - temp
				x[i] = x[i] + temp
			}
		}
	}

	return x
}

オンラインで動作を見る @ GoPlayground

一般的には「既存のライブラリを使う」という選択肢なのだと思います。

しかし、今は亡きツヨツヨだったパイセンの言葉に「車輪の再発明とバカにするけど、車輪を作ったことがなければ、車輪造りの難しさや、ありがたさを知ることはできないし、いざと言う時自分でメンテナンスもできない」というのがあります。

つまり「構造の理解」に対する解像度や「ライブラリのありがたみ」を知るには、一度自分で作ってみないとわからないということです。

そもそも作ってみないと（触らないと）理解できないタイプなので、長年ペンディングにしていた（そっ閉じしていた）チャレンジにトライしてみました。

結論。やっぱり理解したとは言えないことを理解しました。とほほ。

参考文献
- 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