comparison src/vector.c @ 52:45c585bb7996

Rationalised spectral data format. Added spectral_mean et al

 
 #include <fftw3.h>
 
 int xtract_magnitude_spectrum(const float *data, const int N, const void *argv, float *result){
 
-    float *temp, *input, q, sr;
+    float *temp, *input, q;
     size_t bytes;
     int n , M = N >> 1;
     fftwf_plan plan;
 
     temp = (float *)fftwf_malloc(N * sizeof(float));
     input = (float *)malloc(bytes = N * sizeof(float));
     input = memcpy(input, data, bytes);
 
-    q = sr = 0.f;
-
-    sr = *(float *)argv;
-
-    CHECK_SR;
-
-    q = (sr * .5) / M;
+    q = *(float *)argv;
+
+    CHECK_q;
 
     plan = fftwf_plan_r2r_1d(N, input, temp, FFTW_R2HC, FFTW_ESTIMATE);
     
     fftwf_execute(plan);
     
     for(n = 1; n < M; n++){
-        result[M + n] = sqrt(SQ(temp[n]) + SQ(temp[N - n])) / N;
-        result[n] = n * q;
-    }
-    
-    result[M] = fabs(temp[0]) / N;
-    result[0] = q * .5;
+        result[n] = sqrt(SQ(temp[n]) + SQ(temp[N - n])) / N;
+        result[M + n] = n * q;
+    }
+    
+    result[0] = fabs(temp[0]) / N;
+    result[M] = q * .5;
     
     fftwf_destroy_plan(plan);
     fftwf_free(temp);
     free(input);
     

     }
 
     return SUCCESS;
 }
 
-int xtract_peaks(const float *data, const int N, const void *argv, float *result){
-
-    float thresh, max, y, y2, 
-	  y3, p, width, sr,
-	  *input = NULL;
+int xtract_peak_spectrum(const float *data, const int N, const void *argv, float *result){
+
+    float thresh, max, y, y2, y3, p, q, *input = NULL;
     size_t bytes;
     int n = N, M, rv = SUCCESS;
 
-    thresh = max = y = y2 = y3 = p = width = sr = 0.f;
+    thresh = max = y = y2 = y3 = p = q = 0.f;
     
     if(argv != NULL){
         thresh = ((float *)argv)[0];
-        sr = ((float *)argv)[1];
+        q = ((float *)argv)[1];
     }
     else
         rv = BAD_ARGV;
 
     if(thresh < 0 || thresh > 100){
         thresh = 0;
         rv = BAD_ARGV;
     }
 
-    CHECK_SR;
+    CHECK_q;
 
     input = (float *)malloc(bytes = N * sizeof(float));
 
     if(input != NULL)
 	input = memcpy(input, data, bytes);
     else
 	return MALLOC_FAILED;
 
     M = N >> 1;
-    width = sr / N;
 
     while(n--)
         max = MAX(max, input[n]);
     
     thresh *= .01 * max;

     result[M] = 0;
 
     for(n = 1; n < M; n++){
         if(input[n] >= thresh){
             if(input[n] > input[n - 1] && input[n] > input[n + 1]){
-                result[n] = width * (n + (p = .5 * (y = input[n-1] - (y3 = input[n+1])) / (input[n - 1] - 2 * (y2 = input[n]) + input[n + 1])));
-                result[M + n] = y2 - .25 * (y - y3) * p;
+                result[M + n] = q * (n + (p = .5 * (y = input[n-1] - 
+				(y3 = input[n+1])) / (input[n - 1] - 2 * 
+				    (y2 = input[n]) + input[n + 1])));
+                result[n] = y2 - .25 * (y - y3) * p;
             }
             else{
                 result[n] = 0;
                 result[M + n] = 0;
             }

    
     free(input);
     return (rv ? rv : SUCCESS);
 }
 	    
-int xtract_harmonics(const float *data, const int N, const void *argv, float *result){
+int xtract_harmonic_spectrum(const float *data, const int N, const void *argv, float *result){
     
     int n = (N >> 1), M = n; 
 
     const float *freqs, *amps;
     float f0, thresh, ratio, nearest, distance;
 
-    freqs = data;
-    amps = data + n;
+    amps = data;
+    freqs = data + n;
     f0 = *((float *)argv);
     thresh = *((float *)argv+1);
 
     ratio = nearest = distance = 0.f;
 

 	    nearest = round(ratio);
 	    distance = fabs(nearest - ratio);
 	    if(distance > thresh)
 		result[n] = result[M + n] = 0.f;
 	    else {
-		result[n] = freqs[n];
-		result[M + n] = amps[n];
+		result[n] = amps[n];
+		result[M + n] = freqs[n];
 	    }
 	}
 	else
 	    result[n] = result[M + n] = 0.f;
     }