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
}

// 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
}
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?