From 593e1cbe3470f409c19a7f06f5d10c6f5677361a Mon Sep 17 00:00:00 2001
From: Ludovic Pouzenc <ludovic@pouzenc.fr>
Date: Sat, 16 Jun 2012 21:16:46 +0000
Subject: Bon, calcul du niveau sonore dans la plage 200 à 2000Hz en reprennant
 les choses calmement. Il n'est pas impossible que les valeurs soient bonnes
 :P
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

git-svn-id: file:///var/svn/2012-violon-leds/trunk@19 6be1fa4d-33ac-4c33-becc-79fcb3794bb6
---
 tests/test6/Makefile  |  30 +++++++
 tests/test6/compil.sh |   2 -
 tests/test6/compute.c | 138 +++++++++++++++++++++++++++++
 tests/test6/compute.h |  11 +++
 tests/test6/fft.c     | 241 ++++++++++++++++++++++++++++++++++++++++++++++++++
 tests/test6/fft.h     |  19 ++++
 tests/test6/test.raw  |   1 +
 tests/test6/test6.c   |  19 +++-
 8 files changed, 456 insertions(+), 5 deletions(-)
 create mode 100644 tests/test6/Makefile
 delete mode 100755 tests/test6/compil.sh
 create mode 100644 tests/test6/compute.c
 create mode 100644 tests/test6/compute.h
 create mode 100644 tests/test6/fft.c
 create mode 100644 tests/test6/fft.h
 create mode 120000 tests/test6/test.raw

diff --git a/tests/test6/Makefile b/tests/test6/Makefile
new file mode 100644
index 0000000..560ff03
--- /dev/null
+++ b/tests/test6/Makefile
@@ -0,0 +1,30 @@
+CC=gcc
+CFLAGS=-W -Wall -Werror -Wno-error=unused-parameter -g
+LDFLAGS=-Werror -g
+EXEC=test6
+
+#CFLAGS+=$(shell pkg-config --cflags gtk+-2.0 gthread-2.0 libpulse)
+#LDFLAGS+=$(shell pkg-config --libs gtk+-2.0 gthread-2.0 libpulse)
+LDFLAGS+=-lm
+
+SRC= $(wildcard *.c)
+OBJ= $(SRC:.c=.o)
+
+all: $(EXEC)
+
+$(EXEC): $(OBJ)
+	$(CC) -o $@ $^ $(LDFLAGS)
+
+#main.o: hello.h
+
+%.o: %.c
+	$(CC) -o $@ -c $< $(CFLAGS)
+
+.PHONY: clean mrproper
+
+clean:
+	rm *.o
+
+mrproper: clean
+	rm $(EXEC)
+
diff --git a/tests/test6/compil.sh b/tests/test6/compil.sh
deleted file mode 100755
index 2e49576..0000000
--- a/tests/test6/compil.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-#!/bin/bash -xe
-gcc -Wall -g -o test6 test6.c -lm
diff --git a/tests/test6/compute.c b/tests/test6/compute.c
new file mode 100644
index 0000000..017416b
--- /dev/null
+++ b/tests/test6/compute.c
@@ -0,0 +1,138 @@
+#include "compute.h"
+
+#include "fft.h"
+#include <math.h>
+
+#define MIN_SAMPLES 256
+#define MAX_SAMPLES 2048
+
+//#define MAX(a,b) (a>b?a:b)
+//#define MIN(a,b) (a<b?a:b)
+
+//static inline float todB_a(const float *x);
+void compute_spectrum(float *data, int width, float output[PSHalf]);
+
+
+float compute_level(const float *data, size_t nsamples, int rate) {
+
+	size_t i; 
+	float input[MAX_SAMPLES], pwrspec[PSHalf];
+	float value;
+	int f, min_f_index, max_f_index;
+
+	if (nsamples >= MAX_SAMPLES) {
+		printf("WARN : nsamples >= MAX_SAMPLES : %i >= %i\n", nsamples, MAX_SAMPLES);
+		nsamples=MAX_SAMPLES;
+	}
+	if (nsamples < MIN_SAMPLES) {
+		printf("WARN : nsamples < MIN_SAMPLES : %i >= %i\n", nsamples, MIN_SAMPLES);
+		// Replicate with symmetry the sound to obtain an input buffer of the minimal len
+		for (i=0;i<MIN_SAMPLES;i++) {
+			if ( (i/nsamples)%2==1 )
+				input[i]=data[i]; // First channel only
+			else
+				input[i]=data[nsamples-i-1];
+		}
+		nsamples=MIN_SAMPLES;
+	} else {
+		for (i=0;i<nsamples;i++) {
+			input[i]=data[i]; // First channel only
+		}
+	}
+
+	compute_spectrum(input, nsamples, pwrspec);
+
+	// Compute the mean power for 200Hz to 2000Hz band
+	min_f_index=((float)PSHalf)*200.f/(((float)rate)/2.f);
+	max_f_index=((float)PSHalf)*2000.f/(((float)rate)/2.f);
+
+	value=0.f;
+	for (f=min_f_index;f<=max_f_index;f++) {
+		value+=pwrspec[f];
+	}
+	// Mean value
+	value=value/(max_f_index-min_f_index+1);
+
+	return value;
+}
+/*
+static inline float todB_a(const float *x){
+  return (float)((*(int32_t *)x)&0x7fffffff) * 7.17711438e-7f -764.6161886f;
+}
+*/
+// Adapted from Audacity 
+void compute_spectrum(float *data, int width, float output[PSHalf]) {
+
+	int i, start, windows;
+	float temp;
+	float in[PSNumS];
+	float out[PSHalf];
+	float processed[PSHalf]={0.0f};
+
+	start = 0;
+	windows = 0;
+	while (start + PSNumS <= width) {
+		// Windowing : Hanning
+		for (i=0; i<PSNumS; i++)
+			in[i] = data[start+i] *(0.50-0.50*cos(2*M_PI*i/(PSNumS-1)));
+
+		// Returns only the real part of the result
+		PowerSpectrum(in, out);
+
+		for (i=0; i<PSHalf; i++)
+			processed[i] += out[i];
+
+		start += PSHalf;
+		windows++;
+	}
+	// Convert to decibels
+	// But do it safely; -Inf is nobody's friend
+	for (i = 0; i < PSHalf; i++){
+		temp=(processed[i] / PSNumS / windows);
+		if (temp > 0.0)
+			output[i] = 10*log10(temp);
+		else
+			output[i] = 0;
+	}
+}
+
+void audio2hsv_1(int audio_level, int *light_h, int *light_s, int *light_v) {
+	// Dummy code
+	*light_h=-audio_level;
+	*light_s=audio_level;
+	*light_v=65535;
+}
+
+		
+void hsv2rgb(int h, int s, int v, int *r, int *g, int *b) {
+   /*
+    * Purpose:
+    * Convert HSV values to RGB values
+    * All values are in the range [0..65535]
+    */
+   float F, M, N, K;
+   int   I;
+   
+   if ( s == 0 ) {
+      /* 
+       * Achromatic case, set level of grey 
+       */
+      *r = v;
+      *g = v;
+      *b = v;
+   } else {
+      I = (int) h/(65535/6);	/* should be in the range 0..5 */
+      F = h - I;		/* fractional part */
+
+      M = v * (1 - s);
+      N = v * (1 - s * F);
+      K = v * (1 - s * (1 - F));
+
+      if (I == 0) { *r = v; *g = K; *b = M; }
+      if (I == 1) { *r = N; *g = v; *b = M; }
+      if (I == 2) { *r = M; *g = v; *b = K; }
+      if (I == 3) { *r = M; *g = N; *b = v; }
+      if (I == 4) { *r = K; *g = M; *b = v; }
+      if (I == 5) { *r = v; *g = M; *b = N; }
+   }
+}
diff --git a/tests/test6/compute.h b/tests/test6/compute.h
new file mode 100644
index 0000000..d663ba5
--- /dev/null
+++ b/tests/test6/compute.h
@@ -0,0 +1,11 @@
+#ifndef COMPUTE_H
+#define COMPUTE_H
+#include <stdlib.h> // for size_t
+#include <stdio.h> // for printf
+
+float compute_level(const float *data, size_t nsamples, int rate);
+void audio2hsv_1(int audio_level, int *light_h, int *light_s, int *light_v);
+void hsv2rgb(int h, int s, int v, int *r, int *g, int *b);
+
+#endif
+
diff --git a/tests/test6/fft.c b/tests/test6/fft.c
new file mode 100644
index 0000000..dd91d14
--- /dev/null
+++ b/tests/test6/fft.c
@@ -0,0 +1,241 @@
+#include "fft.h"
+#include <math.h>
+#include <stdlib.h>
+
+/* Everything here comes from Audacity 1.3.13
+ (orignally in C++ and with more genericity and functionnality)
+ Original Author : Dominic Mazzoni
+ Licenced under GPL 2.0 (see LICENCE)
+*/
+
+#define MaxFastBits 16
+
+int *gFFTBitTable[MaxFastBits]={NULL};
+
+void InitFFT();
+int NumberOfBitsNeeded(int PowerOfTwo);
+inline int FastReverseBits(int i, int NumBits);
+int ReverseBits(int index, int NumBits);
+
+
+void FFT(int NumSamples,
+		int InverseTransform,
+		float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
+{
+	int NumBits;                 /* Number of bits needed to store indices */
+	int i, j, k, n;
+	int BlockSize, BlockEnd;
+
+	double angle_numerator = 2.0 * M_PI;
+	double tr, ti;                /* temp real, temp imaginary */
+/*
+	if (!IsPowerOfTwo(NumSamples)) {
+		fprintf(stderr, "%d is not a power of two\n", NumSamples);
+		exit(1);
+	}
+*/
+	if (!gFFTBitTable[0])
+		InitFFT();
+
+	if (!InverseTransform)
+		angle_numerator = -angle_numerator;
+
+	NumBits = NumberOfBitsNeeded(NumSamples);
+
+	/*
+	 **   Do simultaneous data copy and bit-reversal ordering into outputs...
+	 */
+	for (i = 0; i < NumSamples; i++) {
+		j = FastReverseBits(i, NumBits);
+		RealOut[j] = RealIn[i];
+		ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
+	}
+
+	/*
+	 **   Do the FFT itself...
+	 */
+
+	BlockEnd = 1;
+	for (BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {
+
+		double delta_angle = angle_numerator / (double) BlockSize;
+
+		double sm2 = sin(-2 * delta_angle);
+		double sm1 = sin(-delta_angle);
+		double cm2 = cos(-2 * delta_angle);
+		double cm1 = cos(-delta_angle);
+		double w = 2 * cm1;
+		double ar0, ar1, ar2, ai0, ai1, ai2;
+
+		for (i = 0; i < NumSamples; i += BlockSize) {
+			ar2 = cm2;
+			ar1 = cm1;
+
+			ai2 = sm2;
+			ai1 = sm1;
+
+			for (j = i, n = 0; n < BlockEnd; j++, n++) {
+				ar0 = w * ar1 - ar2;
+				ar2 = ar1;
+				ar1 = ar0;
+
+				ai0 = w * ai1 - ai2;
+				ai2 = ai1;
+				ai1 = ai0;
+
+				k = j + BlockEnd;
+				tr = ar0 * RealOut[k] - ai0 * ImagOut[k];
+				ti = ar0 * ImagOut[k] + ai0 * RealOut[k];
+
+				RealOut[k] = RealOut[j] - tr;
+				ImagOut[k] = ImagOut[j] - ti;
+
+				RealOut[j] += tr;
+				ImagOut[j] += ti;
+			}
+		}
+		BlockEnd = BlockSize;
+	}
+
+	/*
+	 **   Need to normalize if inverse transform...
+	 */
+
+	if (InverseTransform) {
+		float denom = (float) NumSamples;
+
+		for (i = 0; i < NumSamples; i++) {
+			RealOut[i] /= denom;
+			ImagOut[i] /= denom;
+		}
+	}
+}
+
+void InitFFT() {
+	int len=2;
+	int b, i;
+	for (b=0; b<MaxFastBits; b++) {
+		gFFTBitTable[b]=malloc(len*sizeof(int));
+
+		for (i=0; i<len; i++)
+			gFFTBitTable[b][i] = ReverseBits(i, b+1);
+
+		len <<= 1;
+	}
+}
+
+int NumberOfBitsNeeded(int PowerOfTwo)
+{
+	int i;
+
+/*
+	if (PowerOfTwo < 2) {
+		fprintf(stderr, "Error: FFT called with size %d\n", PowerOfTwo);
+		exit(1);
+	}
+*/
+	for (i = 0;; i++)
+		if (PowerOfTwo & (1 << i))
+			return i;
+}
+
+inline int FastReverseBits(int i, int NumBits)
+{
+	if (NumBits <= MaxFastBits)
+		return gFFTBitTable[NumBits - 1][i];
+	else
+		return ReverseBits(i, NumBits);
+}
+
+int ReverseBits(int index, int NumBits)
+{
+   int i, rev;
+
+   for (i = rev = 0; i < NumBits; i++) {
+      rev = (rev << 1) | (index & 1);
+      index >>= 1;
+   }
+
+   return rev;
+}
+
+
+/*
+ * PowerSpectrum
+ *
+ * This function computes the same as RealFFT, above, but
+ * adds the squares of the real and imaginary part of each
+ * coefficient, extracting the power and throwing away the
+ * phase.
+ *
+ * For speed, it does not call RealFFT, but duplicates some
+ * of its code.
+ */
+
+void PowerSpectrum(float In[PSNumS], float Out[PSHalf]) {
+	int i;
+
+	float theta = M_PI / PSHalf;
+
+	float tmpReal[PSHalf];
+	float tmpImag[PSHalf];
+	float RealOut[PSHalf];
+	float ImagOut[PSHalf];
+
+	for (i=0; i<PSHalf; i++) {
+		tmpReal[i] = In[2*i];
+		tmpImag[i] = In[2*i+1];
+	}
+
+	FFT(PSHalf, 0, tmpReal, tmpImag, RealOut, ImagOut);
+
+	float wtemp = sin(0.5 * theta);
+
+	float wpr = -2.0 * wtemp * wtemp;
+	float wpi = -1.0 * sin(theta);
+	float wr = 1.0 + wpr;
+	float wi = wpi;
+
+	int i3;
+
+	float h1r, h1i, h2r, h2i, rt, it;
+	for (i=1; i < PSHalf/2; i++) {
+
+		i3 = PSHalf - i;
+
+		h1r = 0.5 * (RealOut[i] + RealOut[i3]);
+		h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
+		h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
+		h2i = -0.5 * (RealOut[i] - RealOut[i3]);
+
+		rt = h1r + wr * h2r - wi * h2i;
+		it = h1i + wr * h2i + wi * h2r;
+
+		Out[i] = rt * rt + it * it;
+
+		rt = h1r - wr * h2r + wi * h2i;
+		it = -h1i + wr * h2i + wi * h2r;
+
+		Out[i3] = rt * rt + it * it;
+
+		wr = (wtemp = wr) * wpr - wi * wpi + wr;
+		wi = wi * wpr + wtemp * wpi + wi;
+	}
+
+	rt = (h1r = RealOut[0]) + ImagOut[0];
+	it = h1r - ImagOut[0];
+	Out[0] = rt * rt + it * it;
+	rt = RealOut[PSHalf / 2];
+	it = ImagOut[PSHalf / 2];
+	Out[PSHalf / 2] = rt * rt + it * it;
+}
+
+void DeinitFFT() {
+	int i;
+	if (gFFTBitTable[0]) {
+		for (i=0;i<MaxFastBits;i++) {
+			free(gFFTBitTable[i]);
+		}
+	}
+}
+
diff --git a/tests/test6/fft.h b/tests/test6/fft.h
new file mode 100644
index 0000000..6066d52
--- /dev/null
+++ b/tests/test6/fft.h
@@ -0,0 +1,19 @@
+#ifndef FFT_H
+#define FFT_H
+
+/* Everything here comes from Audacity 1.3.13
+ (orignally in C++ and with more genericity and functionnality)
+ Original Author : Dominic Mazzoni
+ Licenced under GPL 2.0 (see LICENCE)
+*/
+
+#define PSNumS 256
+#define PSHalf 128
+void PowerSpectrum(float In[PSNumS], float Out[PSHalf]);
+
+void FFT(int NumSamples, int InverseTransform,
+	float *RealIn, float *ImagIn, float *RealOut, float *ImagOut);
+void DeinitFFT();
+
+#endif
+
diff --git a/tests/test6/test.raw b/tests/test6/test.raw
new file mode 120000
index 0000000..a266d33
--- /dev/null
+++ b/tests/test6/test.raw
@@ -0,0 +1 @@
+../../../../musique-test/V clean pizz.raw
\ No newline at end of file
diff --git a/tests/test6/test6.c b/tests/test6/test6.c
index abaa3e4..e5a2a12 100644
--- a/tests/test6/test6.c
+++ b/tests/test6/test6.c
@@ -2,6 +2,8 @@
 #include <stdio.h>
 #include <stdlib.h>
 
+#include "compute.h"
+
 typedef void (*cb_processdata_t)(int n, float *);
 
 
@@ -60,10 +62,12 @@ void parse_testfile(cb_processdata_t cb) {
 	FILE *fh=fopen("./test.raw", "r");
 	if (fh==NULL) return;
 
-	n=128;
+	//n=128;
+	n=512;
 	while ( (n=fread(f, sizeof(float), n, fh)) > 0 ) {
 		cb(n,f);
-		n=128+256*(rand()%7);
+		//n=128+256*(rand()%7);
+		n=512;
 	}
 	fclose(fh);
 }
@@ -81,10 +85,19 @@ void process_mean_max(int n, float *f) {
 	printf("%+.3f %+.3f %4i\n", mean, max, n);
 }
 
+void process_level(int n, float *f) {
+	int rate=24000;
+	float level;
+
+	level=compute_level(f, n, rate);
+	printf("%+.3f %4i\n", level, n);
+}
+
 int main() {
 	//test_todb_a();
 	//dump_testfile();
-	parse_testfile(process_mean_max);
+	//parse_testfile(process_mean_max);
+	parse_testfile(process_level);
 
 	return 0;
 }
-- 
cgit v1.2.3