1 /* 2 * Copyright (C) 2010 The Android Open Source Project 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); you may not 5 * use this file except in compliance with the License. You may obtain a copy of 6 * the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT 12 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the 13 * License for the specific language governing permissions and limitations under 14 * the License. 15 */ 16 17 #include "GlitchTest.h" 18 #include "Window.h" 19 #include "Fft.h" 20 21 #include <math.h> 22 23 GlitchTest::GlitchTest(void) : mRe(0), mIm(0), mFt(0), mWind(0) {} 24 25 void GlitchTest::init(float sampleRate, float stimFreq, float onsetThresh, 26 float dbSnrThresh) { 27 cleanup(); // in case Init gets called multiple times. 28 // Analysis parameters 29 float windowDur = 0.025; // Duration of the analysis window in seconds. 30 float frameInterval = 0.01; // Interval in seconds between frame onsets 31 float lowestFreq = 200.0; // Lowest frequency to include in 32 // background noise power integration. 33 mSampleRate = sampleRate; 34 mOnsetThresh = onsetThresh; 35 mDbSnrThresh = dbSnrThresh; 36 mFrameStep = (int)(0.5 + (mSampleRate * frameInterval)); 37 mWindowSize = (int)(0.5 + (mSampleRate * windowDur)); 38 mWind = new Window(mWindowSize); 39 mFt = new Fft; 40 mFt->fftInit(mFt->fftPow2FromWindowSize(mWindowSize)); 41 mFftSize = mFt->fftGetSize(); 42 mRe = new float[mFftSize]; 43 mIm = new float[mFftSize]; 44 float freqInterval = mSampleRate / mFftSize; 45 // We can exclude low frequencies from consideration, since the 46 // phone may have DC offset in the A/D, and there may be rumble in 47 // the testing room. 48 mLowestSpectrumBin = (int)(0.5 + (lowestFreq / freqInterval)); 49 // These are the bin indices within which most of the energy due to 50 // the (windowed) tone should be found. 51 mLowToneBin = (int)(0.5 + (stimFreq / freqInterval)) - 2; 52 mHighToneBin = (int)(0.5 + (stimFreq / freqInterval)) + 2; 53 if (mLowestSpectrumBin >= mLowToneBin) { 54 mLowestSpectrumBin = mHighToneBin + 1; 55 } 56 } 57 58 int GlitchTest::checkToneSnr(short* pcm, int numSamples, float* duration, 59 int* numBadFrames) { 60 *numBadFrames = 0; 61 *duration = 0.0; 62 if (!(mRe && mIm)) { 63 return -1; // not initialized. 64 } 65 int n_frames = 1 + ((numSamples - mWindowSize) / mFrameStep); 66 if (n_frames < 4) { // pathologically short input signal 67 return -2; 68 } 69 *numBadFrames = 0; 70 int onset = -1; 71 int offset = -1; 72 for (int frame = 0; frame < n_frames; ++frame) { 73 int numSpectra = 0; 74 mWind->window(pcm + frame*mFrameStep, mRe, 0.0); 75 realMagSqSpectrum(mRe, mWindowSize, mRe, &numSpectra); 76 int maxLoc = 0; 77 float maxValue = 0.0; 78 findPeak(mRe, mLowestSpectrumBin, numSpectra, &maxLoc, &maxValue); 79 // possible states: (1) before tone onset; (2) within tone 80 // region; (3) after tone offset. 81 if ((onset < 0) && (offset < 0)) { // (1) 82 if ((maxLoc >= mLowToneBin) && (maxLoc <= mHighToneBin) 83 && (maxValue > mOnsetThresh)) { 84 onset = frame; 85 } 86 continue; 87 } 88 if ((onset >= 0) && (offset < 0)) { // (2) 89 if (frame > (onset + 2)) { // let the framer get completely 90 // into the tonal signal 91 double sumNoise = 1.0; // init. to small non-zero vals to 92 // avoid log or divz problems 93 double sumSignal = 1.0; 94 float snr = 0.0; 95 if (maxValue < mOnsetThresh) { 96 offset = frame; 97 *duration = mFrameStep * (offset - onset) / mSampleRate; 98 if (*numBadFrames >= 1) { 99 (*numBadFrames) -= 1; // account for expected glitch at 100 // signal offset 101 } 102 continue; 103 } 104 for (int i = mLowestSpectrumBin; i < mLowToneBin; ++i) { 105 sumNoise += mRe[i]; // note that mRe contains the magnitude 106 // squared spectrum. 107 } 108 for (int i = mLowToneBin; i <= mHighToneBin; ++i) 109 sumSignal += mRe[i]; 110 for (int i = mHighToneBin + 1; i < numSpectra; ++i) { 111 sumNoise += mRe[i]; // Note: mRe has the mag squared spectrum. 112 } 113 snr = 10.0 * log10(sumSignal / sumNoise); 114 if (snr < mDbSnrThresh) 115 (*numBadFrames) += 1; 116 } 117 continue; 118 } 119 if ((onset >= 0) && (offset > 0)) { // (3) 120 if ((maxLoc >= mLowToneBin) && (maxLoc <= mHighToneBin) && 121 (maxValue > mOnsetThresh)) { // tone should not pop up again! 122 (*numBadFrames) += 1; 123 } 124 continue; 125 } 126 } 127 if ((onset >= 0) && (offset > 0)) 128 return 1; // Success. 129 if (onset < 0) { 130 return -3; // Signal onset not found. 131 } 132 return -4; // Signal offset not found. 133 } 134 135 void GlitchTest::cleanup(void) { 136 delete [] mRe; 137 delete [] mIm; 138 delete mFt; 139 delete mWind; 140 mRe = 0; 141 mIm = 0; 142 mWind = 0; 143 mFt = 0; 144 } 145 146 int GlitchTest::realMagSqSpectrum(float* data, int numInput, 147 float* output, int* numOutput) { 148 *numOutput = 0; 149 if ((numInput <= 0) || (numInput > mFftSize)) 150 return 0; 151 int i = 0; 152 for (i = 0; i < numInput; ++i) { 153 mRe[i] = data[i]; 154 mIm[i] = 0.0; 155 } 156 for ( ; i < mFftSize; ++i) { 157 mRe[i] = 0.0; 158 mIm[i] = 0.0; 159 } 160 mFt->fft(mRe, mIm); 161 *numOutput = 1 + (mFftSize / 2); 162 for (i = 0; i < *numOutput; ++i) { 163 output[i] = (mRe[i] * mRe[i]) + (mIm[i] * mIm[i]); 164 } 165 return 1; 166 } 167 168 void GlitchTest::findPeak(float* data, int startSearch, int endSearch, 169 int* maxLoc, float* maxValue) { 170 float amax = data[startSearch]; 171 int loc = startSearch; 172 for (int i = startSearch + 1; i < endSearch; ++i) { 173 if (data[i] > amax) { 174 amax = data[i]; 175 loc = i; 176 } 177 } 178 *maxLoc = loc; 179 *maxValue = 10.0 * log10(amax / mWindowSize); 180 } 181 182
