| /* |
| * Copyright (C) 2010 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. |
| */ |
| |
| #include "GlitchTest.h" |
| #include "Window.h" |
| #include "Fft.h" |
| |
| #include <math.h> |
| |
| GlitchTest::GlitchTest(void) : mRe(0), mIm(0), mFt(0), mWind(0) {} |
| |
| void GlitchTest::init(float sampleRate, float stimFreq, float onsetThresh, |
| float dbSnrThresh) { |
| cleanup(); // in case Init gets called multiple times. |
| // Analysis parameters |
| float windowDur = 0.025; // Duration of the analysis window in seconds. |
| float frameInterval = 0.01; // Interval in seconds between frame onsets |
| float lowestFreq = 200.0; // Lowest frequency to include in |
| // background noise power integration. |
| mSampleRate = sampleRate; |
| mOnsetThresh = onsetThresh; |
| mDbSnrThresh = dbSnrThresh; |
| mFrameStep = (int)(0.5 + (mSampleRate * frameInterval)); |
| mWindowSize = (int)(0.5 + (mSampleRate * windowDur)); |
| mWind = new Window(mWindowSize); |
| mFt = new Fft; |
| mFt->fftInit(mFt->fftPow2FromWindowSize(mWindowSize)); |
| mFftSize = mFt->fftGetSize(); |
| mRe = new float[mFftSize]; |
| mIm = new float[mFftSize]; |
| float freqInterval = mSampleRate / mFftSize; |
| // We can exclude low frequencies from consideration, since the |
| // phone may have DC offset in the A/D, and there may be rumble in |
| // the testing room. |
| mLowestSpectrumBin = (int)(0.5 + (lowestFreq / freqInterval)); |
| // These are the bin indices within which most of the energy due to |
| // the (windowed) tone should be found. |
| mLowToneBin = (int)(0.5 + (stimFreq / freqInterval)) - 2; |
| mHighToneBin = (int)(0.5 + (stimFreq / freqInterval)) + 2; |
| if (mLowestSpectrumBin >= mLowToneBin) { |
| mLowestSpectrumBin = mHighToneBin + 1; |
| } |
| } |
| |
| int GlitchTest::checkToneSnr(short* pcm, int numSamples, float* duration, |
| int* numBadFrames) { |
| *numBadFrames = 0; |
| *duration = 0.0; |
| if (!(mRe && mIm)) { |
| return -1; // not initialized. |
| } |
| int n_frames = 1 + ((numSamples - mWindowSize) / mFrameStep); |
| if (n_frames < 4) { // pathologically short input signal |
| return -2; |
| } |
| *numBadFrames = 0; |
| int onset = -1; |
| int offset = -1; |
| for (int frame = 0; frame < n_frames; ++frame) { |
| int numSpectra = 0; |
| mWind->window(pcm + frame*mFrameStep, mRe, 0.0); |
| realMagSqSpectrum(mRe, mWindowSize, mRe, &numSpectra); |
| int maxLoc = 0; |
| float maxValue = 0.0; |
| findPeak(mRe, mLowestSpectrumBin, numSpectra, &maxLoc, &maxValue); |
| // possible states: (1) before tone onset; (2) within tone |
| // region; (3) after tone offset. |
| if ((onset < 0) && (offset < 0)) { // (1) |
| if ((maxLoc >= mLowToneBin) && (maxLoc <= mHighToneBin) |
| && (maxValue > mOnsetThresh)) { |
| onset = frame; |
| } |
| continue; |
| } |
| if ((onset >= 0) && (offset < 0)) { // (2) |
| if (frame > (onset + 2)) { // let the framer get completely |
| // into the tonal signal |
| double sumNoise = 1.0; // init. to small non-zero vals to |
| // avoid log or divz problems |
| double sumSignal = 1.0; |
| float snr = 0.0; |
| if (maxValue < mOnsetThresh) { |
| offset = frame; |
| *duration = mFrameStep * (offset - onset) / mSampleRate; |
| if (*numBadFrames >= 1) { |
| (*numBadFrames) -= 1; // account for expected glitch at |
| // signal offset |
| } |
| continue; |
| } |
| for (int i = mLowestSpectrumBin; i < mLowToneBin; ++i) { |
| sumNoise += mRe[i]; // note that mRe contains the magnitude |
| // squared spectrum. |
| } |
| for (int i = mLowToneBin; i <= mHighToneBin; ++i) |
| sumSignal += mRe[i]; |
| for (int i = mHighToneBin + 1; i < numSpectra; ++i) { |
| sumNoise += mRe[i]; // Note: mRe has the mag squared spectrum. |
| } |
| snr = 10.0 * log10(sumSignal / sumNoise); |
| if (snr < mDbSnrThresh) |
| (*numBadFrames) += 1; |
| } |
| continue; |
| } |
| if ((onset >= 0) && (offset > 0)) { // (3) |
| if ((maxLoc >= mLowToneBin) && (maxLoc <= mHighToneBin) && |
| (maxValue > mOnsetThresh)) { // tone should not pop up again! |
| (*numBadFrames) += 1; |
| } |
| continue; |
| } |
| } |
| if ((onset >= 0) && (offset > 0)) |
| return 1; // Success. |
| if (onset < 0) { |
| return -3; // Signal onset not found. |
| } |
| return -4; // Signal offset not found. |
| } |
| |
| void GlitchTest::cleanup(void) { |
| delete [] mRe; |
| delete [] mIm; |
| delete mFt; |
| delete mWind; |
| mRe = 0; |
| mIm = 0; |
| mWind = 0; |
| mFt = 0; |
| } |
| |
| int GlitchTest::realMagSqSpectrum(float* data, int numInput, |
| float* output, int* numOutput) { |
| *numOutput = 0; |
| if ((numInput <= 0) || (numInput > mFftSize)) |
| return 0; |
| int i = 0; |
| for (i = 0; i < numInput; ++i) { |
| mRe[i] = data[i]; |
| mIm[i] = 0.0; |
| } |
| for ( ; i < mFftSize; ++i) { |
| mRe[i] = 0.0; |
| mIm[i] = 0.0; |
| } |
| mFt->fft(mRe, mIm); |
| *numOutput = 1 + (mFftSize / 2); |
| for (i = 0; i < *numOutput; ++i) { |
| output[i] = (mRe[i] * mRe[i]) + (mIm[i] * mIm[i]); |
| } |
| return 1; |
| } |
| |
| void GlitchTest::findPeak(float* data, int startSearch, int endSearch, |
| int* maxLoc, float* maxValue) { |
| float amax = data[startSearch]; |
| int loc = startSearch; |
| for (int i = startSearch + 1; i < endSearch; ++i) { |
| if (data[i] > amax) { |
| amax = data[i]; |
| loc = i; |
| } |
| } |
| *maxLoc = loc; |
| *maxValue = 10.0 * log10(amax / mWindowSize); |
| } |
| |