Improve spectrum testing for recording

- Use PSD of each signals as base line for comparison:
  ratio = sqrt(PSD_DUT / PSD_host)
  The ratio is normalized for its mean value for the checked freq range.
- random signal generator changed to have only max frequency

Change-Id: I3639d70c4b6159061ada6ca9ac8f43e48f559275
diff --git a/suite/audio_quality/test_description/dut_recording_spectrum.xml b/suite/audio_quality/test_description/dut_recording_spectrum.xml
index 0f58deb..2bfe661 100644
--- a/suite/audio_quality/test_description/dut_recording_spectrum.xml
+++ b/suite/audio_quality/test_description/dut_recording_spectrum.xml
@@ -17,22 +17,21 @@
 
 <case name="dut_recording_spectrum" version="1.0" description="Check frequency spectrum for recording">
 	<setup>
-		<!-- input: peak amplitude, duration in msec, sampling rate, low frequency, high frequency, output: generated sound-->
-		<process method="script:gen_random" input="consti:10000,consti:30000,consti:44100,consti:50,consti:6000" output="id:sound1" />
-		<!-- sound id="sound1" type="random:10000:30000" / -->
+		<!-- input: peak amplitude, duration in msec, sampling rate, high frequency, output: generated sound-->
+		<process method="script:gen_random" input="consti:10000,consti:30000,consti:44100,consti:6000" output="id:sound1" />
 		<!--  Only for starting client app early. The data is not used -->
 		<sound id="sound2" type="sin:1:1000:2" preload="1"/>
 	</setup>
 	<action>
 		<sequential repeat="1" index="i">
 			<output device="host" id="sound1" gain="100" sync="start" waitforcompletion="0" />
-			<sequential repeat="1" index="j">
+			<sequential repeat="5" index="j">
 				<input device="host" id="host_in_$j" gain="100" time="6000" sync="start" />
 				<input device="DUT" id="dut_in_$j" gain="100" time="4000" sync="start" />
 			</sequential>
 		</sequential>
-		<sequential repeat="1" index="k">
-			<!-- input: host record, device record, samping rate, low frequency in Hz, high frequency in Hz, allowed error for pass in smaller side, allowed error in bigger side%, output: min value in lower side calculated normalized to 1.0, max value in higher side, calculated TF in mannitude only between low f to high f -->
+		<sequential repeat="5" index="k">
+			<!-- input: host record, device record, samping rate, low frequency in Hz, high frequency in Hz, allowed error for pass in smaller side, allowed error in bigger side%, output: min value in lower side calculated normalized to 1.0, max value in higher side, calculated amplitude ratio in mannitude only between low f to high f -->
 			<process method="script:check_spectrum" input="id:host_in_$k,id:dut_in_$k,consti:44100,consti:200,consti:4000,constf:50.0,constf:100.0" output="val:min_val_$k,val:max_val_$k,id:tf_$k" />
 		</sequential>
 	</action>
diff --git a/suite/audio_quality/test_description/processing/check_spectrum.py b/suite/audio_quality/test_description/processing/check_spectrum.py
index 31b5a5e..913b88c 100644
--- a/suite/audio_quality/test_description/processing/check_spectrum.py
+++ b/suite/audio_quality/test_description/processing/check_spectrum.py
@@ -23,7 +23,7 @@
 sys.path.append(sys.path[0])
 import calc_delay
 
-# check if Transfer Function of DUT / Host signal
+# check if amplitude ratio of DUT / Host signal
 #  lies in the given error boundary
 # input: host record
 #        device record,
@@ -34,35 +34,42 @@
 #        allowed error ih positive side for pass
 # output: min value in negative side, normalized to 1.0
 #         max value in positive side
-#         calculated TF in magnitude (DUT / Host)
+#         calculated amplittude ratio in magnitude (DUT / Host)
 
 def do_check_spectrum(hostData, DUTData, samplingRate, fLow, fHigh, margainLow, margainHigh):
     # reduce FFT resolution to have averaging effects
-    N = 512 if (len(hostData) > 512) else len(hostData)
-    iLow = N * fLow / samplingRate
+    N = 256 if (len(hostData) > 512) else len(hostData)
+    iLow = N * fLow / samplingRate + 1 # 1 for DC
     if iLow > (N / 2 - 1):
         iLow = (N / 2 - 1)
-    iHigh = N * fHigh / samplingRate
-    if iHigh > (N / 2):
-        iHigh = N / 2
+    iHigh = N * fHigh / samplingRate + 1 # 1 for DC
+    if iHigh > (N / 2 + 1):
+        iHigh = N / 2 + 1
     print fLow, iLow, fHigh, iHigh, samplingRate
-    hostFFT = abs(fft.fft(hostData, n = N))[iLow:iHigh]
-    dutFFT = abs(fft.fft(DUTData, n = N))[iLow:iHigh]
-    TF = dutFFT / hostFFT
-    TFmean = sum(TF) / len(TF)
-    TF = TF / TFmean # TF normalized to 1
-    positiveMax = abs(max(TF))
-    negativeMin = abs(min(TF))
+
+    Phh, freqs = plt.psd(hostData, NFFT=N, Fs=samplingRate, Fc=0, detrend=plt.mlab.detrend_none,\
+        window=plt.mlab.window_hanning, noverlap=0, pad_to=None, sides='onesided',\
+        scale_by_freq=False)
+    Pdd, freqs = plt.psd(DUTData, NFFT=N, Fs=samplingRate, Fc=0, detrend=plt.mlab.detrend_none,\
+        window=plt.mlab.window_hanning, noverlap=0, pad_to=None, sides='onesided',\
+        scale_by_freq=False)
+    print len(Phh), len(Pdd)
+    print "Phh", abs(Phh[iLow:iHigh])
+    print "Pdd", abs(Pdd[iLow:iHigh])
+    amplitudeRatio = np.sqrt(abs(Pdd[iLow:iHigh]/Phh[iLow:iHigh]))
+    ratioMean = np.mean(amplitudeRatio)
+    amplitudeRatio = amplitudeRatio / ratioMean
+    print "Normialized ratio", amplitudeRatio
+    print "ration mean for normalization", ratioMean
+    positiveMax = abs(max(amplitudeRatio))
+    negativeMin = abs(min(amplitudeRatio))
     passFail = True if (positiveMax < (margainHigh / 100.0 + 1.0)) and\
         ((1.0 - negativeMin) < margainLow / 100.0) else False
-    TFResult = np.zeros(len(TF), dtype=np.int16)
-    for i in range(len(TF)):
-        TFResult[i] = TF[i] * 256 # make fixed point
-    #freq = np.linspace(0.0, fHigh, num=iHigh, endpoint=False)
-    #plt.plot(freq, abs(fft.fft(hostData, n = N))[:iHigh], freq, abs(fft.fft(DUTData, n = N))[:iHigh])
-    #plt.show()
+    RatioResult = np.zeros(len(amplitudeRatio), dtype=np.int16)
+    for i in range(len(amplitudeRatio)):
+        RatioResult[i] = amplitudeRatio[i] * 256 # make fixed point
     print "positiveMax", positiveMax, "negativeMin", negativeMin
-    return (passFail, negativeMin, positiveMax, TFResult)
+    return (passFail, negativeMin, positiveMax, RatioResult)
 
 def check_spectrum(inputData, inputTypes):
     output = []
@@ -125,10 +132,10 @@
     samplingRate = 44100
     fLow = 500
     fHigh = 15000
-    data = getattr(mod, "do_gen_random")(peakAmpl, durationInMSec, samplingRate, fLow, fHigh,\
-                                         stereo=False)
+    data = getattr(mod, "do_gen_random")(peakAmpl, durationInMSec, samplingRate, fHigh,\
+        stereo=False)
     print len(data)
-    (passFail, minVal, maxVal, TF) = do_check_spectrum(data, data, samplingRate, fLow, fHigh,\
-                                                           1.0, 1.0)
-    plt.plot(TF)
+    (passFail, minVal, maxVal, ampRatio) = do_check_spectrum(data, data, samplingRate, fLow, fHigh,\
+        1.0, 1.0)
+    plt.plot(ampRatio)
     plt.show()
diff --git a/suite/audio_quality/test_description/processing/gen_random.py b/suite/audio_quality/test_description/processing/gen_random.py
index 541c60a..305526a 100644
--- a/suite/audio_quality/test_description/processing/gen_random.py
+++ b/suite/audio_quality/test_description/processing/gen_random.py
@@ -20,32 +20,24 @@
 import scipy.fftpack as fft
 import matplotlib.pyplot as plt
 
-# generate random signal
+# generate random signal with max freq
 # Input: peak amplitude,
 #        duration in msec,
 #        sampling rate HZ
-#        low frequency,
 #        high frequency,
 # Output: generated sound (stereo)
 
-def do_gen_random(peakAmpl, durationInMSec, samplingRate, fLow, fHigh, stereo=True):
+def do_gen_random(peakAmpl, durationInMSec, samplingRate, fHigh, stereo=True):
     samples = durationInMSec * samplingRate / 1000
     result = np.zeros(samples * 2 if stereo else samples, dtype=np.int16)
-    #randomSignal = np.random.random_integers(-peakAmpl, peakAmpl, samples)
-    randomSignal = np.random.normal(scale=peakAmpl *2 / 3, size=samples)
+    randomSignal = np.random.normal(scale = peakAmpl * 2 / 3, size=samples)
     fftData = fft.rfft(randomSignal)
     freqSamples = samples/2
-    iLow = 0 #freqSamples * fLow * 2/ samplingRate + 1
-    #if iLow > freqSamples - 1:
-    #    iLow = freqSamples - 1
     iHigh = freqSamples * fHigh * 2 / samplingRate + 1
-    #print len(randomSignal), len(fftData), fLow, iLow, fHigh, iHigh
+    #print len(randomSignal), len(fftData), fLow, fHigh, iHigh
     if iHigh > freqSamples - 1:
         iHigh = freqSamples - 1
     fftData[0] = 0 # DC
-    #for i in range(iLow - 1):
-    #    fftData[2 * i + 1 ] = 0 # Re
-    #    fftData[2 * i + 2 ] = 0 # Imag
     for i in range(iHigh, freqSamples - 1):
         fftData[ 2 * i + 1 ] = 0
         fftData[ 2 * i + 2 ] = 0
@@ -81,15 +73,13 @@
         inputError = True
     if (inputTypes[3] != TYPE_I64):
         inputError = True
-    if (inputTypes[4] != TYPE_I64):
-        inputError = True
     if inputError:
         output.append(RESULT_ERROR)
         output.append(outputData)
         output.append(outputTypes)
         return output
 
-    result = do_gen_random(inputData[0], inputData[1], inputData[2], inputData[3], inputData[4])
+    result = do_gen_random(inputData[0], inputData[1], inputData[2], inputData[3])
 
     output.append(RESULT_OK)
     outputData.append(result)
@@ -103,8 +93,8 @@
     peakAmplitude = 10000
     samplingRate = 44100
     durationInMSec = 10000
-    fLow = 500
+    #fLow = 500
     fHigh = 15000
-    result = do_gen_random(peakAmplitude, durationInMSec, samplingRate, fLow, fHigh)
+    result = do_gen_random(peakAmplitude, durationInMSec, samplingRate, fHigh)
     plt.plot(result)
     plt.show()
diff --git a/suite/audio_quality/test_description/review/review_dut_recording_spectrum.xml b/suite/audio_quality/test_description/review/review_dut_recording_spectrum.xml
new file mode 100644
index 0000000..b8157be
--- /dev/null
+++ b/suite/audio_quality/test_description/review/review_dut_recording_spectrum.xml
@@ -0,0 +1,38 @@
+<?xml version="1.0" encoding="utf-8"?>
+
+<!-- Copyright (C) 2012 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.
+-->
+
+<case name="review_dut_recording_spectrum" version="1.0" description="Review">
+	<setup>
+		<sound id="sound_host_0" type="file:reports/2012_05_14_09_27_40/dut_recording_spectrum/host_in_0.r2m" />
+		<sound id="sound_dut_0" type="file:reports/2012_05_14_09_27_40/dut_recording_spectrum/dut_in_0.r2m" />
+		<sound id="sound_host_1" type="file:reports/2012_05_14_09_27_40/dut_recording_spectrum/host_in_1.r2m" />
+		<sound id="sound_dut_1" type="file:reports/2012_05_14_09_27_40/dut_recording_spectrum/dut_in_1.r2m" />
+		<sound id="sound_host_2" type="file:reports/2012_05_14_09_27_40/dut_recording_spectrum/host_in_2.r2m" />
+		<sound id="sound_dut_2" type="file:reports/2012_05_14_09_27_40/dut_recording_spectrum/dut_in_2.r2m" />
+		<sound id="sound_host_3" type="file:reports/2012_05_14_09_27_40/dut_recording_spectrum/host_in_3.r2m" />
+		<sound id="sound_dut_3" type="file:reports/2012_05_14_09_27_40/dut_recording_spectrum/dut_in_3.r2m" />
+		<sound id="sound_host_4" type="file:reports/2012_05_14_09_27_40/dut_recording_spectrum/host_in_4.r2m" />
+		<sound id="sound_dut_4" type="file:reports/2012_05_14_09_27_40/dut_recording_spectrum/dut_in_4.r2m" />
+	</setup>
+	<action>
+		<sequential repeat="5" index="k">
+			<!-- input: host record, device record, samping rate, low frequency in Hz, high frequency in Hz, allowed error for pass in smaller side, allowed error in bigger side%, output: min value in lower side calculated normalized to 1.0, max value in higher side, calculated TF in mannitude only between low f to high f -->
+			<process method="script:check_spectrum" input="id:sound_host_$k,id:sound_dut_$k,consti:44100,consti:200,consti:4000,constf:50.0,constf:100.0" output="val:min_val_$k,val:max_val_$k,id:tf_$k" />
+		</sequential>
+	</action>
+	<save report="min_val_.*,max_val_.*" />
+</case>