-
Notifications
You must be signed in to change notification settings - Fork 6
/
fft.go
176 lines (161 loc) · 3.31 KB
/
fft.go
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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
// Copyright 2018 The ZikiChombo Authors. All rights reserved. Use of this source
// code is governed by a license that can be found in the License file.
package fft
import (
"math"
"math/cmplx"
)
// R2Size returns the least power of 2 not less than
// n, the size for which fourier transforms use the
// fast radix-2 algorithm.
func R2Size(n int) int {
j := 1
for j < n {
j *= 2
}
return j
}
// Do performs an in-place FFT on d, if possible.
//
// As fft implementations require different slice
// lengths and capacities depending on the length
// of d, the operation Do will only be in-place if
// d has len and cap equal to
//
// New(len(d)).Win(d)
func Do(d []complex128) {
// XXX cap dimensions can break this.
t := New(len(d))
wd := t.Win(d)
t.Do(wd)
if &d[0] != &wd[0] {
copy(d, wd)
}
}
// To performs a FFT on src placing the results
// in dst if dst has approprioate capacity returning
// either dst or the results in a new slice.
func To(dst, src []complex128) []complex128 {
t := New(len(src))
dst = t.ensureSrc(dst)
t.To(dst, src)
return dst
}
// Inv is like Do but for inverse transform
func Inv(d []complex128) {
t := New(len(d))
t.Inv(d)
}
// InvTo is like To but for inverse transform
func InvTo(dst, src []complex128) []complex128 {
t := New(len(src))
dst = t.ensureSrc(dst)
t.InvTo(dst, src)
return dst
}
// Ny gives the index of the first frequency bin at or above the Nyquist limit
// of a FT of size n.
func Ny(n int) int {
m := n / 2
if n&1 != 0 {
m++
}
return m
}
// Dilate changes the frequency basis by a factor of n/m.
// This is a pitch shift relative to the quantized
// frequency domain, hence there is informaion loss,
// and some values may just be clobbered...
func Dilate(d []complex128, p, q int) {
if p == q {
return
}
n := len(d)
h := n/2 - 1
if p > q {
for i := h; i >= 1; i-- {
dst := (i * p) / q
if dst > h {
d[i] = 0i
continue
}
v := d[i]
d[i] = 0i
d[dst] += v
}
for i := n - 1; i > h; i-- {
dst := (i * p) / q
if dst >= n {
d[i] = 0i
continue
}
v := d[i]
d[i] = 0i
d[dst] += v
}
return
}
for i := 1; i <= h; i++ {
dst := (i * p) / q
if dst > h {
d[i] = 0i
continue
}
v := d[i]
d[i] = 0i
d[dst] += v
}
for i := h + 1; i < n; i++ {
dst := (i * p) / q
if dst > n {
d[i] = 0i
continue
}
v := d[i]
d[i] = 0i
d[dst] += v
}
}
// older stuff
// radix 2, in place recursive
func radix2(src []float64, dst []complex128, N, stride int, inv float64) {
if N == 1 {
dst[0] = complex(src[0], 0.0)
return
}
h := N / 2
radix2(src, dst[:h], h, stride*2, inv)
radix2(src[stride:], dst[h:], h, stride*2, inv)
// recombine
for i := 0; i < h; i++ {
ed := dst[i]
od := dst[i+h]
exp := complex(0, inv*math.Pi*2.0*float64(i)/float64(N))
c := cmplx.Exp(exp)
dst[i] = ed + c*od
dst[i+h] = ed - c*od
}
}
// P is prime.
// TBD(wsc) Rader?
func slowft(src []float64, dst []complex128, P, stride int, inv float64) {
if P == 1 {
dst[0] = complex(src[0], 0)
return
}
var a, s, fi float64
var c complex128
lim := P * stride
n := float64(len(dst))
A := inv * 2.0 * math.Pi / n
for i := 0; i < len(dst); i++ {
fi = float64(i)
c = 0i
for j := 0; j < lim; j += stride {
a = A * float64(j/stride) * fi
s = src[j]
c += complex(s*math.Cos(a), s*math.Sin(a))
}
dst[i] = c
}
}