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
(limited to 'tests/test6')
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