package datascience

import (
	"math"
	"math/cmplx"
)

// FFT calcula la Fast Fourier Transform usando el algoritmo Cooley-Tukey radix-2.
// Si la longitud de data no es potencia de 2, se rellena con ceros (zero-padding).
func FFT(data []float64) []complex128 {
	n := len(data)
	if n == 0 {
		return []complex128{}
	}

	// Calcular la siguiente potencia de 2.
	size := nextPow2(n)

	// Convertir a complex128 con zero-padding.
	x := make([]complex128, size)
	for i := 0; i < n; i++ {
		x[i] = complex(data[i], 0)
	}

	fftRecursive(x)
	return x
}

// nextPow2 retorna la menor potencia de 2 >= n.
func nextPow2(n int) int {
	p := 1
	for p < n {
		p <<= 1
	}
	return p
}

// fftRecursive aplica Cooley-Tukey radix-2 DIT in-place.
func fftRecursive(x []complex128) {
	n := len(x)
	if n <= 1 {
		return
	}

	// Separar pares e impares.
	even := make([]complex128, n/2)
	odd := make([]complex128, n/2)
	for i := 0; i < n/2; i++ {
		even[i] = x[2*i]
		odd[i] = x[2*i+1]
	}

	fftRecursive(even)
	fftRecursive(odd)

	for k := 0; k < n/2; k++ {
		t := cmplx.Rect(1, -2*math.Pi*float64(k)/float64(n)) * odd[k]
		x[k] = even[k] + t
		x[k+n/2] = even[k] - t
	}
}