diff options
author | Ludovic Pouzenc <ludovic@pouzenc.fr> | 2012-06-04 21:44:28 +0000 |
---|---|---|
committer | Ludovic Pouzenc <ludovic@pouzenc.fr> | 2012-06-04 21:44:28 +0000 |
commit | 5b385abef33b20c122bddf9cc9e947594e0ebd32 (patch) | |
tree | 2cda5c9fc70fa79ac880b4f744ffb7f65769ed33 /tests/test5/compute.c | |
parent | 9c7c14513dfe3440dad9cf0b1c652dde9e4f1de7 (diff) | |
download | 2012-violon-leds-5b385abef33b20c122bddf9cc9e947594e0ebd32.tar.gz 2012-violon-leds-5b385abef33b20c122bddf9cc9e947594e0ebd32.tar.bz2 2012-violon-leds-5b385abef33b20c122bddf9cc9e947594e0ebd32.zip |
Bon. Partie pulse audio finie je pense.
début de la partie galère sur le "vrai" calcul pour le vu-mètre.
C'est compliqué car si on veut du dbA il faut faire une FFT pour appliquer des poids par fréquence.
Analyse fréquentielle copiée depuis le projet Audacity (adaptée du C++ au C et décimée).
Il y a des tas de petits mallocs pour la FFT et ça pue. D'ailleurs l'exécution de cette version donne un assertion failed sur malloc() que j'avais jamais vu...
git-svn-id: file:///var/svn/2012-violon-leds/trunk@12 6be1fa4d-33ac-4c33-becc-79fcb3794bb6
Diffstat (limited to 'tests/test5/compute.c')
-rw-r--r-- | tests/test5/compute.c | 293 |
1 files changed, 293 insertions, 0 deletions
diff --git a/tests/test5/compute.c b/tests/test5/compute.c index 11a7f81..f6f08c4 100644 --- a/tests/test5/compute.c +++ b/tests/test5/compute.c @@ -1,5 +1,298 @@ #include "compute.h" +#include <stdlib.h> +#include <math.h> + +#define MaxFastBits 16 + +int **gFFTBitTable = NULL; + +gfloat compute_level(const float *data, size_t nsamples, size_t nchan) { + + size_t i; + float level=0; + float *input = malloc(nsamples*sizeof(float)); + float *output = malloc(nsamples*sizeof(float)); + + double rate=44100; //TODO dynamique +/* Just return the max peak + for (i=0;i<nsamples;i+=nchan) { + val=((float *)data)[i]; + //printf("val==%i\n", val); + if (val<0) val=-val; + if (level<val) level=val; + } +*/ + for (i=0;i<nsamples;i++) { + input[i]=data[i*nchan]; + } + + compute_spectrom(input, nsamples, rate, output); + + for (i=0;i<nsamples;i++) { + level+=output[i]; + } + + return level/nsamples; +} + +// From Audacity +void compute_spectrom(float * data, int width, double rate, float *output) { + + int i; + float processed[256]={ 0.0f }; + float in[256]; + float out[256]; + + int start = 0; + int windows = 0; + while (start + 256 <= width) { + for (i=0; i<256; i++) + in[i] = data[start + i]; + + // Hanning + for (i=0; i<256; i++) + in[i] *= 0.50 - 0.50 * cos(2 * M_PI * i / (256 - 1)); + break; + + PowerSpectrum(in, out); + + // Take real part of result + for (i=0; i<256/2; i++) + processed[i] += out[i]; + + start += 256/2; + windows++; + } + // Convert to decibels + // But do it safely; -Inf is nobody's friend + for (i = 0; i < 256/2; i++){ + float temp=(processed[i] / 256 / windows); + if (temp > 0.0) + processed[i] = 10*log10(temp); + else + processed[i] = 0; + } + for(i=0;i<256/2;i++) + output[i] = processed[i]; +} + +/* + * 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, float *Out) +{ + int i; + + float theta = M_PI / 128; + + float tmpReal[128]; + float tmpImag[128]; + float RealOut[128]; + float ImagOut[128]; + + for (i = 0; i < 128; i++) { + tmpReal[i] = In[2 * i]; + tmpImag[i] = In[2 * i + 1]; + } + + FFT(128, 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 < 128 / 2; i++) { + + i3 = 128 - 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[128 / 2]; + it = ImagOut[128 / 2]; + Out[128 / 2] = rt * rt + it * it; +} + +void FFT(int NumSamples, + gboolean 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) + 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() +{ + gFFTBitTable = malloc(MaxFastBits*sizeof(int)); + + int len = 2; + int b, i; + for (b=1; b<=MaxFastBits; b++) { + gFFTBitTable[b-1]=malloc(len*sizeof(int)); + + for (i=0; i<len; i++) + gFFTBitTable[b-1][i] = ReverseBits(i, b); + + 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; +} + + void audio2hsv_1(gint audio_level, gint *light_h, gint *light_s, gint *light_v) { // Dummy code *light_h=-audio_level; |