change LinearityTest algorithm to use RMS of the whole signal
and use linear regression to fit the result
volume level is fixed to 2/3 of max

Change-Id: I1b5b2130a6786c588ce47f1e2306db37539289d6
diff --git a/apps/CtsVerifier/jni/audioquality/Android.mk b/apps/CtsVerifier/jni/audioquality/Android.mk
index 8e827f7..787128a 100644
--- a/apps/CtsVerifier/jni/audioquality/Android.mk
+++ b/apps/CtsVerifier/jni/audioquality/Android.mk
@@ -20,11 +20,13 @@
 
 LOCAL_MODULE_TAGS := optional
 
-
+LOCAL_SHARED_LIBRARIES := \
+	libcutils
 
 LOCAL_MODULE    := libaudioquality
 LOCAL_SRC_FILES := Fft.cpp Window.cpp GlitchTest.cpp MeasureRms.cpp \
    OverflowCheck.cpp LinearityTest.cpp CompareSpectra.cpp \
-   GenerateSinusoid.cpp Wrapper.cpp
+   GenerateSinusoid.cpp Wrapper.cpp LinearityTestRms.cpp
+
 
 include $(BUILD_SHARED_LIBRARY)
diff --git a/apps/CtsVerifier/jni/audioquality/LinearityTest.h b/apps/CtsVerifier/jni/audioquality/LinearityTest.h
index d92566a..1aac4b6 100644
--- a/apps/CtsVerifier/jni/audioquality/LinearityTest.h
+++ b/apps/CtsVerifier/jni/audioquality/LinearityTest.h
@@ -17,6 +17,22 @@
 #ifndef LINEARITY_TEST_H
 #define LINEARITY_TEST_H
 
+/* error codes returned */
+// The input signals or sample counts are missing.
+const int ERROR_INPUT_SIGNAL_MISSING       = -1;
+// The number of input signals is < 2.
+const int ERROR_INPUT_SIGNAL_NUMBERS       = -2;
+// The specified sample rate is <= 4000.0
+const int ERROR_SAMPLE_RATE_TOO_LOW        = -3;
+// The dB step size for the increase in stimulus level is <= 0.0
+const int ERROR_NEGATIVE_STEP_SIZE         = -4;
+// The specified reference stimulus number is out of range.
+const int ERROR_STIMULUS_NUMBER            = -5;
+// One or more of the stimuli is too short in duration.
+const int ERROR_STIMULI_TOO_SHORT          = -6;
+// cannot find linear fitting for the given data
+const int ERROR_LINEAR_FITTING             = -7;
+
 /* There are numSignals int16 signals in pcms.  sampleCounts is an
    integer array of length numSignals containing their respective
    lengths in samples.  They are all sampled at sampleRate.  The pcms
@@ -27,16 +43,15 @@
    normal speaking levels).  The maximum deviation in linearity found
    (in dB) is returned in maxDeviation.  The function returns 1 if
    the measurements could be made, or a negative number that
-   indicates the error, as follows:
-      -1 The input signals or sample counts are missing.
-      -2 The number of input signals is < 2.
-      -3 The specified sample rate is <= 4000.0
-      -4 The dB step size for the increase in stimulus level is <= 0.0/
-      -5 The specified reverence stimulus number is out of range.
-      -6 One or more of the stimuli is too short in duration. */
+   indicates the error, as defined above. */
 int linearityTest(short** pcms, int* sampleCounts, int numSignals,
                   float sampleRate, float dbStepSize,
                   int referenceStim, float* maxDeviation);
 
+/* a version using RMS of the whole signal and linear fitting */
+int linearityTestRms(short** pcms, int* sampleCounts, int numSignals,
+                  float sampleRate, float dbStepSize,
+                  float* maxDeviation);
+
 
 #endif /* LINEARITY_TEST_H */
diff --git a/apps/CtsVerifier/jni/audioquality/LinearityTestRms.cpp b/apps/CtsVerifier/jni/audioquality/LinearityTestRms.cpp
new file mode 100644
index 0000000..cf47ec7
--- /dev/null
+++ b/apps/CtsVerifier/jni/audioquality/LinearityTestRms.cpp
@@ -0,0 +1,285 @@
+/*
+ * Copyright (C) 2011 The Android Open Source Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License. You may obtain a copy of
+ * the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+ * License for the specific language governing permissions and limitations under
+ * the License.
+ */
+
+/* This test accepts a collection of N speech waveforms collected as
+   part of N recognition attempts.  The waveforms are ordered by
+   increasing presentation level.  The test determines the extent to
+   which the peak amplitudes in the waveforms track the change in
+   presentation level.  Failure to track the presentation level within
+   some reasonable margin is an indication of clipping or of automatic
+   gain control in the signal path.
+
+   RMS of each level is used as a parameter for deciding lienairy.
+   For each level, RMS is calculated, and a line fitting into RMS vs level
+   will be calculated. Then for each level, residual error of measurement
+   vs line fitting will be calculated, and the residual error is normalized
+   with each measurement. The test failes if residual error is bigger than
+   2dB.
+
+   This test will be robust to background noise as long as it is persistent.
+   But background noise which appears shortly with enough strength can
+   mess up the result.
+*/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <cutils/log.h>
+#include "LinearityTest.h"
+
+#define LOG_TAG "LinearityTestRms"
+
+// vectorDot, vectorNorm, and solveLeastSquares
+// copied from frameworks/base/libs/ui/Input.cpp
+
+static inline float vectorDot(const float* a, const float* b, uint32_t m) {
+    float r = 0;
+    while (m--) {
+        r += *(a++) * *(b++);
+    }
+    return r;
+}
+
+static inline float vectorNorm(const float* a, uint32_t m) {
+    float r = 0;
+    while (m--) {
+        float t = *(a++);
+        r += t * t;
+    }
+    return sqrtf(r);
+}
+
+/**
+ * Solves a linear least squares problem to obtain a N degree polynomial that fits
+ * the specified input data as nearly as possible.
+ *
+ * Returns true if a solution is found, false otherwise.
+ *
+ * The input consists of two vectors of data points X and Y with indices 0..m-1.
+ * The output is a vector B with indices 0..n-1 that describes a polynomial
+ * that fits the data, such the sum of abs(Y[i] - (B[0] + B[1] X[i] + B[2] X[i]^2 ... B[n] X[i]^n))
+ * for all i between 0 and m-1 is minimized.
+ *
+ * That is to say, the function that generated the input data can be approximated
+ * by y(x) ~= B[0] + B[1] x + B[2] x^2 + ... + B[n] x^n.
+ *
+ * The coefficient of determination (R^2) is also returned to describe the goodness
+ * of fit of the model for the given data.  It is a value between 0 and 1, where 1
+ * indicates perfect correspondence.
+ *
+ * This function first expands the X vector to a m by n matrix A such that
+ * A[i][0] = 1, A[i][1] = X[i], A[i][2] = X[i]^2, ..., A[i][n] = X[i]^n.
+ *
+ * Then it calculates the QR decomposition of A yielding an m by m orthonormal matrix Q
+ * and an m by n upper triangular matrix R.  Because R is upper triangular (lower
+ * part is all zeroes), we can simplify the decomposition into an m by n matrix
+ * Q1 and a n by n matrix R1 such that A = Q1 R1.
+ *
+ * Finally we solve the system of linear equations given by R1 B = (Qtranspose Y)
+ * to find B.
+ *
+ * For efficiency, we lay out A and Q column-wise in memory because we frequently
+ * operate on the column vectors.  Conversely, we lay out R row-wise.
+ *
+ * http://en.wikipedia.org/wiki/Numerical_methods_for_linear_least_squares
+ * http://en.wikipedia.org/wiki/Gram-Schmidt
+ */
+static bool solveLeastSquares(const float* x, const float* y, uint32_t m, uint32_t n,
+        float* outB, float* outDet) {
+#if DEBUG_LEAST_SQUARES
+    LOGD("solveLeastSquares: m=%d, n=%d, x=%s, y=%s", int(m), int(n),
+            vectorToString(x, m).string(), vectorToString(y, m).string());
+#endif
+
+    // Expand the X vector to a matrix A.
+    float a[n][m]; // column-major order
+    for (uint32_t h = 0; h < m; h++) {
+        a[0][h] = 1;
+        for (uint32_t i = 1; i < n; i++) {
+            a[i][h] = a[i - 1][h] * x[h];
+        }
+    }
+#if DEBUG_LEAST_SQUARES
+    LOGD("  - a=%s", matrixToString(&a[0][0], m, n, false /*rowMajor*/).string());
+#endif
+
+    // Apply the Gram-Schmidt process to A to obtain its QR decomposition.
+    float q[n][m]; // orthonormal basis, column-major order
+    float r[n][n]; // upper triangular matrix, row-major order
+    for (uint32_t j = 0; j < n; j++) {
+        for (uint32_t h = 0; h < m; h++) {
+            q[j][h] = a[j][h];
+        }
+        for (uint32_t i = 0; i < j; i++) {
+            float dot = vectorDot(&q[j][0], &q[i][0], m);
+            for (uint32_t h = 0; h < m; h++) {
+                q[j][h] -= dot * q[i][h];
+            }
+        }
+
+        float norm = vectorNorm(&q[j][0], m);
+        if (norm < 0.000001f) {
+            // vectors are linearly dependent or zero so no solution
+#if DEBUG_LEAST_SQUARES
+            LOGD("  - no solution, norm=%f", norm);
+#endif
+            return false;
+        }
+
+        float invNorm = 1.0f / norm;
+        for (uint32_t h = 0; h < m; h++) {
+            q[j][h] *= invNorm;
+        }
+        for (uint32_t i = 0; i < n; i++) {
+            r[j][i] = i < j ? 0 : vectorDot(&q[j][0], &a[i][0], m);
+        }
+    }
+#if DEBUG_LEAST_SQUARES
+    LOGD("  - q=%s", matrixToString(&q[0][0], m, n, false /*rowMajor*/).string());
+    LOGD("  - r=%s", matrixToString(&r[0][0], n, n, true /*rowMajor*/).string());
+
+    // calculate QR, if we factored A correctly then QR should equal A
+    float qr[n][m];
+    for (uint32_t h = 0; h < m; h++) {
+        for (uint32_t i = 0; i < n; i++) {
+            qr[i][h] = 0;
+            for (uint32_t j = 0; j < n; j++) {
+                qr[i][h] += q[j][h] * r[j][i];
+            }
+        }
+    }
+    LOGD("  - qr=%s", matrixToString(&qr[0][0], m, n, false /*rowMajor*/).string());
+#endif
+
+    // Solve R B = Qt Y to find B.  This is easy because R is upper triangular.
+    // We just work from bottom-right to top-left calculating B's coefficients.
+    for (uint32_t i = n; i-- != 0; ) {
+        outB[i] = vectorDot(&q[i][0], y, m);
+        for (uint32_t j = n - 1; j > i; j--) {
+            outB[i] -= r[i][j] * outB[j];
+        }
+        outB[i] /= r[i][i];
+    }
+#if DEBUG_LEAST_SQUARES
+    LOGD("  - b=%s", vectorToString(outB, n).string());
+#endif
+
+    // Calculate the coefficient of determination as 1 - (SSerr / SStot) where
+    // SSerr is the residual sum of squares (squared variance of the error),
+    // and SStot is the total sum of squares (squared variance of the data).
+    float ymean = 0;
+    for (uint32_t h = 0; h < m; h++) {
+        ymean += y[h];
+    }
+    ymean /= m;
+
+    float sserr = 0;
+    float sstot = 0;
+    for (uint32_t h = 0; h < m; h++) {
+        float err = y[h] - outB[0];
+        float term = 1;
+        for (uint32_t i = 1; i < n; i++) {
+            term *= x[h];
+            err -= term * outB[i];
+        }
+        sserr += err * err;
+        float var = y[h] - ymean;
+        sstot += var * var;
+    }
+    *outDet = sstot > 0.000001f ? 1.0f - (sserr / sstot) : 1;
+#if DEBUG_LEAST_SQUARES
+    LOGD("  - sserr=%f", sserr);
+    LOGD("  - sstot=%f", sstot);
+    LOGD("  - det=%f", *outDet);
+#endif
+    return true;
+}
+
+/* calculate RMS of given sample with numSamples of length */
+float calcRms(short* pcm, int numSamples)
+{
+    float energy = 0.0f;
+    for(int i = 0; i < numSamples; i++) {
+        float sample = (float)pcm[i];
+        energy += (sample * sample);
+    }
+    return sqrtf(energy);
+}
+
+/* There are numSignals int16 signals in pcms.  sampleCounts is an
+   integer array of length numSignals containing their respective
+   lengths in samples.  They are all sampled at sampleRate.  The pcms
+   are ordered by increasing stimulus level.  The level steps between
+   successive stimuli were of size dbStepSize dB.
+   The maximum deviation in linearity found
+   (in dB) is returned in maxDeviation.  The function returns 1 if
+   the measurements could be made, or a negative number that
+   indicates the error, as defined in LinearityTest.h */
+int linearityTestRms(short** pcms, int* sampleCounts, int numSignals,
+                  float sampleRate, float dbStepSize,
+                  float* maxDeviation) {
+    if (!(pcms && sampleCounts)) {
+        return ERROR_INPUT_SIGNAL_MISSING;
+    }
+    if (numSignals < 2) {
+      return ERROR_INPUT_SIGNAL_NUMBERS;
+    }
+    if (sampleRate <= 4000.0) {
+      return ERROR_SAMPLE_RATE_TOO_LOW;
+    }
+    if (dbStepSize <= 0.0) {
+        return ERROR_NEGATIVE_STEP_SIZE;
+    }
+
+    float* levels = new float[numSignals];
+    levels[0] = 1.0f;
+    float stepInMag = powf(10.0f, dbStepSize/20.0f);
+    for(int i = 1; i < numSignals; i++) {
+        levels[i] = levels[i - 1] * stepInMag;
+    }
+
+    float* rmsValues = new float[numSignals];
+    for (int i = 0; i < numSignals; i++) {
+        rmsValues[i] = calcRms(pcms[i], sampleCounts[i]);
+    }
+    const int NO_COEFFS = 2; // only line fitting
+    float coeffs[NO_COEFFS];
+    float outDet;
+    if(!solveLeastSquares(levels, rmsValues, numSignals, NO_COEFFS,
+                          coeffs, &outDet)) {
+        LOGI(" solveLeastSquares fails with det %f", outDet);
+        return ERROR_LINEAR_FITTING;
+    }
+    LOGI(" coeffs offset %f linear %f", coeffs[0], coeffs[1]);
+    float maxDev = 0.0f;
+    for(int i = 0; i < numSignals; i++) {
+        float residue = coeffs[0] + coeffs[1] * levels[i] - rmsValues[i];
+        // to make devInDb positive, add measured value itself
+        // then normalize
+        float devInDb = 20.0f * log10f((fabs(residue) + rmsValues[i])
+                                       / rmsValues[i]);
+        LOGI(" %d-th residue %f dB", i, devInDb);
+        if (devInDb > maxDev) {
+            maxDev = devInDb;
+        }
+    }
+    *maxDeviation = maxDev;
+
+    delete[] levels;
+    delete[] rmsValues;
+
+    return 1;
+}
diff --git a/apps/CtsVerifier/jni/audioquality/Wrapper.cpp b/apps/CtsVerifier/jni/audioquality/Wrapper.cpp
index bee15c6..438f0d89 100644
--- a/apps/CtsVerifier/jni/audioquality/Wrapper.cpp
+++ b/apps/CtsVerifier/jni/audioquality/Wrapper.cpp
@@ -57,6 +57,12 @@
             JNIEnv *env, jobject obj,
             jobjectArray jpcms,
             jfloat sampleRate, jfloat dbStepSize, jint referenceStim);
+
+    JNIEXPORT jfloat JNICALL
+        Java_com_android_cts_verifier_audioquality_Native_linearityTestRms(
+            JNIEnv *env, jobject obj,
+            jobjectArray jpcms,
+            jfloat sampleRate, jfloat dbStepSize);
 };
 
 /* Returns an array of sinusoidal samples.
@@ -241,3 +247,41 @@
 
     return maxDeviation;
 }
+
+
+/* Return maximum deviation from linearity in dB.
+   On failure returns:
+      -1.0 The input signals or sample counts are missing.
+      -2.0 The number of input signals is < 2.
+      -3.0 The specified sample rate is <= 4000.0
+      -4.0 The dB step size for the increase in stimulus level is <= 0.0
+      -6.0 One or more of the stimuli is too short in duration.
+*/
+JNIEXPORT jfloat JNICALL
+    Java_com_android_cts_verifier_audioquality_Native_linearityTestRms(
+        JNIEnv *env, jobject obj,
+        jobjectArray jpcms,
+        jfloat sampleRate, jfloat dbStepSize) {
+    int numSignals = env->GetArrayLength(jpcms);
+    int *sampleCounts = new int[numSignals];
+    short **pcms = new shortPtr[numSignals];
+    jshortArray ja;
+    for (int i = 0; i < numSignals; i++) {
+        ja = (jshortArray) env->GetObjectArrayElement(jpcms, i);
+        sampleCounts[i] = env->GetArrayLength(ja);
+        pcms[i] = new short[sampleCounts[i]];
+        env->GetShortArrayRegion(ja, 0, sampleCounts[i], pcms[i]);
+    }
+
+    float maxDeviation = -1.0;
+    int ret = linearityTestRms(pcms, sampleCounts, numSignals,
+            sampleRate, dbStepSize, &maxDeviation);
+    delete[] sampleCounts;
+    for (int i = 0; i < numSignals; i++) {
+        delete[] pcms[i];
+    }
+    delete[] pcms;
+    if (ret < 1) return ret;
+
+    return maxDeviation;
+}
diff --git a/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/AudioQualityVerifierActivity.java b/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/AudioQualityVerifierActivity.java
index 755b62d..af389a6 100644
--- a/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/AudioQualityVerifierActivity.java
+++ b/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/AudioQualityVerifierActivity.java
@@ -134,7 +134,7 @@
         fillAdapter();
         mList.setAdapter(mAdapter);
         mList.setOnItemClickListener(this);
-        checkNotSilent();
+        adjustVolume();
     }
 
     @Override
@@ -142,18 +142,18 @@
         super.onResume();
         mAdapter.notifyDataSetChanged(); // Update List UI
         setVolumeControlStream(AudioManager.STREAM_MUSIC);
-        checkNotSilent();
+        adjustVolume();
     }
 
-    private void checkNotSilent() {
+    private void adjustVolume() {
         AudioManager mgr = (AudioManager) getSystemService(Context.AUDIO_SERVICE);
         mgr.setStreamMute(PLAYBACK_STREAM, false);
         int volume = mgr.getStreamVolume(PLAYBACK_STREAM);
         int max = mgr.getStreamMaxVolume(PLAYBACK_STREAM);
+        int target = (max * 2) / 3;
         Log.i(TAG, "Volume " + volume + ", max " + max);
-        if (volume <= max / 10) {
-            // Volume level is silent or very quiet; increase to two-thirds
-            mgr.setStreamVolume(PLAYBACK_STREAM, (max * 2) / 3, AudioManager.FLAG_SHOW_UI);
+        if (volume != target) {
+            mgr.setStreamVolume(PLAYBACK_STREAM, target, AudioManager.FLAG_SHOW_UI);
         }
     }
 
@@ -188,6 +188,7 @@
     // Implements AdapterView.OnItemClickListener
     public void onItemClick(AdapterView<?> parent, View view, int position, long id) {
         if (mRunningExperiment) return;
+        adjustVolume();
         runExperiment(position, false);
     }
 
@@ -207,6 +208,7 @@
 
     // Implements View.OnClickListener:
     public void onClick(View v) {
+        adjustVolume();
         if (v == mCalibrateButton) {
             Intent intent = new Intent(this, CalibrateVolumeActivity.class);
             startActivity(intent);
diff --git a/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/Native.java b/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/Native.java
index 87b11d1..cdb3e25 100644
--- a/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/Native.java
+++ b/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/Native.java
@@ -31,6 +31,8 @@
             float sampleRate);
     public native float linearityTest(short[][] pcms,
         float sampleRate, float dbStepSize, int referenceStim);
+    public native float linearityTestRms(short[][] pcms,
+        float sampleRate, float dbStepSize);
 
     // The following indexes must match those in wrapper.cc
     public static final int MEASURE_RMS_RMS = 0;
diff --git a/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/experiments/GainLinearityExperiment.java b/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/experiments/GainLinearityExperiment.java
index 05f1602..225df1a 100644
--- a/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/experiments/GainLinearityExperiment.java
+++ b/apps/CtsVerifier/src/com/android/cts/verifier/audioquality/experiments/GainLinearityExperiment.java
@@ -66,9 +66,8 @@
         for (int i = 0; i < LEVELS; i++) {
             pcms[i] = Utils.byteToShortArray(record[i]);
         }
-        // We specify the middle stimulus (LEVELS / 2) as the "reference":
-        float deviation = mNative.linearityTest(pcms, AudioQualityVerifierActivity.SAMPLE_RATE,
-                DB_STEP_SIZE, LEVELS / 2);
+        float deviation = mNative.linearityTestRms(pcms, AudioQualityVerifierActivity.SAMPLE_RATE,
+                DB_STEP_SIZE);
         if (deviation < 0.0f) {
             setScore(getString(R.string.aq_fail));
             setReport(String.format(getString(R.string.aq_linearity_report_error), deviation));