Merge pull request #53 from seanlikeskites/tristimulusWork

Possible Corrections to Tristimulus Functions

2 files changed, 60 insertions, 43 deletions

diff --git a/src/scalar.c b/src/scalar.c
index e6c6687..74d6cee 100644
--- a/src/scalar.c
+++ b/src/scalar.c
@@ -300,20 +300,21 @@ int xtract_irregularity_j(const double *data, const int N, const void *argv, dou
 
 int xtract_tristimulus_1(const double *data, const int N, const void *argv, double *result)
 {
+    int n = N >> 1, h = 0, i;
+    double den = 0.0, p1 = 0.0, fund = 0.0, temp = 0.0;
+    const double *freqs;
 
-    int n = N;
-
-    double den, p1, temp;
-
-    den = p1 = temp = 0.0;
+    fund = *(double *)argv;
+    freqs = data + n;
 
-    for(n = 0; n < N; n++)
+    for(i = 0; i < n; i++)
     {
-        if((temp = data[n]))
+        if((temp = data[i]))
         {
             den += temp;
-            if(!p1)
-                p1 = temp;
+            h = round(freqs[i] / fund);
+            if(h == 1)
+                p1 += temp;
         }
     }
 
@@ -332,23 +333,36 @@ int xtract_tristimulus_1(const double *data, const int N, const void *argv, doub
 int xtract_tristimulus_2(const double *data, const int N, const void *argv, double *result)
 {
 
-    int n = N;
+    int n = N >> 1, h = 0, i;
+    double den, p2, p3, p4, ps, fund, temp;
+    den = p2 = p3 = p4 = ps = fund = temp = 0.0;
+    const double *freqs;
 
-    double den, p2, p3, p4, ps, temp;
-
-    den = p2 = p3 = p4 = ps = temp = 0.0;
+    fund = *(double *)argv;
+    freqs = data + n;
 
-    for(n = 0; n < N; n++)
+    for(i = 0; i < n; i++)
     {
-        if((temp = data[n]))
+        if((temp = data[i]))
         {
             den += temp;
-            if(!p2)
-                p2 = temp;
-            else if(!p3)
-                p3 = temp;
-            else if(!p4)
-                p4 = temp;
+            h = round(freqs[i] / fund);
+            switch (h)
+            {
+                case 2:
+                    p2 += temp;
+                    break;
+
+                case 3:
+                    p3 += temp;
+                    break;
+
+                case 4:
+                    p4 += temp;
+
+                default:
+                    break;
+            }
         }
     }
 
@@ -369,21 +383,21 @@ int xtract_tristimulus_2(const double *data, const int N, const void *argv, doub
 
 int xtract_tristimulus_3(const double *data, const int N, const void *argv, double *result)
 {
+    int n = N >> 1, h = 0, i;
+    double den = 0.0, num = 0.0, fund = 0.0, temp = 0.0;
+    const double *freqs;
 
-    int n = N, count = 0;
-
-    double den, num, temp;
-
-    den = num = temp = 0.0;
+    fund = *(double *)argv;
+    freqs = data + n;
 
-    for(n = 0; n < N; n++)
+    for(i = 0; i < n; i++)
     {
-        if((temp = data[n]))
+        if((temp = data[i]))
         {
             den += temp;
-            if(count >= 5)
+            h = round(freqs[i] / fund);
+            if(h >= 5)
                 num += temp;
-            count++;
         }
     }
 
@@ -622,8 +636,7 @@ int xtract_rms_amplitude(const double *data, const int N, const void *argv, doub
 
 int xtract_spectral_inharmonicity(const double *data, const int N, const void *argv, double *result)
 {
-
-    int n = N >> 1;
+    int n = N >> 1, h = 0;
     double num = 0.0, den = 0.0, fund;
     const double *freqs, *amps;
 
@@ -635,7 +648,8 @@ int xtract_spectral_inharmonicity(const double *data, const int N, const void *a
     {
         if(amps[n])
         {
-            num += fabs(freqs[n] - n * fund) * XTRACT_SQ(amps[n]);
+            h = round(freqs[n] / fund);
+            num += fabs(freqs[n] - h * fund) * XTRACT_SQ(amps[n]);
             den += XTRACT_SQ(amps[n]);
         }
     }
@@ -655,16 +669,19 @@ int xtract_power(const double *data, const int N, const void *argv, double *resu
 
 int xtract_odd_even_ratio(const double *data, const int N, const void *argv, double *result)
 {
+    int n = N >> 1, h = 0;
+    double odd = 0.0, even = 0.0, fund, temp;
+    const double *freqs;
 
-    int M = (N >> 1), n;
-
-    double odd = 0.0, even = 0.0,  temp;
+    fund = *(double *)argv;
+    freqs = data + n;
 
-    for(n = 0; n < M; n++)
+    while(n--)
     {
         if((temp = data[n]))
         {
-            if(XTRACT_IS_ODD(n))
+            h = round(freqs[n] / fund);
+            if(XTRACT_IS_ODD(h))
             {
                 odd += temp;
             }
diff --git a/xtract/xtract_scalar.h b/xtract/xtract_scalar.h
index 1c79cc9..8e24284 100644
--- a/xtract/xtract_scalar.h
+++ b/xtract/xtract_scalar.h
@@ -182,9 +182,9 @@ int xtract_irregularity_j(const double *data, const int N, const void *argv, dou
 
 /** \brief Calculate the Tristimulus of an input vector using a method described by Pollard and Jansson (1982)
  * 
- * \param *data: a pointer to the first element in an array of doubles representing the magnitude coefficients of the harmonic spectrum of an audio vector e.g. a pointer to the first half of the array pointed to by *result from xtract_harmonics(). The amplitudes of the peak spectrum (e.g. *result from xtract_peak_spectrum()) can be used if one wishes to consider all partials not just harmonics.
+ * \param *data: a pointer to the first element in an array of doubles representing a harmonic spectrum of size N/2, and a frequency spectrum of size N/2 (This is the output format of xtract_harmonic_spectrum())
  * \param N: the number of elements to be considered
- * \param *argv: a pointer to NULL
+ * \param *argv: a pointer to a double representing the fundamental frequency of the input vector.
  * \param *result: the tristimulus of N values from the array pointed to by *data
  *
  * These three functions provide the first, second and third order tristimulus formulae
@@ -324,9 +324,9 @@ int xtract_power(const double *data, const int N, const void *argv, double *resu
 /* Odd to even harmonic ratio */
 /** \brief Extract the Odd to even harmonic ratio of an input vector 
  * 
- * \param *data: a pointer to the first element in an array of doubles 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 *data: a pointer to the first element in an array of doubles representing a harmonic spectrum of size N/2, and a frequency spectrum of size N/2 (This is the output format of xtract_harmonic_spectrum())
+ * \param N: the number of elements to be considered. 
+ * \param *argv: a pointer to a double representing the fundamental frequency of the input vector.
  * \param *result: the even/odd harmonic ratio of N values from the array pointed to by *data
  */
 int xtract_odd_even_ratio(const double *data, const int N, const void *argv, double *result);