Implementation of the Cooley-Tukey Radix2 DIT fast Fourier Transform

dft.go

package dft

import (
	"math"
	"math/cmplx"
)

func Radix2DIT(data []complex128) []complex128 {
	l := len(data)
	l-- // power of two >= l
	l |= l >> 1
	l |= l >> 2
	l |= l >> 4
	l |= l >> 8
	l |= l >> 16
	l++

	out := make([]complex128, l) // Allocate output data

	// increase size of input data to make it a pow of 2
	in := append(data, make([]complex128, len(data)-l)...)

	// Perform the DFT
	ditdft(in, out, 0, 1, l, 0)

	return out
}

func ditdft(in, out []complex128, index, jump, n, outdex int) {
	if n == 1 {
		// Trivial case
		out[outdex] = in[index]
		return
	}

	// Split data in half. Calculate DFT(evens) and DFT(odds)
	ditdft(in, out, index, jump<<1, n>>1, outdex)
	ditdft(in, out, index+jump, jump<<1, n>>1, outdex+n>>1)

	// Rejoin data
	// DFT(in)[k]     = DFT(evens)[k] + e^(-2i*pi*k/n) * DFT(odds)[k]
	// DFT(in)[k+n/2] = DFT(evens)[k] - e^(-2i*pi*k/n) * DFT(odds)[k]
	for k := 0; k < n/2; k++ {
		// DFT(evens)[k]
		even := out[outdex+k]
		// e^(-2i*pi*k/n) * DFT(odds)[k]
		odd := out[outdex+k+n>>1] * cmplx.Rect(1, -2*math.Pi*float64(k)/float64(n))

		out[outdex+k] = even + odd
		out[outdex+k+n>>1] = even - odd
	}
}