Add FFTW bindings, use it in the pipeline

.gitignore

@@ -29,6 +29,7 @@ po/stamp-it
 samples/
 src/*.c
 !src/spek-audio.c
+!src/spek-fft.c
 src/*.o
 src/*.stamp
 src/spek

configure.ac

@@ -24,11 +24,12 @@ AM_PROG_VALAC([0.7.0])
 AC_PROG_INSTALL
 AC_PROG_INTLTOOL([0.35])
 
-pkg_modules="gtk+-2.0 >= 2.14.0 gstreamer-0.10 >= 0.10.17 libavformat libavcodec"
-SPEK_PACKAGES="--pkg gtk+-2.0 --pkg gstreamer-0.10"
+pkg_modules="gtk+-2.0 >= 2.14.0 gstreamer-0.10 >= 0.10.17 libavformat libavcodec fftw3f"
 PKG_CHECK_MODULES(SPEK, [$pkg_modules])
 AC_SUBST(SPEK_CFLAGS)
 AC_SUBST(SPEK_LIBS)
+
+SPEK_PACKAGES="--pkg gtk+-2.0 --pkg gstreamer-0.10"
 AC_SUBST(SPEK_PACKAGES)
 
 AC_CHECK_LIB(m, log10f)

src/Makefile.am

@@ -3,6 +3,7 @@ bin_PROGRAMS = spek
 spek_SOURCES = \
 	spek.vala \
 	spek-audio.c \
+	spek-fft.c \
 	spek-pipeline.vala \
 	spek-ruler.vala \
 	spek-source.vala \
@@ -21,6 +22,7 @@ VALAFLAGS = \
 	--vapidir=$(srcdir)/../vapi \
 	--pkg config \
 	--pkg spek-audio \
+	--pkg spek-fft \
 	@SPEK_PACKAGES@
 
 spek_LDADD = \
@@ -28,4 +30,5 @@ spek_LDADD = \
 	$(IGE_MAC_LIBS)
 
 EXTRA_DIST = \
-	spek-audio.h
+	spek-audio.h \
+	spek-fft.h

src/spek-fft.c

@@ -0,0 +1,55 @@
+/* spek-fft.c
+ *
+ * Copyright (C) 2010  Alexander Kojevnikov <alexander@kojevnikov.com>
+ *
+ * Spek is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Spek is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Spek.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <math.h>
+
+#include "spek-fft.h"
+
+SpekFftPlan * spek_fft_plan_new (gint n) {
+	SpekFftPlan *p = g_new0 (SpekFftPlan, 1);
+	p->input = (gfloat *) fftwf_malloc (sizeof (gfloat) * n);
+	p->output = (gfloat *) fftwf_malloc (sizeof (gfloat) * (n / 2 + 1));
+	p->result = (fftwf_complex *) fftwf_malloc (sizeof (fftwf_complex) * (n / 2 + 1));
+	p->n = n;
+	p->plan = fftwf_plan_dft_r2c_1d (n, p->input, p->result, FFTW_ESTIMATE);
+	return p;
+}
+
+void spek_fft_execute (SpekFftPlan *p) {
+	int i;
+	int bands = p->n / 2 + 1;
+
+	fftwf_execute (p->plan);
+
+	/* Calculate magnitudes */
+	for (i = 0; i < bands; i++) {
+		gfloat val;
+		val = p->result[i][0] * p->result[i][0] + p->result[i][1] * p->result[i][1];
+		val /= p->n * p->n;
+		val = 10.0 * log10f (val);
+		p->output[i] = val;
+	}
+}
+
+void spek_fft_destroy (SpekFftPlan *p) {
+	fftwf_destroy_plan (p->plan);
+	fftwf_free (p->result);
+	fftwf_free (p->output);
+	fftwf_free (p->input);
+	g_free (p);
+}

src/spek-fft.h

@@ -0,0 +1,45 @@
+/* spek-fft.h
+ *
+ * Copyright (C) 2010  Alexander Kojevnikov <alexander@kojevnikov.com>
+ *
+ * Spek is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Spek is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Spek.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef __SPEK_FFT_H__
+#define __SPEK_FFT_H__
+
+#include <glib.h>
+#include <fftw3.h>
+
+typedef struct {
+	/* Internal data */
+	fftwf_plan plan;
+	gint n;
+	fftwf_complex *result;
+
+	/* Exposed properties */
+	gfloat *input;
+	gfloat *output;
+} SpekFftPlan;
+
+/* Allocate buffers and create a new FFT plan */
+SpekFftPlan * spek_fft_plan_new (gint n);
+
+/* Execute the FFT on plan->input */
+void spek_fft_execute (SpekFftPlan *p);
+
+/* Destroy the plan and de-allocate buffers */
+void spek_fft_destroy (SpekFftPlan *p);
+
+#endif

src/spek-pipeline.vala

@@ -77,9 +77,13 @@ namespace Spek {
 		public void start () {
 			int nfft = 2 * bands - 2;
 			var input = new float[nfft];
+			var output = new float[bands];
 			int pos = 0;
 			int frames = 0;
 			int size;
+			// TODO: make it a static class variable and experiment with FFTW_MEASURE.
+			var fft = new Fft.Plan (nfft);
+
 			while ((size = cx.read (this.buffer)) > 0) {
 				uint8 *buffer = (uint8 *) this.buffer;
 				var block_size = cx.width * cx.channels;
@@ -90,7 +94,26 @@ namespace Spek {
 					pos = (pos + 1) % nfft;
 					frames++;
 
-					// TODO
+					// If we have enough frames for an FFT or we
+					// have all frames required for the interval run
+					// an FFT. In the last case we probably take the
+					// FFT of frames that we already handled.
+					if (frames % nfft == 0
+						// TODO: error correction
+//						|| ((spectrum->accumulated_error < GST_SECOND
+//							 && spectrum->num_frames == spectrum->frames_per_interval)
+//							||
+//							(spectrum->accumulated_error >= GST_SECOND
+//							 && spectrum->num_frames - 1 == spectrum->frames_per_interval))
+						) {
+						for (int i = 0; i < nfft; i++) {
+							fft.input[i] = input[(pos + i) % nfft];
+						}
+						fft.execute ();
+						for (int i = 0; i < bands; i++) {
+							output[i] += fft.output[i];
+						}
+					}
 				}
 				assert (size == 0);
 			}

src/spek-spectrogram.vala

@@ -91,6 +91,7 @@ namespace Spek {
 
 			// TODO
 			var pipeline = new Pipeline (file_name, BANDS, samples, THRESHOLD, data_cb);
+			pipeline.start ();
 			print ("\n%s:\n%s\n", file_name, pipeline.description);
 
 			queue_draw ();

vapi/Makefile.am

@@ -1,6 +1,7 @@
 noinst_DATA = \
 	config.vapi \
-	spek-audio.vapi
+	spek-audio.vapi \
+	spek-fft.vapi
 
 EXTRA_DIST = \
 	$(noinst_DATA)

vapi/spek-fft.vapi

@@ -0,0 +1,16 @@
+[CCode (cprefix = "SpekFft", lower_case_cprefix = "spek_fft_", cheader_filename = "spek-fft.h")]
+namespace Spek.Fft {
+	[Compact]
+	[CCode (free_function = "spek_fft_destroy")]
+	public class Plan {
+		[CCode (array_length = false)]
+		public unowned float[] input;
+		[CCode (array_length = false)]
+		public unowned float[] output;
+
+		[CCode (cname = "spek_fft_plan_new")]
+		public Plan (int n);
+		[CCode (cname = "spek_fft_execute")]
+		public void execute ();
+	}
+}