-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfft.odin
111 lines (100 loc) · 2.46 KB
/
fft.odin
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
package main
import "base:intrinsics"
import "base:runtime"
import "core:c"
import "core:fmt"
import "core:math"
import "core:mem"
audioPeriod: int : 3600 / 2
currentLeftChannel: [dynamic]complex64
currentRightChannel: [dynamic]complex64
currentPeriod: int = 0
AudioVizualizerStates :: enum {
Filling,
Copying,
}
vizualizerState: AudioVizualizerStates = .Filling
ResetVizualizerState :: proc() {
// Ensure that the visualizer is reset when we stop
currentPeriod = 0
currentLeftChannel = {}
currentRightChannel = {}
vizualizerState = .Filling
}
AudioProcessFFT :: proc "c" (buffer: rawptr, frames: c.uint) {
context = runtime.default_context()
if buffer == nil || frames == 0 {
return
}
if mem.check_zero_ptr(buffer, int(currentStream.frameCount)) {
fmt.printfln("Buffer is zero")
return
}
if currentPeriod >= audioPeriod {
vizualizerState = .Copying
}
if vizualizerState == .Copying {
return
}
#no_bounds_check {
fs := mem.slice_ptr(cast(^[2]f32)(buffer), int(frames))
for i in 0 ..< int(frames) {
if fs[i][0] != 0 {
append(¤tLeftChannel, fs[i][0])
}
if fs[i][1] != 0 {
append(¤tRightChannel, fs[i][1])
}
}
currentPeriod += int(frames)
}
}
fft :: proc(buffer: []complex64) -> []f32 {
data := buffer
// DFT
#no_bounds_check {
N: uint = uint(len(data))
k: uint = N
n: uint
thetaT: f64 = math.PI / f64(N)
phiT: complex64 = complex(math.cos(thetaT), -math.sin(thetaT))
T: complex64
for k > 1 {
n = k
k >>= 1
phiT = phiT * phiT
T = 1.0
for l: uint = 0; l < k; l = l + 1 {
for a: uint = l; a < N; a += n {
b: uint = a + k
if b < N { // Ensure 'b' is within valid range
t: complex64 = data[a] - data[b]
data[a] += data[b]
data[b] = t * T
}
}
T *= phiT
}
}
// Decimate
m: uint = (uint)(math.log2(f64(N)))
for a: uint = 0; a < N; a = a + 1 {
b: uint = a
// Reverse bits
b = (((b & 0xaaaaaaaa) >> 1) | ((b & 0x55555555) << 1))
b = (((b & 0xcccccccc) >> 2) | ((b & 0x33333333) << 2))
b = (((b & 0xf0f0f0f0) >> 4) | ((b & 0x0f0f0f0f) << 4))
b = (((b & 0xff00ff00) >> 8) | ((b & 0x00ff00ff) << 8))
b = ((b >> 16) | (b << 16)) >> (32 - m)
if b < N && b > a { // Ensure 'b' is within valid range and prevent redundant swaps
t: complex64 = data[a]
data[a] = data[b]
data[b] = t
}
}
// Normalize
f: complex64 = 1.0 / math.sqrt(f32(N))
for i: uint = 0; i < N; i = i + 1 {data[i] *= f}
return transmute([]f32)(data)
}
}