Mercurial > hg > libxtract

changeset 113:72a9a393d5bd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
- Fixed bugs in xtract_flatness(), or at least added necessary documentation and error checking to avoid problems - Added xtract_is_denormal() helper function and XTRACT_DENORMAL_FOUND return code - Replaced all instances of log, sqrt, exp etc. with respective floating point counterparts (logf etc.) - Added check for architecture endianness to configure script - Bug fix to PD example, now no longer crashes if no arguments are given - Minor documentation updates
author Jamie Bullock <jamie@postlude.co.uk>
date Fri, 15 Feb 2008 12:43:13 +0000
parents a76501dc5307
children f5040ed4e555
files ChangeLog TODO configure.in examples/puredata/xtract~.c src/descriptors.c src/helper.c src/init.c src/libxtract.c src/scalar.c src/vector.c xtract/libxtract.h xtract/xtract_helper.h xtract/xtract_scalar.h xtract/xtract_vector.h
diffstat 14 files changed, 279 insertions(+), 103 deletions(-) [+]
line wrap: on
line diff
--- a/ChangeLog	Wed Jan 02 04:10:21 2008 +0000
+++ b/ChangeLog	Fri Feb 15 12:43:13 2008 +0000
@@ -1,3 +1,16 @@
+2008-2-15 Jamie Bullock <jamie@postlude.co.uk>
+    * version 0.5.9
+	* Fixed bugs in xtract_flatness(), or at least added necessary
+	documentation and error checking to avoid problems
+	* Added xtract_is_denormal() helper function and XTRACT_DENORMAL_FOUND
+	return code
+	* Replaced all instances of log, sqrt, exp etc. with respective
+	floating point counterparts (logf etc.)
+	* Added check for architecture endianness to configure script
+	* Bug fix to PD example, now no longer crashes if no arguments are
+	given
+	* Minor documentation updates
+
 2007-10-09 Dan Stowell
     * version 0.5.1
         * Fixed bug in xtract_init_mfcc() which accidentally missed filling in the very top filter frequency value
--- a/TODO	Wed Jan 02 04:10:21 2008 +0000
+++ b/TODO	Fri Feb 15 12:43:13 2008 +0000
@@ -15,3 +15,6 @@
 Use the fftw guru interface to create multipurpose global plans from xtract_fft_init()
 Add LPC via Schur algorithm
 If argv uses an integer, then it should point to a variable of type int. At the moment floats are used to store ints.
+Write macro to truncate numbers <= 0 prior to taking log and integrate into functions across the library
+Add xtract_chords() - gets pitch classes from spectrum either as MIDI note numbers or pitch class and octave class - e.g. 'F#4' etc.
+Add xtract_subbands() - gets subband power or subband magnitude for given bands
--- a/configure.in	Wed Jan 02 04:10:21 2008 +0000
+++ b/configure.in	Fri Feb 15 12:43:13 2008 +0000
@@ -4,7 +4,7 @@
 # Increment for feature additions and enhancements
 m4_define(libxtract_minor_version, 5)
 # Increment for fixes 
-m4_define(libxtract_fix_version, 0)
+m4_define(libxtract_fix_version, 9)
 
 m4_define(libxtract_version, libxtract_major_version.libxtract_minor_version.libxtract_fix_version)
 		
@@ -17,6 +17,7 @@
 AC_PROG_CC
 AC_PROG_LIBTOOL
 AC_PROG_INSTALL
+AC_C_BIGENDIAN
 AC_PATH_PROG(PKG_CONFIG, pkg-config, no)
 AC_ENABLE_STATIC(no)
 AC_ENABLE_SHARED(yes)
@@ -125,7 +126,6 @@
 
 AM_CONDITIONAL(BUILD_SIMPLETEST, test "x${simpletest}" = 'xtrue')
 
-
 dnl Are we building the PD examples?
 if [[ "$pd_example" = "true" ]] ; then
     PD_SOURCES="xtract~.c"
@@ -217,6 +217,14 @@
 
 AM_CONDITIONAL(BUILD_FFT, test "x${fft}" = 'xtrue')
 
+dnl Check for architecture endian-ness
+#AC_C_BIGENDIAN(bigendian=true, bigendian=false, bigendian=undefined)
+#if [[ "$is_bigendian" = "false" ]] ; then
+#    AC_DEFINE([WORDS_BIGENDIAN], [0], [Architecture is big endian])
+#else
+#    AC_DEFINE([WORDS_BIGENDIAN], [1], [Architecture is not big endian])
+#fi
+
 
 dnl ------------------------------------------
 dnl ----  do some magic to gues the host opsys
--- a/examples/puredata/xtract~.c	Wed Jan 02 04:10:21 2008 +0000
+++ b/examples/puredata/xtract~.c	Fri Feb 15 12:43:13 2008 +0000
@@ -168,6 +168,10 @@
         else
             x->feature_name = atom_getsymbol(argv);
     }
+    else {
+	post("xtract~: No arguments given");
+        return (void *)x;
+    }
     if(argc > 1){
         if(x->is_subframe)
             x->feature_name = atom_getsymbol(argv+1);
@@ -233,8 +237,6 @@
 	if(strcmp(author, "") && year)	
 	    post("xtract~: %s(%d)", author, year);
     }	
-    else
-	post("xtract~: No arguments given");
     
     /* Adjust frame size if we are using subframe features */
     if(x->is_subframe)
--- a/src/descriptors.c	Wed Jan 02 04:10:21 2008 +0000
+++ b/src/descriptors.c	Fri Feb 15 12:43:13 2008 +0000
@@ -530,7 +530,7 @@
                 break;
             case XTRACT_ODD_EVEN_RATIO:
                 strcpy(name, "odd_even_ratio");
-                strcpy(p_name, "Odd/Even Harmonic Ratio");
+                strcpy(p_name, "Odd/even Harmonic Ratio");
                 strcpy(desc, 
                         "Extract the odd-to-even harmonic ratio of a spectrum");
                 strcpy(p_desc, 
--- a/src/helper.c	Wed Jan 02 04:10:21 2008 +0000
+++ b/src/helper.c	Fri Feb 15 12:43:13 2008 +0000
@@ -21,8 +21,18 @@
 
 /* helper.c: helper functions. */
 
+#include <config.h>
+
+#include <stdio.h>
+
 #include "xtract/libxtract.h"
 
+#ifdef WORDS_BIGENDIAN
+#define INDEX 0
+#else
+#define INDEX 1
+#endif
+
 int xtract_windowed(const float *data, const int N, const void *argv, float *result){
 
     int n;
@@ -45,8 +55,7 @@
     float *result1,
           *result2;
 
-    int n, i,
-        rv;
+    int n, rv;
 
     n = N >> 1;
 
@@ -63,3 +72,12 @@
     return rv;
 
 }
+
+inline int xtract_is_denormal(double const d){
+    if(sizeof(d) != 2 * sizeof(int))
+        fprintf(stderr, "libxtract: Error: xtract_is_denormal() detects inconsistent wordlength for type 'double'\n");
+
+    int l = ((int *)&d)[INDEX];
+    return (l&0x7ff00000) == 0 && d!=0; //Check for 0 may not be necessary
+}
+
--- a/src/init.c	Wed Jan 02 04:10:21 2008 +0000
+++ b/src/init.c	Fri Feb 15 12:43:13 2008 +0000
@@ -149,7 +149,7 @@
 
     input = output = NULL;
 
-    fprintf(stderr, "Optimisation level: %d\n", XTRACT_FFT_OPTIMISATION_LEVEL);
+    //fprintf(stderr, "Optimisation level: %d\n", XTRACT_FFT_OPTIMISATION_LEVEL);
 
     if(XTRACT_FFT_OPTIMISATION_LEVEL == 0)
         optimisation = FFTW_ESTIMATE;
--- a/src/libxtract.c	Wed Jan 02 04:10:21 2008 +0000
+++ b/src/libxtract.c	Fri Feb 15 12:43:13 2008 +0000
@@ -48,6 +48,7 @@
     xtract_rolloff,
     xtract_loudness,
     xtract_flatness,
+    xtract_flatness_db,
     xtract_tonality,
     xtract_crest,
     xtract_noisiness,
--- a/src/scalar.c	Wed Jan 02 04:10:21 2008 +0000
+++ b/src/scalar.c	Fri Feb 15 12:43:13 2008 +0000
@@ -21,12 +21,16 @@
 
 /* scalar.c: defines functions that extract a feature as a single value from an input vector */
 
-#include "xtract/libxtract.h"
-#include "xtract_macros_private.h"
-#include "math.h"
+#include <config.h>
+
 #include <stdlib.h>
 #include <string.h>
 #include <stdio.h>
+#include <math.h>
+
+#include "xtract/libxtract.h"
+#include "xtract/xtract_helper.h"
+#include "xtract_macros_private.h"
 
 #ifndef powf
     #define powf pow
@@ -36,8 +40,20 @@
     #define expf exp
 #endif
 
+#ifndef sqrtf
+    #define sqrtf sqrt
+#endif
+
+#ifndef fabsf
+    #define fabsf fabs
+#endif
+
+
 void test(void){
-    printf("Hello world");
+    printf("Hello world\n");
+#ifdef WORDS_BIGENDIAN
+    printf("Big endian!\n");
+#endif
 }
 
 int xtract_mean(const float *data, const int N, const void *argv, float *result){
@@ -58,8 +74,10 @@
 
     int n = N;
 
+    *result += 0.f;
+
     while(n--)
-	*result += pow(data[n] - *(float *)argv, 2);
+	*result += powf(data[n] - *(float *)argv, 2);
 
     *result = *result / (N - 1);
 
@@ -68,7 +86,7 @@
 
 int xtract_standard_deviation(const float *data, const int N, const void *argv, float *result){
 
-    *result = sqrt(*(float *)argv);
+    *result = sqrtf(*(float *)argv);
 
     return XTRACT_SUCCESS;
 }
@@ -77,8 +95,10 @@
 
     int n = N;
 
+    *result += 0.f;
+
     while(n--)
-	*result += fabs(data[n] - *(float *)argv);
+	*result += fabsf(data[n] - *(float *)argv);
 
     *result /= N;
 
@@ -91,9 +111,11 @@
 
     float temp = 0.f;
 
+    *result += 0.f;
+
     while(n--){
 	temp = (data[n] - ((float *)argv)[0]) / ((float *)argv)[1];
-	*result += pow(temp, 3);
+	*result += powf(temp, 3);
     }
 
     *result /= N;
@@ -106,11 +128,13 @@
 
     int n = N;
 
-    float temp;
+    float temp = 0.f;
+
+    *result += 0.f;
 
     while(n--){
 	temp = (data[n] - ((float *)argv)[0]) / ((float *)argv)[1];
-	*result += pow(temp, 4);
+	*result += powf(temp, 4);
     }
 
     *result /= N;
@@ -134,7 +158,10 @@
 	A += amps[n];
     }
 
-    *result = FA / A;
+    if(A == 0.f)
+        *result = 0.f;
+    else
+        *result = FA / A;
 
     return XTRACT_SUCCESS;
 }
@@ -156,6 +183,8 @@
     amps = data;
     freqs = data + m;
 
+    *result += 0.f;
+
     while(m--){
 	A += amps[m];
 	*result += powf((freqs[m]  - *(float *)argv) * amps[m], 2);
@@ -168,7 +197,7 @@
 
 int xtract_spectral_standard_deviation(const float *data, const int N, const void *argv, float *result){
 
-    *result = sqrt(*(float *)argv); 
+    *result = sqrtf(*(float *)argv); 
 
     return XTRACT_SUCCESS;
 }
@@ -184,9 +213,11 @@
     amps = data;
     freqs = data + m;
 
+    *result += 0.f;
+
     while(m--){
 	A += amps[m];
-	*result += fabs((amps[m] * freqs[m]) - *(float *)argv);
+	*result += fabsf((amps[m] * freqs[m]) - *(float *)argv);
     }
 
     *result /= A;
@@ -205,11 +236,13 @@
     amps = data;
     freqs = data + m;
 
+    *result += 0.f;
+
     while(m--){
 	A += amps[m];
 	temp = ((amps[m] * freqs[m]) - 
 		((float *)argv)[0]) / ((float *)argv)[1];
-	*result += pow(temp, 3);
+	*result += powf(temp, 3);
     }
 
     *result /= A;
@@ -228,11 +261,13 @@
     amps = data;
     freqs = data + m;
 
+    *result += 0.f;
+
     while(m--){
 	A += amps[m];
 	temp = ((amps[m] * freqs[m]) - 
 		((float *)argv)[0]) / ((float *)argv)[1];
-	*result += pow(temp, 4);
+	*result += powf(temp, 4);
     }
 
     *result /= A;
@@ -246,10 +281,10 @@
     int n,
 	M = N - 1;
 	
-	*result = 0.f;
+    *result = 0.f;
 	
     for(n = 1; n < M; n++)
-	*result += fabs(data[n] - (data[n-1] + data[n] + data[n+1]) / 3);
+	*result += fabsf(data[n] - (data[n-1] + data[n] + data[n+1]) / 3.f);
 
     return XTRACT_SUCCESS;
 }
@@ -261,8 +296,8 @@
     double num = 0.f, den = 0.f;
 
     while(n--){
-	num += pow(data[n] - data[n+1], 2);
-	den += pow(data[n], 2);
+	num += powf(data[n] - data[n+1], 2);
+	den += powf(data[n], 2);
     }
 
     *result = (float)(num / den);
@@ -286,18 +321,23 @@
 	}
     }
 
-    *result = p1 / den;
-
-    return XTRACT_SUCCESS;
+    if(den == 0.f || p1 == 0.f){
+        *result = 0.f;
+        return XTRACT_NO_RESULT;
+    }
+    else{
+        *result = p1 / den;
+        return XTRACT_SUCCESS;
+    }
 }
 
 int xtract_tristimulus_2(const float *data, const int N, const void *argv, float *result){
 
     int n = N;
 
-    float den, p2, p3, p4, temp;
+    float den, p2, p3, p4, ps, temp;
 
-    den = p2 = p3 = p4 = temp = 0.f;
+    den = p2 = p3 = p4 = ps = temp = 0.f;
 
     for(n = 0; n < N; n++){
 	if((temp = data[n])){
@@ -311,9 +351,17 @@
 	}
     }
 
-    *result = (p2 + p3 + p4)  / den;
+    ps = p2 + p3 + p4;
 
-    return XTRACT_SUCCESS;
+    if(den == 0.f || ps == 0.f){
+        *result = 0.f;
+        return XTRACT_NO_RESULT;
+    }
+    else{
+        *result = ps / den;
+        return XTRACT_SUCCESS;
+    }
+
 }
 
 int xtract_tristimulus_3(const float *data, const int N, const void *argv, float *result){
@@ -333,9 +381,14 @@
 	}
     }
 
-    *result = num / den;
-
-    return XTRACT_SUCCESS;
+    if(den == 0.f || num == 0.f){
+        *result = 0.f;
+        return XTRACT_NO_RESULT;
+    }
+    else{
+        *result = num / den;
+        return XTRACT_SUCCESS;
+    }
 }
 
 int xtract_smoothness(const float *data, const int N, const void *argv, float *result){
@@ -347,14 +400,18 @@
     input = (float *)malloc(N * sizeof(float));
     memcpy(input, data, N * sizeof(float));
 
-    if (input[0] <= 0) input[0] = 1;
+    if (input[0] <= 0) 
+        input[0] = XTRACT_LOG_LIMIT;
+    if (input[1] <= 0) 
+        input[1] = XTRACT_LOG_LIMIT;
 
     M = N - 1;
 
     for(n = 1; n < M; n++){ 
-	if(input[n] <= 0) input[n] = 1;
-	*result += fabs(20 * log(input[n]) - (20 * log(input[n-1]) + 
-		    20 * log(input[n]) + 20 * log(input[n+1])) / 3);
+	if(input[n+1] <= 0) 
+            input[n+1] = XTRACT_LOG_LIMIT;
+	*result += fabsf(20.f * logf(input[n]) - (20.f * logf(input[n-1]) + 
+		    20.f * logf(input[n]) + 20.f * logf(input[n+1])) / 3.f);
     }
 
     free(input);
@@ -377,7 +434,7 @@
 	den += data[n];
     }
 
-    *result = sqrt(num / den);
+    *result = sqrtf(num / den);
 
     return XTRACT_SUCCESS;
 }
@@ -389,7 +446,7 @@
     for(n = 1; n < N; n++)
 	if(data[n] * data[n-1] < 0) (*result)++;
 
-    *result /= N;
+    *result /= (float)N;
 
     return XTRACT_SUCCESS;
 }
@@ -419,6 +476,8 @@
 
     int n = N, rv;
 
+    *result = 0.f;
+
     if(n > XTRACT_BARK_BANDS){
 	n = XTRACT_BARK_BANDS; 
 	rv = XTRACT_BAD_VECTOR_SIZE; 
@@ -434,30 +493,56 @@
 
 int xtract_flatness(const float *data, const int N, const void *argv, float *result){
 
-    int n;
+    int n, count, denormal_found;
 
     double num, den, temp; 
 
-    den = data[0];
-    num = (data[0] == 0.f ? 1.f : data[0]);
+    num = 1.f;
+    den = temp = 0.f;
 
-    for(n = 1; n < N; n++){
-	if((temp = data[n]) != 0.f) {
-	    num *= temp;
-	    den += temp;
-	}
+    denormal_found = 0;
+    count = 0;
+
+    for(n = 0; n < N; n++){
+        if((temp = data[n]) != 0.f) {
+            if (xtract_is_denormal(num)){
+                denormal_found = 1;
+                break;
+            }
+            num *= temp;
+            den += temp;
+            count++;
+        }
     }
 
-    num = powf(num, 1.f / N);
-    den /= N;
+    if(!count){
+        *result = 0.f;
+        return XTRACT_NO_RESULT;
+    }
 
-    if(num < XTRACT_VERY_SMALL_NUMBER) 
-	num = XTRACT_VERY_SMALL_NUMBER;
+    num = powf(num, 1.f / (float)N);
+    den /= (float)N;
 
-    if(den < XTRACT_VERY_SMALL_NUMBER) 
-	den = XTRACT_VERY_SMALL_NUMBER;
 
-    *result = 10 * log10(num / den);
+    *result = (float) (num / den);
+
+    if(denormal_found)
+        return XTRACT_DENORMAL_FOUND;
+    else
+        return XTRACT_SUCCESS;
+    
+}
+
+int xtract_flatness_db(const float *data, const int N, const void *argv, float *result){
+
+    float flatness_db;
+
+    flatness_db = *(float *)argv;
+
+    if (flatness_db <= 0) 
+        flatness_db = XTRACT_LOG_LIMIT;
+
+    *result = 10 * log10f(flatness_db);
 
     return XTRACT_SUCCESS;
 
@@ -465,13 +550,11 @@
 
 int xtract_tonality(const float *data, const int N, const void *argv, float *result){
 
-    float sfmdb, sfm;
+    float sfmdb;
 
-    sfm = *(float *)argv;
+    sfmdb = *(float *)argv;
 
-    sfmdb = sfm / -60.f;
-
-    *result = XTRACT_MIN(sfmdb, 1);
+    *result = XTRACT_MIN(sfmdb / -60.f, 1);
 
     return XTRACT_SUCCESS;
 }
@@ -512,9 +595,11 @@
 
     int n = N;
 
+    *result = 0.f;
+
     while(n--) *result += XTRACT_SQ(data[n]);
 
-    *result = sqrt(*result / N);
+    *result = sqrtf(*result / (float)N);
 
     return XTRACT_SUCCESS;
 }
@@ -531,7 +616,7 @@
 
     while(n--){
 	if(amps[n]){
-	    num += fabs(freqs[n] - n * fund) * XTRACT_SQ(amps[n]);
+	    num += fabsf(freqs[n] - n * fund) * XTRACT_SQ(amps[n]);
 	    den += XTRACT_SQ(amps[n]);
 	}
     }
@@ -552,30 +637,37 @@
 
     int M = (N >> 1), n;
 
-    float num = 0.f, den = 0.f,  temp;
+    float odd = 0.f, even = 0.f,  temp;
 
     for(n = 0; n < M; n++){
 	if((temp = data[n])){
 	    if(XTRACT_IS_ODD(n)){
-		num += temp;
+		odd += temp;
 	    }
 	    else{
-		den += temp;
+		even += temp;
 	    }
 	}
     }
 
-    *result = num / den;
-
-    return XTRACT_SUCCESS;
+    if(odd == 0.f || even == 0.f){
+        *result = 0.f;
+        return XTRACT_NO_RESULT;
+    }
+    else {
+        *result = odd / even;
+        return XTRACT_SUCCESS;
+    }
 }
 
 int xtract_sharpness(const float *data, const int N, const void *argv, float *result){
 
     int n = N, rv;
-    float sl, g, temp; /* sl = specific loudness */
+    float sl, g; /* sl = specific loudness */
+    double temp; 
 
-    sl = g = temp = 0.f;
+    sl = g = 0.f;
+    temp = 0.f;
 
     if(n > XTRACT_BARK_BANDS) 
 	rv = XTRACT_BAD_VECTOR_SIZE; 
@@ -589,7 +681,8 @@
 	temp += n * g * sl;
     }
 
-    *result = 0.11 * temp / N;
+    temp = 0.11 * temp / (float)N;
+    *result = (float)temp;
 
     return rv;
 
@@ -655,6 +748,8 @@
 
     int n = N;
 
+    *result = 0.f;
+
     while(n--)
 	*result += *data++;
 
@@ -665,6 +760,8 @@
 int xtract_nonzero_count(const float *data, const int N, const void *argv, float *result){
 
     int n = N;
+    
+    *result += 0.f;
 
     while(n--)
 	*result += (*data++ ? 1 : 0);
@@ -690,7 +787,7 @@
     while(n--) coeffs2[n] = coeffs3[n] = 1;
 
     M = N >> 1;
-    L = N / 3;
+    L = N / 3.f;
 
     while(M--){
 	m = M << 1;
@@ -698,7 +795,7 @@
 
 	if(M < L){
 	    l = M * 3;
-	    coeffs3[M] = (data[l] + data[l+1] + data[l+2]) / 3;
+	    coeffs3[M] = (data[l] + data[l+1] + data[l+2]) / 3.f;
 	}
     }
 
@@ -792,12 +889,12 @@
 
     /* Estimate fundamental freq */
     for (n = 1; n < M; n++)
-	err_tau_1 = err_tau_1 + fabs(input[n] - input[n+1]);
+	err_tau_1 = err_tau_1 + fabsf(input[n] - input[n+1]);
     /* FIX: this doesn't pose too much load if it returns 'early', but if it can't find f0, load can be significant for larger block sizes M^2 iterations! */  
     for (tau = 2; tau < M; tau++){
 	err_tau_x = 0;
 	for (n = 1; n < M; n++){
-	    err_tau_x = err_tau_x + fabs(input[n] - input[n+tau]);
+	    err_tau_x = err_tau_x + fabsf(input[n] - input[n+tau]);
 	}
 	if (err_tau_x < err_tau_1) {
 	    f0 = sr / (tau + (err_tau_x / err_tau_1));
--- a/src/vector.c	Wed Jan 02 04:10:21 2008 +0000
+++ b/src/vector.c	Fri Feb 15 12:43:13 2008 +0000
@@ -37,6 +37,22 @@
     }
 #endif
 
+#ifndef powf
+    #define powf pow
+#endif
+
+#ifndef expf
+    #define expf exp
+#endif
+
+#ifndef sqrtf
+    #define sqrtf sqrt
+#endif
+
+#ifndef fabsf
+    #define fabsf fabs
+#endif
+
 #ifdef XTRACT_FFT
 
 #include <fftw3.h>
@@ -45,11 +61,10 @@
 
 int xtract_spectrum(const float *data, const int N, const void *argv, float *result){
 
-    float *input, *rfft, q, temp, max;
+    float *input, *rfft, q, temp, max, NxN;
     size_t bytes;
     int n,
         m,
-        NxN, 
         M, 
         vector, 
         withDC, 
@@ -89,7 +104,7 @@
 	    for(n = 1; n < M; n++){
 		if ((temp = XTRACT_SQ(rfft[n]) + 
 			    XTRACT_SQ(rfft[N - n])) > XTRACT_LOG_LIMIT)
-		    temp = log(sqrt(temp) / N);
+		    temp = logf(sqrtf(temp) / (float)N);
 		else
 		    temp = XTRACT_LOG_LIMIT_DB;
 
@@ -130,7 +145,7 @@
 	    for(n = 1; n < M; n++){
 		if ((temp = XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) > 
 			XTRACT_LOG_LIMIT)
-		    temp = log(temp / NxN);
+		    temp = logf(temp / NxN);
 		else
 		    temp = XTRACT_LOG_LIMIT_DB; 		
 
@@ -161,8 +176,8 @@
 		    result[M + m] = n * q;
                 }
 
-                result[m] = sqrt(XTRACT_SQ(rfft[n]) + 
-                        XTRACT_SQ(rfft[N - n])) / N; 
+                result[m] = sqrtf(XTRACT_SQ(rfft[n]) + 
+                        XTRACT_SQ(rfft[N - n])) / (float)N; 
                 max = result[m] > max ? result[m] : max;
 	    }
 	    break;
@@ -254,7 +269,7 @@
         for(n = 0; n < N; n++){
             result[filter] += data[n] * f->filters[filter][n];
         }
-        result[filter] = log(result[filter] < XTRACT_LOG_LIMIT ? XTRACT_LOG_LIMIT : result[filter]);
+        result[filter] = logf(result[filter] < XTRACT_LOG_LIMIT ? XTRACT_LOG_LIMIT : result[filter]);
     }
 
     xtract_dct(result, f->n_filters, NULL, result);
@@ -334,13 +349,13 @@
     float md, temp;
 
     while(n--){
-       md = 0;
+       md = 0.f;
         for(i = 0; i < N - n; i++){
             temp = data[i] - data[i + n];
 			temp = (temp < 0 ? -temp : temp);
 			md += temp;
         }
-        result[n] = md / N;
+        result[n] = md / (float)N;
     }
 
     return XTRACT_SUCCESS;
@@ -353,12 +368,12 @@
     float sd;
 
     while(n--){
-       sd = 0;
+       sd = 0.f;
         for(i = 0; i < N - n; i++){
             /*sd = 1;*/
             sd += XTRACT_SQ(data[i] - data[i + n]);
         }
-        result[n] = sd / N;
+        result[n] = sd / (float)N;
     }
 
     return XTRACT_SUCCESS;
@@ -488,13 +503,12 @@
     M = L * 2; /* The length of *result */
     ref = result;
     lpc = result+L;
-/*
+
     if(error == 0.0){
-        for(i = 0; i < M; i++)
-            result[i] = 0.f;
+        memset(result, 0, M * sizeof(float));
         return XTRACT_NO_RESULT;
     }
-*/
+
     memset(result, 0, M * sizeof(float));
 
     for (i = 0; i < L; i++) {
--- a/xtract/libxtract.h	Wed Jan 02 04:10:21 2008 +0000
+++ b/xtract/libxtract.h	Fri Feb 15 12:43:13 2008 +0000
@@ -96,6 +96,7 @@
     XTRACT_ROLLOFF,
     XTRACT_LOUDNESS,
     XTRACT_FLATNESS,
+    XTRACT_FLATNESS_DB,
     XTRACT_TONALITY,
     XTRACT_CREST,
     XTRACT_NOISINESS,
@@ -165,7 +166,8 @@
     XTRACT_MALLOC_FAILED,
     XTRACT_BAD_ARGV,
     XTRACT_BAD_VECTOR_SIZE,
-    XTRACT_NO_RESULT,
+    XTRACT_DENORMAL_FOUND,
+    XTRACT_NO_RESULT, /* This usually occurs when the correct calculation cannot take place because required data is missing or would result in a NaN or infinity/-infinity. Under these curcumstances 0.f is usually given by *result */
     XTRACT_FEATURE_NOT_IMPLEMENTED
 };
 
--- a/xtract/xtract_helper.h	Wed Jan 02 04:10:21 2008 +0000
+++ b/xtract/xtract_helper.h	Fri Feb 15 12:43:13 2008 +0000
@@ -63,6 +63,9 @@
  */
 int xtract_features_from_subframes(const float *data, const int N, const int feature, const void *argv, float *result);
 
+/** \brief Test whether a number is denormal */
+inline int xtract_is_denormal(double const d);
+
 /** @} */
 
 #ifdef __cplusplus
--- a/xtract/xtract_scalar.h	Wed Jan 02 04:10:21 2008 +0000
+++ b/xtract/xtract_scalar.h	Fri Feb 15 12:43:13 2008 +0000
@@ -239,21 +239,34 @@
  */
 int xtract_loudness(const float *data, const int N, const void *argv, float *result);
 
-/** \brief Extract the spectral flatness measure of an input vector using a method described by Tristan Jehan (2005)
+/** \brief Extract the spectral flatness measure of an input vector, where the flatness measure (SFM) is defined as the ratio of the geometric mean to the arithmetic mean of a magnitude spectrum.
+ *
+ * \note The computation method used here is the most efficient by a significant margin, but suffers from precision problems due to the multiplication operationin the geometric mean calculation. This is particularly accute for larger values of N (>=256). However, as noted by Peeters (2003), the SFM should generally be computed on a small number of 'bands' rather than on the complete magnitude spectrum. It is therefore highly recommended that xtract_bands() is used prior to calling xtract_flatness().
  * 
- * \param *data: a pointer to the first element in an array of floats representing the magnitude coefficients from the spectrum of an audio vector, (e.g. the first half of the array pointed to by *result from xtract_spectrum().
- * \param N: the number of elements to be considered
+ * \param *data: a pointer to the first element in an array of floats representing the magnitude coefficients from the spectrum of an audio vector, (e.g. the first half of the array pointed to by *result from xtract_spectrum(). Alternatively the magnitudes from a number of 'subbands' can be used by using *result from xtract_bands().
+ * \param N: the number of *data array elements to be considered
  * \param *argv: a pointer to NULL
- * \param *result: the spectral flatness of N values from the array pointed to by *data
+ * \param *result: the flatness of N values from the array pointed to by *data
  */
 int xtract_flatness(const float *data, const int N, const void *argv, float *result);
 
+/** \brief Extract the LOG spectral flatness measure of an input vector
+ *
+ * \param *data: a pointer to NULL.
+ * \param N: not used - can safely be set to 0.
+ * \param *argv: a pointer to a float represnting spectral flatness.
+ * \param *result: the LOG spectral flatness of N values from the array pointed to by *data
+ *
+ *  flatness_db = 10 * log10(flatness)
+ *
+ */
+int xtract_flatness_db(const float *data, const int N, const void *argv, float *result);
 
-/** \brief Extract the tonality factor of an input vector using a method described by Tristan Jehan (2005)
+/** \brief Extract the tonality factor of an input vector using a method described by Peeters 2003
  * 
- * \param *data: not used.
- * \param N: not used
- * \param *argv: a pointer to the spectral flatness measure of an audio vector (e.g. the output from xtract_flatness)
+ * \param *data: a pointer to NULL.
+ * \param N: not used - can safely be set to 0.
+ * \param *argv: a pointer to the LOG spectral flatness measure of an audio vector (e.g. the output from xtract_flatness_db)
  * \param *result: the tonality factor of N values from the array pointed to by *data
  */
 int xtract_tonality(const float *data, const int N, const void *argv, float *result);
@@ -309,7 +322,7 @@
  * \param *data: a pointer to the first element in an array of floats representing the amplitudes of the harmonic spectrum of an audio vector. It is sufficient to pass in a pointer to the first half of the array pointed to by *result from xtract_harmonic_spectrum().
  * \param N: the number of elements to be considered. If using the array pointed to by *result from xtract_harmonics, N should equal half the total array size i.e., just the amplitudes of the peaks.
  * \param *argv: a pointer to NULL
- * \param *result: the odd/even harmonic ratio of N values from the array pointed to by *data
+ * \param *result: the even/odd harmonic ratio of N values from the array pointed to by *data
  */
 int xtract_odd_even_ratio(const float *data, const int N, const void *argv, float *result);
 
--- a/xtract/xtract_vector.h	Wed Jan 02 04:10:21 2008 +0000
+++ b/xtract/xtract_vector.h	Fri Feb 15 12:43:13 2008 +0000
@@ -43,6 +43,8 @@
  *
  * The magnitude/power coefficients are scaled to the range 0-1 so that for a given coefficient x, 0 <= x <= 1
  *
+ * \note Before calling xtract_spectrum(), the FFT must be initialised by calling xtract_init_fft(N, XTRACT_SPECTRUM)
+ *
  */
 int xtract_spectrum(const float *data, const int N, const void *argv, float *result);