0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

【Golang】FFT と IFFT の実働サンプル

Last updated at Posted at 2025-01-04

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
}

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

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

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

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

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

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?