comparison src/vector.c @ 111:c7e1fe63eedb

- Re-factoring in xtract_spectrum and fixed normalisation bug
- Fixed bug in xtract_lnorm

 
 int xtract_spectrum(const float *data, const int N, const void *argv, float *result){
 
     float *input, *rfft, q, temp, max;
     size_t bytes;
-    int n, 
+    int n,
+        m,
         NxN, 
         M, 
         vector, 
         withDC, 
         argc,

 		if ((temp = XTRACT_SQ(rfft[n]) + 
 			    XTRACT_SQ(rfft[N - n])) > XTRACT_LOG_LIMIT)
 		    temp = log(sqrt(temp) / N);
 		else
 		    temp = XTRACT_LOG_LIMIT_DB;
-		if(withDC) {
-		    result[n] = 
-			/*Normalise*/
-			(temp + XTRACT_DB_SCALE_OFFSET) / 
-			XTRACT_DB_SCALE_OFFSET; 
-		    result[M + n + 1] = n * q;
-		}
-		else {
-		    result[n - 1] =
-			(temp + XTRACT_DB_SCALE_OFFSET) / 
-			XTRACT_DB_SCALE_OFFSET; 
-		    result[M + n - 1] = n * q;
-		}
-                max = result[n] > max ? result[n] : max;
+
+                if(withDC){
+                    m = n;
+		    result[M + m + 1] = n * q;
+                }
+                else{
+                    m = n - 1;
+		    result[M + m] = n * q;
+                }
+                
+                result[m] = 
+                    /* Scaling */
+                    (temp + XTRACT_DB_SCALE_OFFSET) / 
+                    XTRACT_DB_SCALE_OFFSET; 
+		
+                max = result[m] > max ? result[m] : max;
 	    }
 	    break;
 
 	case XTRACT_POWER_SPECTRUM:
 	    for(n = 1; n < M; n++){
-		if(withDC){
-		    result[n] = (XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) 
-			/ NxN;
-		    result[M + n  + 1] = n * q;
-		}
-		else {
-		    result[n - 1] = 
-			(XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) / NxN;
-		    result[M + n - 1] = n * q;
-		}
-                max = result[n] > max ? result[n] : max;
+                if(withDC){
+                    m = n;
+		    result[M + m + 1] = n * q;
+                }
+                else{
+                    m = n - 1;
+		    result[M + m] = n * q;
+                }
+                result[m] = (XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) / NxN;
+                max = result[m] > max ? result[m] : max;
 	    }
 	    break;
 
 	case XTRACT_LOG_POWER_SPECTRUM:
 	    for(n = 1; n < M; n++){
 		if ((temp = XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) > 
 			XTRACT_LOG_LIMIT)
 		    temp = log(temp / NxN);
 		else
 		    temp = XTRACT_LOG_LIMIT_DB; 		
-		if(withDC){
-		    result[n] = (temp + XTRACT_DB_SCALE_OFFSET) / 
+
+                if(withDC){
+                    m = n;
+		    result[M + m + 1] = n * q;
+                }
+                else{
+                    m = n - 1;
+		    result[M + m] = n * q;
+                }
+
+                result[m] = (temp + XTRACT_DB_SCALE_OFFSET) / 
 			XTRACT_DB_SCALE_OFFSET; 
-		    result[M + n + 1] = n * q;
-		}
-		else {
-		    result[n - 1] = (temp + XTRACT_DB_SCALE_OFFSET) / 
-			XTRACT_DB_SCALE_OFFSET; 
-		    result[M + n - 1] = n * q;
-		}
-                max = result[n] > max ? result[n] : max;
+                max = result[m] > max ? result[m] : max;
 	    }
 	    break;
 
 	default:
 	    /* MAGNITUDE_SPECTRUM */
 	    for(n = 1; n < M; n++){
-		if(withDC){
-		    result[n] = sqrt(XTRACT_SQ(rfft[n]) + 
-			    XTRACT_SQ(rfft[N - n])) / N; 
-		    result[M + n + 1] = n * q;
-		}
-		else {
-		    result[n - 1] = sqrt(XTRACT_SQ(rfft[n]) + 
-			    XTRACT_SQ(rfft[N - n])) / N; 
-		    result[M + n - 1] = n * q;
-		}
-                max = result[n] > max ? result[n] : max;
+                if(withDC){
+                    m = n;
+		    result[M + m + 1] = n * q;
+                }
+                else{
+                    m = n - 1;
+		    result[M + m] = n * q;
+                }
+
+                result[m] = sqrt(XTRACT_SQ(rfft[n]) + 
+                        XTRACT_SQ(rfft[N - n])) / N; 
+                max = result[m] > max ? result[m] : max;
 	    }
 	    break;
     }
     
     if(withDC){

         max = result[0] > max ? result[0] : max;
 	/* The Nyquist */ 
 	result[M] = XTRACT_SQ(rfft[M]);
 	result[N + 1] = q * M;
         max = result[M] > max ? result[M] : max;
-        M++; /* So we normalise the Nyquist (below) */
     }
     else {
 	/* The Nyquist */ 
 	result[M - 1] = (float)XTRACT_SQ(rfft[M]);
 	result[N - 1] = q * M;
         max = result[M - 1] > max ? result[M - 1] : max;
     }
-    
 
     if(normalise){
         for(n = 0; n < M; n++)
             result[n] /= max;
     }