From 145abd9383de744d3c7e477a7a6cc3d4e530d815 Mon Sep 17 00:00:00 2001 From: Louis Dalibard Date: Thu, 16 Jan 2025 21:47:06 +0100 Subject: [PATCH] yes --- vendor/github.com/mjibson/go-dsp/LICENSE | 13 + .../mjibson/go-dsp/dsputils/compare.go | 96 ++++++++ .../mjibson/go-dsp/dsputils/dsputils.go | 115 +++++++++ .../mjibson/go-dsp/dsputils/matrix.go | 216 +++++++++++++++++ .../mjibson/go-dsp/fft/bluestein.go | 94 ++++++++ vendor/github.com/mjibson/go-dsp/fft/fft.go | 224 ++++++++++++++++++ .../github.com/mjibson/go-dsp/fft/radix2.go | 199 ++++++++++++++++ 7 files changed, 957 insertions(+) create mode 100644 vendor/github.com/mjibson/go-dsp/LICENSE create mode 100644 vendor/github.com/mjibson/go-dsp/dsputils/compare.go create mode 100644 vendor/github.com/mjibson/go-dsp/dsputils/dsputils.go create mode 100644 vendor/github.com/mjibson/go-dsp/dsputils/matrix.go create mode 100644 vendor/github.com/mjibson/go-dsp/fft/bluestein.go create mode 100644 vendor/github.com/mjibson/go-dsp/fft/fft.go create mode 100644 vendor/github.com/mjibson/go-dsp/fft/radix2.go diff --git a/vendor/github.com/mjibson/go-dsp/LICENSE b/vendor/github.com/mjibson/go-dsp/LICENSE new file mode 100644 index 0000000..d412027 --- /dev/null +++ b/vendor/github.com/mjibson/go-dsp/LICENSE @@ -0,0 +1,13 @@ +Copyright (c) 2011 Matt Jibson + +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. diff --git a/vendor/github.com/mjibson/go-dsp/dsputils/compare.go b/vendor/github.com/mjibson/go-dsp/dsputils/compare.go new file mode 100644 index 0000000..8cab075 --- /dev/null +++ b/vendor/github.com/mjibson/go-dsp/dsputils/compare.go @@ -0,0 +1,96 @@ +/* + * Copyright (c) 2011 Matt Jibson + * + * 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 +} diff --git a/vendor/github.com/mjibson/go-dsp/dsputils/dsputils.go b/vendor/github.com/mjibson/go-dsp/dsputils/dsputils.go new file mode 100644 index 0000000..8ea3be7 --- /dev/null +++ b/vendor/github.com/mjibson/go-dsp/dsputils/dsputils.go @@ -0,0 +1,115 @@ +/* + * Copyright (c) 2011 Matt Jibson + * + * 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 +} diff --git a/vendor/github.com/mjibson/go-dsp/dsputils/matrix.go b/vendor/github.com/mjibson/go-dsp/dsputils/matrix.go new file mode 100644 index 0000000..ac07855 --- /dev/null +++ b/vendor/github.com/mjibson/go-dsp/dsputils/matrix.go @@ -0,0 +1,216 @@ +/* + * Copyright (c) 2011 Matt Jibson + * + * 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) +} diff --git a/vendor/github.com/mjibson/go-dsp/fft/bluestein.go b/vendor/github.com/mjibson/go-dsp/fft/bluestein.go new file mode 100644 index 0000000..998fb98 --- /dev/null +++ b/vendor/github.com/mjibson/go-dsp/fft/bluestein.go @@ -0,0 +1,94 @@ +/* + * Copyright (c) 2011 Matt Jibson + * + * 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] +} diff --git a/vendor/github.com/mjibson/go-dsp/fft/fft.go b/vendor/github.com/mjibson/go-dsp/fft/fft.go new file mode 100644 index 0000000..fee48f9 --- /dev/null +++ b/vendor/github.com/mjibson/go-dsp/fft/fft.go @@ -0,0 +1,224 @@ +/* + * Copyright (c) 2011 Matt Jibson + * + * 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 +} diff --git a/vendor/github.com/mjibson/go-dsp/fft/radix2.go b/vendor/github.com/mjibson/go-dsp/fft/radix2.go new file mode 100644 index 0000000..34be593 --- /dev/null +++ b/vendor/github.com/mjibson/go-dsp/fft/radix2.go @@ -0,0 +1,199 @@ +/* + * Copyright (c) 2011 Matt Jibson + * + * 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 +}