FFT now uses a lookup-table for sin and cos values

On x86, this makes the program approx. 5 times faster.

Makefile

@@ -3,8 +3,8 @@ CFLAGS+=-O3 -Wall -march=native -pedantic -std=c99 -D_POSIX_C_SOURCE=20120607L -
 LIBS=-lm -lpthread -lrt
 
 TARGET=musiclight2
-SOURCE=main.c fft.c utils.c ws2801.c
-DEPS=config.h fft.h utils.h ws2801.h
+SOURCE=main.c fft.c utils.c ws2801.c lut.c
+DEPS=config.h fft.h utils.h ws2801.h lut.h
 
 OBJ=$(patsubst %.c, %.o, $(SOURCE))
 

config.h

@@ -19,7 +19,7 @@
 #define PORT 2703
 
 // FFT transformation parameters
-#define FFT_EXPONENT 10
+#define FFT_EXPONENT 10 // ATTENTION: when you change this, run gen_lut.py with this value as argument
 #define BLOCK_LEN (1 << FFT_EXPONENT) // 2^FFT_EXPONENT
 #define SAMPLE_RATE 44100
 #define DATALEN (BLOCK_LEN / 2)

fft.c

@@ -14,6 +14,7 @@
 
 #include "config.h"
 
+#include "lut.h"
 #include "fft.h"
 
 value_type hanning_buffer[BLOCK_LEN];
@@ -68,7 +69,6 @@ void fft_transform(sample *samples, value_type *resultRe, value_type *resultIm)
   int left, right;
 
   value_type x_left_re, x_left_im, x_right_re, x_right_im;
-  value_type param;
   value_type sinval, cosval;
 
   // re-arrange the input array according to the lookup table
@@ -105,12 +105,10 @@ void fft_transform(sample *samples, value_type *resultRe, value_type *resultIm)
         x_right_re = resultRe[right];
         x_right_im = resultIm[right];
 
-        // precalculate the parameter for sin and cos
-        param = -M_PI * element / (1 << layer);
-
-        // precalculate sinus and cosinus values for param
-        sinval = sin(param);
-        cosval = cos(param);
+        // use lookup table to get sinus and cosinus values for param
+        //param = -M_PI * element / (1 << layer);
+        sinval = lookup_sin(layer, element);
+        cosval = lookup_cos(layer, element);
 
         // combine the values according to a butterfly diagram
         resultRe[left]  = x_right_re + x_left_re * cosval - x_left_im * sinval;

gen_lut.py

@@ -0,0 +1,68 @@
+#!/usr/bin/env python3
+
+import sys
+from math import *
+
+preamble = """// This file was auto-generated using gen_lut.py
+
+#include "lut.h"
+
+"""
+
+postamble = """
+
+value_type lookup_sin(int layer, int element) {
+  return sin_lut[layer][element];
+}
+
+value_type lookup_cos(int layer, int element) {
+  return cos_lut[layer][element];
+}
+"""
+
+if len(sys.argv) < 2:
+  print("Argument required: FFT_EXPONENT")
+  exit(1)
+
+fft_exponent = int(sys.argv[1])
+
+with open("lut.c", "w") as ofile:
+  ofile.write(preamble)
+
+  # generate the sin() lookup table
+  for layer in range(1, fft_exponent+1):
+    num_elements = (1 << layer)
+    ofile.write("value_type sin_lut%i[%i] = {" % (layer, num_elements))
+
+    ofile.write("0")
+
+    for element in range(1, num_elements):
+      ofile.write(", %.10f" % sin(-pi * element / num_elements));
+
+    ofile.write("};\n\n")
+
+  ofile.write("value_type *sin_lut[%i] = {sin_lut1" % fft_exponent);
+
+  for i in range(2, fft_exponent+1):
+    ofile.write(", sin_lut" + str(i));
+  ofile.write("};\n");
+
+  # generate the cos() lookup table
+  for layer in range(1, fft_exponent+1):
+    num_elements = (1 << layer)
+    ofile.write("value_type cos_lut%i[%i] = {" % (layer, num_elements))
+
+    ofile.write("0")
+
+    for element in range(1, num_elements):
+      ofile.write(", %.10f" % cos(-pi * element / num_elements));
+
+    ofile.write("};\n\n")
+
+  ofile.write("value_type *cos_lut[%i] = {cos_lut1" % fft_exponent);
+
+  for i in range(2, fft_exponent+1):
+    ofile.write(", cos_lut" + str(i));
+  ofile.write("};\n");
+
+  ofile.write(postamble)

lut.h

@@ -0,0 +1,9 @@
+#ifndef LUT_H
+#define LUT_H
+
+#include "config.h"
+
+value_type lookup_sin(int layer, int element);
+value_type lookup_cos(int layer, int element);
+
+#endif // LUT_H