From d284e8161ebd63d7748cef33314b9f4947bfffd1 Mon Sep 17 00:00:00 2001
From: rbn42 <bl100@students.waikato.ac.nz>
Date: Fri, 22 Sep 2017 23:59:45 +1200
Subject: [PATCH] glfft

---
 panon/glsl/fft.glsl            |  87 ++++++++++++++++++++++++++
 panon/visualizer/__init__.py   | 108 ++++-----------------------------
 panon/visualizer/fallback.py   |   3 +-
 panon/visualizer/glfft.py      |  33 ++++++++++
 panon/visualizer/glspectrum.py |  80 ++++++++++++++++++++++++
 panon/visualizer/source.py     |  23 +++++++
 panon/visualizer/spectrum.py   |  83 +++++++++++++++++++++++++
 7 files changed, 320 insertions(+), 97 deletions(-)
 create mode 100644 panon/glsl/fft.glsl
 create mode 100644 panon/visualizer/glfft.py
 create mode 100644 panon/visualizer/glspectrum.py
 create mode 100644 panon/visualizer/source.py
 create mode 100644 panon/visualizer/spectrum.py

diff --git a/panon/glsl/fft.glsl b/panon/glsl/fft.glsl
new file mode 100644
index 0000000..8f7c1a4
--- /dev/null
+++ b/panon/glsl/fft.glsl
@@ -0,0 +1,87 @@
+#version 430
+// based on the algorithm described in http://research.microsoft.com/pubs/70576/tr-2008-62.pdf
+#define SIZE 1024
+#define SIZE2 4 
+#define PI 3.14159265358979323844
+layout(local_size_x = SIZE) in;
+layout(std430) buffer;
+//layout(binding = 0, r32f) writeonly uniform image2D dest_texture;
+layout(binding = 1) readonly buffer Input {
+    float input_data[SIZE*SIZE2];
+};
+layout (std430, binding = 2) writeonly buffer Output {
+    float v2[SIZE*SIZE2];
+};
+
+
+uniform uint real_size;
+
+shared float values1[SIZE*SIZE2];
+shared float values2[SIZE*SIZE2];
+void synchronize()
+{
+    memoryBarrierShared();
+    barrier();
+}
+vec2 getvalue(uint index){
+    float x=values1[index];
+    float y=values2[index];
+    return vec2(x,y);
+}
+
+void setvalue(uint index,vec2 value){
+    values1[index]=value.x;
+    values2[index]=value.y;
+}
+
+void
+fft_pass(int ns, int source,uint i)
+{
+    uint base = (i/ns)*(ns/2);
+    uint offs = i%(ns/2);
+
+    uint i0 = base + offs;
+    uint i1 = i0 + real_size/2;
+    
+    vec2 v0 = getvalue(i0*2+source);
+    vec2 v1 = getvalue(i1*2+source);
+    
+    float a = -2.*PI*float(i)/ns;
+
+    float t_re = cos(a);
+    float t_im = sin(a);
+
+    setvalue(i*2+source ^ 1 , v0 + vec2(dot(vec2(t_re, -t_im), v1), dot(vec2(t_im, t_re), v1)));
+}
+
+void main()
+{
+    uint i = gl_LocalInvocationID.x*SIZE2;
+    for(uint i2=0;i2<SIZE2;i2++){
+        uint index=i+i2;
+        if(index>=real_size)
+            break;
+        setvalue(index*2+0,  vec2(input_data[index], 0.));
+    }
+    synchronize();
+
+    int source = 0;
+
+    for (int n = 2; n <= SIZE; n *= 2) {
+        for(uint i2=0;i2<SIZE2;i2++){
+            uint index=i+i2;
+            if(index>=real_size)
+                break;
+            fft_pass(n, source,index);
+        }   
+        source ^= 1;
+        synchronize();
+    }
+
+    for(uint i2=0;i2<SIZE2;i2++){
+        uint index=i+i2;
+        if(index>=real_size)
+            break;
+        v2[index]=length(getvalue(index*2+source));
+    }
+}
diff --git a/panon/visualizer/__init__.py b/panon/visualizer/__init__.py
index 26f3805..6b056c2 100644
--- a/panon/visualizer/__init__.py
+++ b/panon/visualizer/__init__.py
@@ -1,31 +1,17 @@
 import gi
 gi.require_version('Gtk', '3.0')
-from gi.repository import Gtk, GObject, Gdk
+from gi.repository import Gtk, GObject
 import cairo
-import pyaudio
-import numpy as np
 from .. import helper
 from .. import config
 from .fallback import VisualizerCairo
 from .opengl import VisualizerGL
+from .source import Source
+from .spectrum import Spectrum
 from queue import Queue
 from threading import Thread
 
 
-def record_pyaudio(fps, channel_count, sample_rate):
-    p = pyaudio.PyAudio()
-    stream = p.open(format=pyaudio.paInt16,
-                    channels=channel_count,
-                    rate=sample_rate,
-                    input=True)
-    stop = False
-    while not stop:
-        size = stream.get_read_available()
-        stop = yield np.fromstring(stream.read(size), 'int16')
-    stream.close()
-    yield
-
-
 class Visualizer(Gtk.EventBox):
     stop = False
     stop_gen_data = False
@@ -41,19 +27,20 @@ def __init__(self, background_color, fps=60, channel_count=2, sample_rate=44100,
         super(Visualizer, self).__init__()
         self.sample_rate = sample_rate
         self.background_color = helper.color(background_color)
-        self.history = [[]] * 8
         self.data_queue = Queue(3)
-        self.min_sample = 10
-        self.max_sample = self.min_sample
+
         self.padding = padding
-        self.buffer_size = sample_rate // fps * channel_count
         self.fps = fps
         self.channel_count = channel_count
         self.sample_rate = sample_rate
-        self.sample = record_pyaudio(fps, channel_count, sample_rate)
+        self.sample = Source(channel_count, sample_rate)
+        buffer_size = sample_rate // fps * channel_count
+        self.spectrum = Spectrum(
+            self.sample, buffer_size, config.visualizer_decay)
         GObject.timeout_add(1000 // fps, self.tick)
 
-        self.use_opengl=use_opengl
+        self.use_opengl = use_opengl
+
         if use_opengl:
             Thread(target=self.run).start()
             self.da = VisualizerGL(self.getData)
@@ -80,12 +67,11 @@ def do_button_release_event(self, widget, event=None):
         if event and event.button == 1:
             if not self.stop:
                 self.da.stop()
-                self.sample.send(True)
+                self.sample.stop()
                 self.stop = True
             else:
                 self.stop = False
-                self.sample = record_pyaudio(
-                    self.fps, self.channel_count, self.sample_rate)
+                self.sample.start()
                 self.da.start()
             return True
         else:
@@ -99,72 +85,4 @@ def getData(self):
         if self.use_opengl:
             return self.data_queue.get()
         else:
-            return self.__getData()
-
-    def __getData(self):
-        #fft = np.absolute(np.fft.rfft(data, n=len(data)))/len(data)
-
-        data = next(self.sample)
-        self.history.append(data)
-        if sum([len(d) for d in self.history[1:]]) > self.buffer_size * 8:
-            self.history.pop(0)
-
-        data_history = np.concatenate(self.history)
-        fft_freq = []
-
-        def fun(start, end,  rel):
-            size = self.buffer_size
-            if rel > 20:
-                start, end = int(start), int(end)
-                rel = int(rel)
-                d = data_history[-size * rel:].reshape((rel, size))
-                d = np.mean(d, axis=0)
-            else:
-                start = int(start * rel)
-                end = int(end * rel)
-                size = int(size * rel)
-                d = data_history[-size:]
-
-            fft = np.absolute(np.fft.rfft(d, n=size))
-            end = min(len(fft) // 2, end)
-            fft_freq.insert(0, fft[start:end])
-            fft_freq.append(fft[len(fft) - end:len(fft) - start])
-        # higher resolution and latency for lower frequency
-
-        sections = 8
-        r = 0.6
-        rels = 8 * r**np.arange(sections)
-        start = 0
-        sections = []
-        for rel, freq_width in zip(rels, len(data) * 1 / rels / sum(1 / rels) // 4):
-            if rel > 2:
-                freq_width *= rel
-                pass
-            sections.append((start, start + freq_width, rel))
-            start += freq_width
-        sections.reverse()
-        for start, end, rel in sections:
-            #fun(start, end, rel)
-            pass
-
-        #fun(400, len(data),  0.3)
-        #fun(300, 400,  0.5)
-        #fun(200,300 ,  0.75)
-        #fun(150, 200,  1)
-        fun(110, 150,  2)
-        fun(80, 110,  3)
-        fun(50, 80,  4)
-        fun(30, 50,  5)
-        fun(10, 30,  6)
-        fun(0, 10,  8)
-
-        fft = np.concatenate(fft_freq)
-
-        exp = 2
-        retain = (1 - config.visualizer_decay)**exp
-        decay = 1 - retain
-
-        vol = self.min_sample + np.mean(fft ** exp)
-        self.max_sample = self.max_sample * retain + vol * decay
-        bins = fft / self.max_sample ** (1 / exp)
-        return bins
+            return self.spectrum.getData()
diff --git a/panon/visualizer/fallback.py b/panon/visualizer/fallback.py
index 4c85115..e566e5f 100644
--- a/panon/visualizer/fallback.py
+++ b/panon/visualizer/fallback.py
@@ -1,8 +1,7 @@
 import cairo
 import gi
 gi.require_version('Gtk', '3.0')
-from gi.repository import Gtk, GObject, Gdk
-
+from gi.repository import Gtk, Gdk
 from .. import helper
 
 
diff --git a/panon/visualizer/glfft.py b/panon/visualizer/glfft.py
new file mode 100644
index 0000000..6af22a6
--- /dev/null
+++ b/panon/visualizer/glfft.py
@@ -0,0 +1,33 @@
+from .. import glsl
+import ModernGL
+
+
+class GLFFT:
+    def __init__(self):
+        ctx = ModernGL.create_standalone_context()
+        compute_shader = ctx.compute_shader(glsl.load('fft.glsl'))
+        size = 1024 * 4
+        empty = b'\00' * size * 4
+        buf1 = ctx.buffer(empty)
+        empty = b'\00' * size * 4
+        buf2 = ctx.buffer(empty)
+        buf1.bind_to_storage_buffer(1)
+        buf2.bind_to_storage_buffer(2)
+        #compute_shader.uniforms['mul'].value = 100.0
+
+        self.compute_shader = compute_shader
+        #self.compute_shader.uniforms['mul'].value = 100.0
+        self.ctx = ctx
+        self.buf1 = buf1
+        self.buf2 = buf2
+
+    def compute(self, data):
+        data=data[-1024*4*4:]
+        self.compute_shader.uniforms['real_size'].value = len(data)//4
+        #self.compute_shader.uniforms['mul'].value = 100.0
+        self.buf1.write(data)
+        self.compute_shader.run()
+        return self.buf2.read()[:len(data)]
+
+    def destroy(self):
+        self.ctx.release()
diff --git a/panon/visualizer/glspectrum.py b/panon/visualizer/glspectrum.py
new file mode 100644
index 0000000..02ed182
--- /dev/null
+++ b/panon/visualizer/glspectrum.py
@@ -0,0 +1,80 @@
+import numpy as np
+
+
+class GLSpectrum:
+    def __init__(self, sample, buffer_size, decay):
+        self.sample = sample
+        self.decay = decay
+        self.history = [[]] * 8
+        self.buffer_size = buffer_size
+        self.min_sample = 10
+        self.max_sample = self.min_sample
+
+        from .glfft import GLFFT
+        self.glfft = GLFFT()
+
+    def getData(self):
+        data = self.sample.read()
+        data = np.fromstring(data, 'int16')
+
+        self.history.append(data)
+        if sum([len(d) for d in self.history[1:]]) > self.buffer_size * 8:
+            self.history.pop(0)
+
+        data_history = np.concatenate(self.history)
+        fft_freq = []
+
+        def fun(start, end,  rel):
+            size = self.buffer_size
+            if rel > 20:
+                start, end = int(start), int(end)
+                rel = int(rel)
+                d = data_history[-size * rel:].reshape((rel, size))
+                d = np.mean(d, axis=0)
+            else:
+                start = int(start * rel)
+                end = int(end * rel)
+                size = int(size * rel)
+                d = data_history[-size:]
+
+            fft = self.glfft.compute(d.astype('float32').tobytes())
+            fft = np.frombuffer(fft, dtype='float32')
+            end = min(len(fft) // 2, end)
+            fft_freq.insert(0, fft[start:end])
+            fft_freq.append(fft[len(fft) - end:len(fft) - start])
+        # higher resolution and latency for lower frequency
+
+        sections = 8
+        r = 0.6
+        rels = 8 * r**np.arange(sections)
+        start = 0
+        sections = []
+        for rel, freq_width in zip(rels, len(data) * 1 / rels / sum(1 / rels) // 4):
+            if rel > 2:
+                freq_width *= rel
+                pass
+            sections.append((start, start + freq_width, rel))
+            start += freq_width
+        sections.reverse()
+        for start, end, rel in sections:
+            #fun(start, end, rel)
+            pass
+
+        fun(110, 150, 2)
+#        fun(0, 110, 3)
+        fun(80, 110, 3)
+        fun(50, 80, 4)
+        fun(30, 50, 5)
+        fun(10, 30, 6)
+        fun(0, 10, 8)
+
+        fft = np.concatenate(fft_freq)
+
+        exp = 2
+        retain = (1 - self.decay)**exp
+        decay = 1 - retain
+
+        vol = self.min_sample + np.mean(fft ** exp)
+        self.max_sample = self.max_sample * retain + vol * decay
+        bins = fft / self.max_sample ** (1 / exp)
+        return bins
diff --git a/panon/visualizer/source.py b/panon/visualizer/source.py
new file mode 100644
index 0000000..305e464
--- /dev/null
+++ b/panon/visualizer/source.py
@@ -0,0 +1,23 @@
+import pyaudio
+
+
+class Source:
+    def __init__(self,  channel_count, sample_rate):
+        self.channel_count = channel_count
+        self.sample_rate = sample_rate
+
+        self.start()
+
+    def read(self):
+        size = self.stream.get_read_available()
+        return self.stream.read(size)
+
+    def stop(self):
+        self.stream.close()
+
+    def start(self):
+        p = pyaudio.PyAudio()
+        self.stream = p.open(format=pyaudio.paInt16,
+                             channels=self.channel_count,
+                             rate=self.sample_rate,
+                             input=True)
diff --git a/panon/visualizer/spectrum.py b/panon/visualizer/spectrum.py
new file mode 100644
index 0000000..429152c
--- /dev/null
+++ b/panon/visualizer/spectrum.py
@@ -0,0 +1,83 @@
+import numpy as np
+
+
+class Spectrum:
+    def __init__(self, sample, buffer_size, decay):
+        self.sample = sample
+        self.decay = decay
+        self.history = [[]] * 8
+        self.buffer_size = buffer_size
+        self.min_sample = 10
+        self.max_sample = self.min_sample
+
+        if False:
+            from .glfft import GLFFT
+            self.glfft = GLFFT(512)
+
+    def getData(self):
+        data = self.sample.read()
+        data = np.fromstring(data, 'int16')
+
+        self.history.append(data)
+        if sum([len(d) for d in self.history[1:]]) > self.buffer_size * 8:
+            self.history.pop(0)
+
+        data_history = np.concatenate(self.history)
+        fft_freq = []
+
+        def fun(start, end,  rel):
+            size = self.buffer_size
+            if rel > 20:
+                start, end = int(start), int(end)
+                rel = int(rel)
+                d = data_history[-size * rel:].reshape((rel, size))
+                d = np.mean(d, axis=0)
+            else:
+                start = int(start * rel)
+                end = int(end * rel)
+                size = int(size * rel)
+                d = data_history[-size:]
+
+            if False:
+                fft = self.glfft.compute(d.astype('float32')[:512].tobytes())
+                fft = np.frombuffer(fft, dtype='float32')
+            else:
+                fft = np.absolute(np.fft.rfft(d, n=size))
+            end = min(len(fft) // 2, end)
+            fft_freq.insert(0, fft[start:end])
+            fft_freq.append(fft[len(fft) - end:len(fft) - start])
+        # higher resolution and latency for lower frequency
+
+        sections = 8
+        r = 0.6
+        rels = 8 * r**np.arange(sections)
+        start = 0
+        sections = []
+        for rel, freq_width in zip(rels, len(data) * 1 / rels / sum(1 / rels) // 4):
+            if rel > 2:
+                freq_width *= rel
+                pass
+            sections.append((start, start + freq_width, rel))
+            start += freq_width
+        sections.reverse()
+        for start, end, rel in sections:
+            #fun(start, end, rel)
+            pass
+
+        fun(110, 150, 2)
+        fun(80, 110, 3)
+        fun(50, 80, 4)
+        fun(30, 50, 5)
+        fun(10, 30, 6)
+        fun(0, 10, 8)
+
+        fft = np.concatenate(fft_freq)
+
+        exp = 2
+        retain = (1 - self.decay)**exp
+        decay = 1 - retain
+
+        vol = self.min_sample + np.mean(fft ** exp)
+        self.max_sample = self.max_sample * retain + vol * decay
+        bins = fft / self.max_sample ** (1 / exp)
+        return bins