-
Notifications
You must be signed in to change notification settings - Fork 2
/
dft.go
77 lines (63 loc) · 1.84 KB
/
dft.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
package gdsp
import (
"math"
"math/cmplx"
)
// DFT performs a discrete Fourier transform on the complex-valued input vector
// and returns the result. Pass true for the forward parameter to perform a forward
// transform and false for an inverse transform.
func DFT(input VectorComplex, forward bool) VectorComplex {
N := len(input)
theta := 2.0 * math.Pi / float64(N)
coeff := 1.0
if forward {
coeff = -1.0
}
output := MakeVectorComplex(0.0, N)
for i := 0; i < N; i++ {
for j := 0; j < N; j++ {
x := float64(i * j)
realP := math.Cos(theta * x)
imagP := coeff * math.Sin(theta*x)
output[i] += input[j] * complex(realP, imagP)
}
if !forward {
output[i] /= complex(float64(N), 0.0)
}
}
return output
}
// FFT performs a discrete Fourier transform on the complex-valued input vector
// using the Cooley-Turkey FFT algorithm. For an inverse FFT, see the IFFT function.
func FFT(input VectorComplex) VectorComplex {
if len(input) == 1 {
return input
}
if len(input)%2 != 0 {
return DFT(input, true)
}
var evenInput VectorComplex
var oddInput VectorComplex
for i := 0; i < len(input); i += 2 {
evenInput = append(evenInput, input[i])
oddInput = append(oddInput, input[i+1])
}
evenDFT := FFT(evenInput)
oddDFT := FFT(oddInput)
for k := 0; k < len(input)/2; k++ {
x := float64(k) / float64(len(input))
t := evenDFT[k]
ec := cmplx.Exp(complex(0.0, -2.0*math.Pi*x))
evenDFT[k] = t + ec*oddDFT[k]
oddDFT[k] = t - ec*oddDFT[k]
}
return append(evenDFT, oddDFT...)
}
// IFFT performs an inverse discrete Fourier transform on the complex-valued input
// vector using the Cooley-Turkey FFT algorithm. For a forward FFT, see the FFT
// function.
func IFFT(input VectorComplex) VectorComplex {
inputConjugate := input.Conj()
fft := FFT(inputConjugate)
return VSDivC(fft.Conj(), complex(float64(len(input)), 0.0))
}