1 /*
2 * fft.c
3 * Copyright 2011 John Lindgren
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions, and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions, and the following disclaimer in the documentation
13 * provided with the distribution.
14 *
15 * This software is provided "as is" and without any warranty, express or
16 * implied. In no event shall the authors be liable for any damages arising from
17 * the use of this software.
18 */
19
20 // this version has few changes compared to the original audacious fft.c
21 // please find the original file in audacious
22
23 #ifdef HAVE_CONFIG_H
24 # include <config.h>
25 #endif
26 #include "deadbeef.h"
27 #include <math.h>
28 #include <complex.h>
29
30 #if __FreeBSD_version < 902000
31 # define cexpf(x) (expf(crealf(x))*(cosf(cimagf(x))+sinf(cimagf(x))*I))
32 #endif
33
34 #define N (DDB_FREQ_BANDS * 2)
35
36 static float hamming[N]; /* hamming window, scaled to sum to 1 */
37 static int reversed[N]; /* bit-reversal table */
38 static float complex roots[N / 2]; /* N-th roots of unity */
39 static int generated = 0;
40 static float LOGN; /* log N (base 2) */
41
42 /* Reverse the order of the lowest LOGN bits in an integer. */
43
bit_reverse(int x)44 static int bit_reverse (int x)
45 {
46 int y = 0;
47
48 for (int n = LOGN; n --; )
49 {
50 y = (y << 1) | (x & 1);
51 x >>= 1;
52 }
53
54 return y;
55 }
56
57 #ifndef HAVE_LOG2
log2(float x)58 static inline float log2(float x) {return (float)log(x)/M_LN2;}
59 #endif
60
61 /* Generate lookup tables. */
62
generate_tables(void)63 static void generate_tables (void)
64 {
65 if (generated)
66 return;
67
68 LOGN = log2(N);
69 for (int n = 0; n < N; n ++)
70 hamming[n] = 1 - 0.85 * cosf (2 * M_PI * n / N);
71 for (int n = 0; n < N; n ++)
72 reversed[n] = bit_reverse (n);
73 for (int n = 0; n < N / 2; n ++)
74 roots[n] = cexpf (2 * M_PI * I * n / N);
75
76 generated = 1;
77 }
78
do_fft(float complex a[N])79 static void do_fft (float complex a[N])
80 {
81 int half = 1; /* (2^s)/2 */
82 int inv = N / 2; /* N/(2^s) */
83
84 /* loop through steps */
85 while (inv)
86 {
87 /* loop through groups */
88 for (int g = 0; g < N; g += half << 1)
89 {
90 /* loop through butterflies */
91 for (int b = 0, r = 0; b < half; b ++, r += inv)
92 {
93 float complex even = a[g + b];
94 float complex odd = roots[r] * a[g + half + b];
95 a[g + b] = even + odd;
96 a[g + half + b] = even - odd;
97 }
98 }
99
100 half <<= 1;
101 inv >>= 1;
102 }
103 }
104
105 void
calc_freq(float * data,float * freq)106 calc_freq (float *data, float *freq) {
107 generate_tables ();
108
109 // fft code shamelessly stolen from audacious
110 // thanks, John
111 float complex a[N];
112 for (int n = 0; n < N; n ++) {
113 a[reversed[n]] = data[n] * hamming[n];
114 }
115 do_fft(a);
116
117 for (int n = 0; n < N / 2 - 1; n ++)
118 freq[n] = 2 * cabsf (a[1 + n]) / N;
119 freq[N / 2 - 1] = cabsf(a[N / 2]) / N;
120 }
121