go/xyosc
1
0
Fork 0
mirror of https://github.com/make-42/xyosc synced 2025-01-18 18:57:10 +01:00
Code Issues Actions Packages Projects Releases Wiki Activity

yes

Browse Source
This commit is contained in:
Louis Dalibard 2025-01-16 21:47:06 +01:00
parent 35f0b32420
commit 145abd9383
7 changed files with 957 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
vendor/github.com/mjibson/go-dsp/LICENSE generated vendored Normal file
View File

@ -0,0 +1,13 @@
Copyright (c) 2011 Matt Jibson <matt.jibson@gmail.com>
Permission to use, copy, modify, and distribute this software for any
purpose with or without fee is hereby granted, provided that the above
copyright notice and this permission notice appear in all copies.
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

96
vendor/github.com/mjibson/go-dsp/dsputils/compare.go generated vendored Normal file
View File

@ -0,0 +1,96 @@
/*
* Copyright (c) 2011 Matt Jibson <matt.jibson@gmail.com>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
package dsputils
import (
"math"
)
const (
closeFactor = 1e-8
)
// PrettyClose returns true if the slices a and b are very close, else false.
func PrettyClose(a, b []float64) bool {
if len(a) != len(b) {
return false
}
for i, c := range a {
if !Float64Equal(c, b[i]) {
return false
}
}
return true
}
// PrettyCloseC returns true if the slices a and b are very close, else false.
func PrettyCloseC(a, b []complex128) bool {
if len(a) != len(b) {
return false
}
for i, c := range a {
if !ComplexEqual(c, b[i]) {
return false
}
}
return true
}
// PrettyClose2 returns true if the matrixes a and b are very close, else false.
func PrettyClose2(a, b [][]complex128) bool {
if len(a) != len(b) {
return false
}
for i, c := range a {
if !PrettyCloseC(c, b[i]) {
return false
}
}
return true
}
// PrettyClose2F returns true if the matrixes a and b are very close, else false.
func PrettyClose2F(a, b [][]float64) bool {
if len(a) != len(b) {
return false
}
for i, c := range a {
if !PrettyClose(c, b[i]) {
return false
}
}
return true
}
// ComplexEqual returns true if a and b are very close, else false.
func ComplexEqual(a, b complex128) bool {
r_a := real(a)
r_b := real(b)
i_a := imag(a)
i_b := imag(b)
return Float64Equal(r_a, r_b) && Float64Equal(i_a, i_b)
}
// Float64Equal returns true if a and b are very close, else false.
func Float64Equal(a, b float64) bool {
return math.Abs(a-b) <= closeFactor || math.Abs(1-a/b) <= closeFactor
}

115
vendor/github.com/mjibson/go-dsp/dsputils/dsputils.go generated vendored Normal file
View File

@ -0,0 +1,115 @@
/*
* Copyright (c) 2011 Matt Jibson <matt.jibson@gmail.com>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
// Package dsputils provides functions useful in digital signal processing.
package dsputils
import (
"math"
)
// ToComplex returns the complex equivalent of the real-valued slice.
func ToComplex(x []float64) []complex128 {
y := make([]complex128, len(x))
for n, v := range x {
y[n] = complex(v, 0)
}
return y
}
// IsPowerOf2 returns true if x is a power of 2, else false.
func IsPowerOf2(x int) bool {
return x&(x-1) == 0
}
// NextPowerOf2 returns the next power of 2 >= x.
func NextPowerOf2(x int) int {
if IsPowerOf2(x) {
return x
}
return int(math.Pow(2, math.Ceil(math.Log2(float64(x)))))
}
// ZeroPad returns x with zeros appended to the end to the specified length.
// If len(x) >= length, x is returned, otherwise a new array is returned.
func ZeroPad(x []complex128, length int) []complex128 {
if len(x) >= length {
return x
}
r := make([]complex128, length)
copy(r, x)
return r
}
// ZeroPadF returns x with zeros appended to the end to the specified length.
// If len(x) >= length, x is returned, otherwise a new array is returned.
func ZeroPadF(x []float64, length int) []float64 {
if len(x) >= length {
return x
}
r := make([]float64, length)
copy(r, x)
return r
}
// ZeroPad2 returns ZeroPad of x, with the length as the next power of 2 >= len(x).
func ZeroPad2(x []complex128) []complex128 {
return ZeroPad(x, NextPowerOf2(len(x)))
}
// ToComplex2 returns the complex equivalent of the real-valued matrix.
func ToComplex2(x [][]float64) [][]complex128 {
y := make([][]complex128, len(x))
for n, v := range x {
y[n] = ToComplex(v)
}
return y
}
// Segment returns segs equal-length slices that are segments of x with noverlap% of overlap.
// The returned slices are not copies of x, but slices into it.
// Trailing entries in x that connot be included in the equal-length segments are discarded.
// noverlap is a percentage, thus 0 <= noverlap <= 1, and noverlap = 0.5 is 50% overlap.
func Segment(x []complex128, segs int, noverlap float64) [][]complex128 {
lx := len(x)
// determine step, length, and overlap
var overlap, length, step, tot int
for length = lx; length > 0; length-- {
overlap = int(float64(length) * noverlap)
tot = segs*(length-overlap) + overlap
if tot <= lx {
step = length - overlap
break
}
}
if length == 0 {
panic("too many segments")
}
r := make([][]complex128, segs)
s := 0
for n := range r {
r[n] = x[s : s+length]
s += step
}
return r
}

216
vendor/github.com/mjibson/go-dsp/dsputils/matrix.go generated vendored Normal file
View File

@ -0,0 +1,216 @@
/*
* Copyright (c) 2011 Matt Jibson <matt.jibson@gmail.com>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
package dsputils
// Matrix is a multidimensional matrix of arbitrary size and dimension.
// It cannot be resized after creation. Arrays in any axis can be set or fetched.
type Matrix struct {
list []complex128
dims, offsets []int
}
// MakeMatrix returns a new Matrix populated with x having dimensions dims.
// For example, to create a 3-dimensional Matrix with 2 components, 3 rows, and 4 columns:
// MakeMatrix([]complex128 {
// 1, 2, 3, 4,
// 5, 6, 7, 8,
// 9, 0, 1, 2,
//
// 3, 4, 5, 6,
// 7, 8, 9, 0,
// 4, 3, 2, 1},
// []int {2, 3, 4})
func MakeMatrix(x []complex128, dims []int) *Matrix {
length := 1
offsets := make([]int, len(dims))
for i := len(dims) - 1; i >= 0; i-- {
if dims[i] < 1 {
panic("invalid dimensions")
}
offsets[i] = length
length *= dims[i]
}
if len(x) != length {
panic("incorrect dimensions")
}
dc := make([]int, len(dims))
copy(dc, dims)
return &Matrix{x, dc, offsets}
}
// MakeMatrix2 is a helper function to convert a 2-d array to a matrix.
func MakeMatrix2(x [][]complex128) *Matrix {
dims := []int{len(x), len(x[0])}
r := make([]complex128, dims[0]*dims[1])
for n, v := range x {
if len(v) != dims[1] {
panic("ragged array")
}
copy(r[n*dims[1]:(n+1)*dims[1]], v)
}
return MakeMatrix(r, dims)
}
// Copy returns a new copy of m.
func (m *Matrix) Copy() *Matrix {
r := &Matrix{m.list, m.dims, m.offsets}
r.list = make([]complex128, len(m.list))
copy(r.list, m.list)
return r
}
// MakeEmptyMatrix creates an empty Matrix with given dimensions.
func MakeEmptyMatrix(dims []int) *Matrix {
x := 1
for _, v := range dims {
x *= v
}
return MakeMatrix(make([]complex128, x), dims)
}
// offset returns the index in the one-dimensional array
func (s *Matrix) offset(dims []int) int {
if len(dims) != len(s.dims) {
panic("incorrect dimensions")
}
i := 0
for n, v := range dims {
if v > s.dims[n] {
panic("incorrect dimensions")
}
i += v * s.offsets[n]
}
return i
}
func (m *Matrix) indexes(dims []int) []int {
i := -1
for n, v := range dims {
if v == -1 {
if i >= 0 {
panic("only one dimension index allowed")
}
i = n
} else if v >= m.dims[n] {
panic("dimension out of bounds")
}
}
if i == -1 {
panic("must specify one dimension index")
}
x := 0
for n, v := range dims {
if v >= 0 {
x += m.offsets[n] * v
}
}
r := make([]int, m.dims[i])
for j := range r {
r[j] = x + m.offsets[i]*j
}
return r
}
// Dimensions returns the dimension array of the Matrix.
func (m *Matrix) Dimensions() []int {
r := make([]int, len(m.dims))
copy(r, m.dims)
return r
}
// Dim returns the array of any given index of the Matrix.
// Exactly one value in dims must be -1. This is the array dimension returned.
// For example, using the Matrix documented in MakeMatrix:
// m.Dim([]int {1, 0, -1}) = []complex128 {3, 4, 5, 6}
// m.Dim([]int {0, -1, 2}) = []complex128 {3, 7, 1}
// m.Dim([]int {-1, 1, 3}) = []complex128 {8, 0}
func (s *Matrix) Dim(dims []int) []complex128 {
inds := s.indexes(dims)
r := make([]complex128, len(inds))
for n, v := range inds {
r[n] = s.list[v]
}
return r
}
func (m *Matrix) SetDim(x []complex128, dims []int) {
inds := m.indexes(dims)
if len(x) != len(inds) {
panic("incorrect array length")
}
for n, v := range inds {
m.list[v] = x[n]
}
}
// Value returns the value at the given index.
// m.Value([]int {1, 2, 3, 4}) is equivalent to m[1][2][3][4].
func (s *Matrix) Value(dims []int) complex128 {
return s.list[s.offset(dims)]
}
// SetValue sets the value at the given index.
// m.SetValue(10, []int {1, 2, 3, 4}) is equivalent to m[1][2][3][4] = 10.
func (s *Matrix) SetValue(x complex128, dims []int) {
s.list[s.offset(dims)] = x
}
// To2D returns the 2-D array equivalent of the Matrix.
// Only works on Matrixes of 2 dimensions.
func (m *Matrix) To2D() [][]complex128 {
if len(m.dims) != 2 {
panic("can only convert 2-D Matrixes")
}
r := make([][]complex128, m.dims[0])
for i := 0; i < m.dims[0]; i++ {
r[i] = make([]complex128, m.dims[1])
copy(r[i], m.list[i*m.dims[1]:(i+1)*m.dims[1]])
}
return r
}
// PrettyClose returns true if the Matrixes are very close, else false.
// Comparison done using dsputils.PrettyCloseC().
func (m *Matrix) PrettyClose(n *Matrix) bool {
// todo: use new slice equality comparison
for i, v := range m.dims {
if v != n.dims[i] {
return false
}
}
return PrettyCloseC(m.list, n.list)
}

94
vendor/github.com/mjibson/go-dsp/fft/bluestein.go generated vendored Normal file
View File

@ -0,0 +1,94 @@
/*
* Copyright (c) 2011 Matt Jibson <matt.jibson@gmail.com>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
package fft
import (
"math"
"sync"
"github.com/mjibson/go-dsp/dsputils"
)
var (
bluesteinLock sync.RWMutex
bluesteinFactors = map[int][]complex128{}
bluesteinInvFactors = map[int][]complex128{}
)
func getBluesteinFactors(input_len int) ([]complex128, []complex128) {
bluesteinLock.RLock()
if hasBluesteinFactors(input_len) {
defer bluesteinLock.RUnlock()
return bluesteinFactors[input_len], bluesteinInvFactors[input_len]
}
bluesteinLock.RUnlock()
bluesteinLock.Lock()
defer bluesteinLock.Unlock()
if !hasBluesteinFactors(input_len) {
bluesteinFactors[input_len] = make([]complex128, input_len)
bluesteinInvFactors[input_len] = make([]complex128, input_len)
var sin, cos float64
for i := 0; i < input_len; i++ {
if i == 0 {
sin, cos = 0, 1
} else {
sin, cos = math.Sincos(math.Pi / float64(input_len) * float64(i*i))
}
bluesteinFactors[input_len][i] = complex(cos, sin)
bluesteinInvFactors[input_len][i] = complex(cos, -sin)
}
}
return bluesteinFactors[input_len], bluesteinInvFactors[input_len]
}
func hasBluesteinFactors(idx int) bool {
return bluesteinFactors[idx] != nil
}
// bluesteinFFT returns the FFT calculated using the Bluestein algorithm.
func bluesteinFFT(x []complex128) []complex128 {
lx := len(x)
a := dsputils.ZeroPad(x, dsputils.NextPowerOf2(lx*2-1))
la := len(a)
factors, invFactors := getBluesteinFactors(lx)
for n, v := range x {
a[n] = v * invFactors[n]
}
b := make([]complex128, la)
for i := 0; i < lx; i++ {
b[i] = factors[i]
if i != 0 {
b[la-i] = factors[i]
}
}
r := Convolve(a, b)
for i := 0; i < lx; i++ {
r[i] *= invFactors[i]
}
return r[:lx]
}

224
vendor/github.com/mjibson/go-dsp/fft/fft.go generated vendored Normal file
View File

@ -0,0 +1,224 @@
/*
* Copyright (c) 2011 Matt Jibson <matt.jibson@gmail.com>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
// Package fft provides forward and inverse fast Fourier transform functions.
package fft
import (
"github.com/mjibson/go-dsp/dsputils"
)
// FFTReal returns the forward FFT of the real-valued slice.
func FFTReal(x []float64) []complex128 {
return FFT(dsputils.ToComplex(x))
}
// IFFTReal returns the inverse FFT of the real-valued slice.
func IFFTReal(x []float64) []complex128 {
return IFFT(dsputils.ToComplex(x))
}
// IFFT returns the inverse FFT of the complex-valued slice.
func IFFT(x []complex128) []complex128 {
lx := len(x)
r := make([]complex128, lx)
// Reverse inputs, which is calculated with modulo N, hence x[0] as an outlier
r[0] = x[0]
for i := 1; i < lx; i++ {
r[i] = x[lx-i]
}
r = FFT(r)
N := complex(float64(lx), 0)
for n := range r {
r[n] /= N
}
return r
}
// Convolve returns the convolution of x ∗ y.
func Convolve(x, y []complex128) []complex128 {
if len(x) != len(y) {
panic("arrays not of equal size")
}
fft_x := FFT(x)
fft_y := FFT(y)
r := make([]complex128, len(x))
for i := 0; i < len(r); i++ {
r[i] = fft_x[i] * fft_y[i]
}
return IFFT(r)
}
// FFT returns the forward FFT of the complex-valued slice.
func FFT(x []complex128) []complex128 {
lx := len(x)
// todo: non-hack handling length <= 1 cases
if lx <= 1 {
r := make([]complex128, lx)
copy(r, x)
return r
}
if dsputils.IsPowerOf2(lx) {
return radix2FFT(x)
}
return bluesteinFFT(x)
}
var (
worker_pool_size = 0
)
// SetWorkerPoolSize sets the number of workers during FFT computation on multicore systems.
// If n is 0 (the default), then GOMAXPROCS workers will be created.
func SetWorkerPoolSize(n int) {
if n < 0 {
n = 0
}
worker_pool_size = n
}
// FFT2Real returns the 2-dimensional, forward FFT of the real-valued matrix.
func FFT2Real(x [][]float64) [][]complex128 {
return FFT2(dsputils.ToComplex2(x))
}
// FFT2 returns the 2-dimensional, forward FFT of the complex-valued matrix.
func FFT2(x [][]complex128) [][]complex128 {
return computeFFT2(x, FFT)
}
// IFFT2Real returns the 2-dimensional, inverse FFT of the real-valued matrix.
func IFFT2Real(x [][]float64) [][]complex128 {
return IFFT2(dsputils.ToComplex2(x))
}
// IFFT2 returns the 2-dimensional, inverse FFT of the complex-valued matrix.
func IFFT2(x [][]complex128) [][]complex128 {
return computeFFT2(x, IFFT)
}
func computeFFT2(x [][]complex128, fftFunc func([]complex128) []complex128) [][]complex128 {
rows := len(x)
if rows == 0 {
panic("empty input array")
}
cols := len(x[0])
r := make([][]complex128, rows)
for i := 0; i < rows; i++ {
if len(x[i]) != cols {
panic("ragged input array")
}
r[i] = make([]complex128, cols)
}
for i := 0; i < cols; i++ {
t := make([]complex128, rows)
for j := 0; j < rows; j++ {
t[j] = x[j][i]
}
for n, v := range fftFunc(t) {
r[n][i] = v
}
}
for n, v := range r {
r[n] = fftFunc(v)
}
return r
}
// FFTN returns the forward FFT of the matrix m, computed in all N dimensions.
func FFTN(m *dsputils.Matrix) *dsputils.Matrix {
return computeFFTN(m, FFT)
}
// IFFTN returns the forward FFT of the matrix m, computed in all N dimensions.
func IFFTN(m *dsputils.Matrix) *dsputils.Matrix {
return computeFFTN(m, IFFT)
}
func computeFFTN(m *dsputils.Matrix, fftFunc func([]complex128) []complex128) *dsputils.Matrix {
dims := m.Dimensions()
t := m.Copy()
r := dsputils.MakeEmptyMatrix(dims)
for n := range dims {
dims[n] -= 1
}
for n := range dims {
d := make([]int, len(dims))
copy(d, dims)
d[n] = -1
for {
r.SetDim(fftFunc(t.Dim(d)), d)
if !decrDim(d, dims) {
break
}
}
r, t = t, r
}
return t
}
// decrDim decrements an element of x by 1, skipping all -1s, and wrapping up to d.
// If a value is 0, it will be reset to its correspending value in d, and will carry one from the next non -1 value to the right.
// Returns true if decremented, else false.
func decrDim(x, d []int) bool {
for n, v := range x {
if v == -1 {
continue
} else if v == 0 {
i := n
// find the next element to decrement
for ; i < len(x); i++ {
if x[i] == -1 {
continue
} else if x[i] == 0 {
x[i] = d[i]
} else {
x[i] -= 1
return true
}
}
// no decrement
return false
} else {
x[n] -= 1
return true
}
}
return false
}

199
vendor/github.com/mjibson/go-dsp/fft/radix2.go generated vendored Normal file
View File

@ -0,0 +1,199 @@
/*
* Copyright (c) 2011 Matt Jibson <matt.jibson@gmail.com>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
package fft
import (
"math"
"runtime"
"sync"
)
var (
radix2Lock sync.RWMutex
radix2Factors = map[int][]complex128{
4: {complex(1, 0), complex(0, -1), complex(-1, 0), complex(0, 1)},
}
)
// EnsureRadix2Factors ensures that all radix 2 factors are computed for inputs
// of length input_len. This is used to precompute needed factors for known
// sizes. Generally should only be used for benchmarks.
func EnsureRadix2Factors(input_len int) {
getRadix2Factors(input_len)
}
func getRadix2Factors(input_len int) []complex128 {
radix2Lock.RLock()
if hasRadix2Factors(input_len) {
defer radix2Lock.RUnlock()
return radix2Factors[input_len]
}
radix2Lock.RUnlock()
radix2Lock.Lock()
defer radix2Lock.Unlock()
if !hasRadix2Factors(input_len) {
for i, p := 8, 4; i <= input_len; i, p = i<<1, i {
if radix2Factors[i] == nil {
radix2Factors[i] = make([]complex128, i)
for n, j := 0, 0; n < i; n, j = n+2, j+1 {
radix2Factors[i][n] = radix2Factors[p][j]
}
for n := 1; n < i; n += 2 {
sin, cos := math.Sincos(-2 * math.Pi / float64(i) * float64(n))
radix2Factors[i][n] = complex(cos, sin)
}
}
}
}
return radix2Factors[input_len]
}
func hasRadix2Factors(idx int) bool {
return radix2Factors[idx] != nil
}
type fft_work struct {
start, end int
}
// radix2FFT returns the FFT calculated using the radix-2 DIT Cooley-Tukey algorithm.
func radix2FFT(x []complex128) []complex128 {
lx := len(x)
factors := getRadix2Factors(lx)
t := make([]complex128, lx) // temp
r := reorderData(x)
var blocks, stage, s_2 int
jobs := make(chan *fft_work, lx)
wg := sync.WaitGroup{}
num_workers := worker_pool_size
if (num_workers) == 0 {
num_workers = runtime.GOMAXPROCS(0)
}
idx_diff := lx / num_workers
if idx_diff < 2 {
idx_diff = 2
}
worker := func() {
for work := range jobs {
for nb := work.start; nb < work.end; nb += stage {
if stage != 2 {
for j := 0; j < s_2; j++ {
idx := j + nb
idx2 := idx + s_2
ridx := r[idx]
w_n := r[idx2] * factors[blocks*j]
t[idx] = ridx + w_n
t[idx2] = ridx - w_n
}
} else {
n1 := nb + 1
rn := r[nb]
rn1 := r[n1]
t[nb] = rn + rn1
t[n1] = rn - rn1
}
}
wg.Done()
}
}
for i := 0; i < num_workers; i++ {
go worker()
}
defer close(jobs)
for stage = 2; stage <= lx; stage <<= 1 {
blocks = lx / stage
s_2 = stage / 2
for start, end := 0, stage; ; {
if end-start >= idx_diff || end == lx {
wg.Add(1)
jobs <- &fft_work{start, end}
if end == lx {
break
}
start = end
}
end += stage
}
wg.Wait()
r, t = t, r
}
return r
}
// reorderData returns a copy of x reordered for the DFT.
func reorderData(x []complex128) []complex128 {
lx := uint(len(x))
r := make([]complex128, lx)
s := log2(lx)
var n uint
for ; n < lx; n++ {
r[reverseBits(n, s)] = x[n]
}
return r
}
// log2 returns the log base 2 of v
// from: http://graphics.stanford.edu/~seander/bithacks.html#IntegerLogObvious
func log2(v uint) uint {
var r uint
for v >>= 1; v != 0; v >>= 1 {
r++
}
return r
}
// reverseBits returns the first s bits of v in reverse order
// from: http://graphics.stanford.edu/~seander/bithacks.html#BitReverseObvious
func reverseBits(v, s uint) uint {
var r uint
// Since we aren't reversing all the bits in v (just the first s bits),
// we only need the first bit of v instead of a full copy.
r = v & 1
s--
for v >>= 1; v != 0; v >>= 1 {
r <<= 1
r |= v & 1
s--
}
return r << s
}