1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
package detect
import "math"
// NextPowerOfTwo returns the smallest power of two at or above n.
func NextPowerOfTwo(n int) int {
p := 1
for p < n {
p <<= 1
}
return p
}
// Transform computes an in-place radix-2 Cooley-Tukey FFT (or its inverse) of a complex signal
// given as parallel real and imaginary slices. len(re) must be a power of two.
//
// Twiddle factors are computed once per (stage, k) and reused across every block at that
// stage, rather than once per butterfly. The Kotlin version recomputes them per butterfly --
// (n/2)*log2(n) sin+cos calls against this version's n-1 -- which measured as not mattering
// there, because the JVM compiles those to hardware intrinsics and memory traffic dominates
// anyway. Done properly here because there is no reason not to.
func Transform(re, im []float64, inverse bool) {
n := len(re)
if n <= 1 {
return
}
if n&(n-1) != 0 {
panic("detect: FFT length must be a power of two")
}
bitReverse(re, im)
sign := -1.0
if inverse {
sign = 1.0
}
for size := 2; size <= n; size <<= 1 {
half := size / 2
step := sign * 2 * math.Pi / float64(size)
for k := 0; k < half; k++ {
wRe := math.Cos(step * float64(k))
wIm := math.Sin(step * float64(k))
for start := 0; start < n; start += size {
i, j := start+k, start+k+half
oddRe := re[j]*wRe - im[j]*wIm
oddIm := re[j]*wIm + im[j]*wRe
evenRe, evenIm := re[i], im[i]
re[i], im[i] = evenRe+oddRe, evenIm+oddIm
re[j], im[j] = evenRe-oddRe, evenIm-oddIm
}
}
}
if inverse {
for i := range re {
re[i] /= float64(n)
im[i] /= float64(n)
}
}
}
func bitReverse(re, im []float64) {
n := len(re)
j := 0
for i := 0; i < n; i++ {
if i < j {
re[i], re[j] = re[j], re[i]
im[i], im[j] = im[j], im[i]
}
bit := n >> 1
for bit != 0 && j&bit != 0 {
j ^= bit
bit >>= 1
}
j ^= bit
}
}
|
