1 /*
2  * Copyright (C) 2022 The Android Open Source Project
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of 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,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 package com.android.cts.verifier.audio.analyzers;
17 
18 import java.util.Random;
19 
20 /**
21  * Output a steady sine wave and analyze the return signal.
22  *
23  * Use a cosine transform to measure the predicted magnitude and relative phase of the
24  * looped back sine wave. Then generate a predicted signal and compare with the actual signal.
25  *
26  * Derived from oboetester::BaseSineAnalyzer
27  */
28 public class BaseSineAnalyzer implements SignalAnalyzer {
29     @SuppressWarnings("unused")
30     static final String TAG = "BaseSineAnalyzer";
31 
32     int  mSinePeriod = 1; // this will be set before use
33     double  mInverseSinePeriod = 1.0;
34     double  mPhaseIncrement = 0.0;
35     double  mOutputPhase = 0.0;
36     double  mOutputAmplitude = 0.75;
37     double  mPreviousPhaseOffset = 0.0;
38     double  mPhaseTolerance = 2 * Math.PI  / 48;
39 
40     double mMagnitude = 0.0;
41     double mMaxMagnitude = 0.0;
42 
43     double mPhaseErrorSum;
44     double mPhaseErrorCount;
45     double mPhaseJitter = 0.0;
46 
47     int mPhaseCount = 0;
48 
49     // If this jumps around then we are probably just hearing noise.
50     double  mPhaseOffset = 0.0;
51     int mFramesAccumulated = 0;
52     double  mSinAccumulator = 0.0;
53     double  mCosAccumulator = 0.0;
54     double  mScaledTolerance = 0.0;
55     double  mTolerance = 0.10; // scaled from 0.0 to 1.0
56     int mInputChannel = 0;
57     int mOutputChannel = 0;
58 
59     static final int DEFAULT_SAMPLERATE = 48000;
60     static final int MILLIS_PER_SECOND = 1000;  // by definition
61     static final int MAX_LATENCY_MILLIS = 1000;  // arbitrary and generous
62     static final int TARGET_GLITCH_FREQUENCY = 1000;
63     static final double MIN_REQUIRED_MAGNITUDE = 0.001;
64     static final int TYPICAL_SAMPLE_RATE = 48000;
65     static final double MAX_SINE_FREQUENCY = 1000.0;
66     static final double FRAMES_PER_CYCLE = TYPICAL_SAMPLE_RATE / MAX_SINE_FREQUENCY;
67     static final double PHASE_PER_BIN = 2.0 * Math.PI / FRAMES_PER_CYCLE;
68     static final double MAX_ALLOWED_JITTER = 2.0 * PHASE_PER_BIN;
69     // Start by failing then let good results drive us into a pass value.
70     static final double INITIAL_JITTER = 2.0 * MAX_ALLOWED_JITTER;
71     // A coefficient of 0.0 is no filtering. 0.9999 is extreme low pass.
72     static final double JITTER_FILTER_COEFFICIENT = 0.8;
73 
74     int mSampleRate = DEFAULT_SAMPLERATE;
75 
76     double mNoiseAmplitude = 0.00; // Used to experiment with warbling caused by DRC.
77     Random mWhiteNoise = new Random();
78 
79     MagnitudePhase mMagPhase = new MagnitudePhase();
80 
81     InfiniteRecording mInfiniteRecording = new InfiniteRecording(64 * 1024);
82 
83     enum RESULT_CODE {
84         RESULT_OK,
85         ERROR_NOISY,
86         ERROR_VOLUME_TOO_LOW,
87         ERROR_VOLUME_TOO_HIGH,
88         ERROR_CONFIDENCE,
89         ERROR_INVALID_STATE,
90         ERROR_GLITCHES,
91         ERROR_NO_LOCK
92     };
93 
BaseSineAnalyzer()94     public BaseSineAnalyzer() {
95         // Add a little bit of noise to reduce blockage by speaker protection and DRC.
96         mNoiseAmplitude = 0.02;
97     };
98 
getSampleRate()99     public int getSampleRate() {
100         return mSampleRate;
101     }
102 
103     /**
104      * Set the assumed sample rate for the analysis
105      * @param sampleRate
106      */
setSampleRate(int sampleRate)107     public void setSampleRate(int sampleRate) {
108         mSampleRate = sampleRate;
109         updatePhaseIncrement();
110     }
111 
112     /**
113      * @return output frequency that will have an integer period on input
114      */
getAdjustedFrequency()115     public double getAdjustedFrequency() {
116         updatePhaseIncrement();
117         return mInverseSinePeriod * getSampleRate();
118     }
119 
setInputChannel(int inputChannel)120     public void setInputChannel(int inputChannel) {
121         mInputChannel = inputChannel;
122     }
123 
getInputChannel()124     public int getInputChannel() {
125         return mInputChannel;
126     }
127 
setOutputChannel(int outputChannel)128     public void setOutputChannel(int outputChannel) {
129         mOutputChannel = outputChannel;
130     }
131 
getOutputChannel()132     public int getOutputChannel() {
133         return mOutputChannel;
134     }
135 
setNoiseAmplitude(double noiseAmplitude)136     public void setNoiseAmplitude(double noiseAmplitude) {
137         mNoiseAmplitude = noiseAmplitude;
138     }
139 
getNoiseAmplitude()140     public double getNoiseAmplitude() {
141         return mNoiseAmplitude;
142     }
143 
setMagnitude(double magnitude)144     void setMagnitude(double magnitude) {
145         mMagnitude = magnitude;
146         mScaledTolerance = mMagnitude * mTolerance;
147     }
148 
getTolerance()149     public double getTolerance() {
150         return mTolerance;
151     }
152 
setTolerance(double tolerance)153     public void setTolerance(double tolerance) {
154         mTolerance = tolerance;
155     }
156 
getMagnitude()157     public double getMagnitude() {
158         return mMagnitude;
159     }
160 
getMaxMagnitude()161     public double getMaxMagnitude() {
162         return mMaxMagnitude;
163     }
164 
getPhaseOffset()165     public double getPhaseOffset() {
166         return mPhaseOffset;
167     }
168 
getOutputPhase()169     public double getOutputPhase() {
170         return mOutputPhase;
171     }
172 
getPhaseJitter()173     public double getPhaseJitter() {
174         return mPhaseJitter;
175     }
176 
177     // reset the sine wave detector
resetAccumulator()178     void resetAccumulator() {
179         mFramesAccumulated = 0;
180         mSinAccumulator = 0.0;
181         mCosAccumulator = 0.0;
182     }
183 
184     class MagnitudePhase {
185         public double mMagnitude;
186         public double mPhase;
187     }
188 
189     /**
190      * Calculate the magnitude of the component of the input signal
191      * that matches the analysis frequency.
192      * Also calculate the phase that we can use to create a
193      * signal that matches that component.
194      * The phase will be between -PI and +PI.
195      */
calculateMagnitudePhase(MagnitudePhase magphase)196     double calculateMagnitudePhase(MagnitudePhase magphase) {
197         if (mFramesAccumulated == 0) {
198             return 0.0;
199         }
200         double sinMean = mSinAccumulator / mFramesAccumulated;
201         double cosMean = mCosAccumulator / mFramesAccumulated;
202 
203         double magnitude = 2.0 * Math.sqrt((sinMean * sinMean) + (cosMean * cosMean));
204         magphase.mPhase = Math.atan2(cosMean, sinMean);
205         return magphase.mMagnitude = magnitude;
206     }
207 
208     // advance and wrap phase
incrementOutputPhase()209     void incrementOutputPhase() {
210         mOutputPhase += mPhaseIncrement;
211         if (mOutputPhase > Math.PI) {
212             mOutputPhase -= (2.0 * Math.PI);
213         }
214     }
215 
calculatePhaseError(double p1, double p2)216     double calculatePhaseError(double p1, double p2) {
217         double diff = p1 - p2;
218         // Wrap around the circle.
219         while (diff > Math.PI) {
220             diff -= 2 * Math.PI;
221         }
222         while (diff < -Math.PI) {
223             diff += 2 * Math.PI;
224         }
225         return diff;
226     }
227 
getAveragePhaseError()228     double getAveragePhaseError() {
229         // If we have no measurements then return maximum possible phase jitter
230         // to avoid dividing by zero.
231         return (mPhaseErrorCount > 0) ? (mPhaseErrorSum / mPhaseErrorCount) : Math.PI;
232     }
233 
234     /**
235      * Perform sin/cos analysis on each sample.
236      * Measure magnitude and phase on every period.
237      * @param sample
238      * @param referencePhase
239      * @return true if magnitude and phase updated
240      */
transformSample(float sample, double referencePhase)241     boolean transformSample(float sample, double referencePhase) {
242         // Track incoming signal and slowly adjust magnitude to account
243         // for drift in the DRC or AGC.
244         mSinAccumulator += ((double) sample) * Math.sin(referencePhase);
245         mCosAccumulator += ((double) sample) * Math.cos(referencePhase);
246         mFramesAccumulated++;
247 
248         incrementOutputPhase();
249 
250         // Must be a multiple of the period or the calculation will not be accurate.
251         if (mFramesAccumulated == mSinePeriod) {
252             final double coefficient = 0.1;
253 
254             double magnitude = calculateMagnitudePhase(mMagPhase);
255             mPhaseOffset = mMagPhase.mPhase;
256             // One pole averaging filter.
257             setMagnitude((mMagnitude * (1.0 - coefficient)) + (magnitude * coefficient));
258             return true;
259         } else {
260             return false;
261         }
262     }
263 
264     /**
265      * @param frameData contains microphone data with sine signal feedback
266      * @param channelCount
267      */
processInputFrame(float[] frameData, int offset)268     RESULT_CODE processInputFrame(float[] frameData, int offset) {
269         RESULT_CODE result = RESULT_CODE.RESULT_OK;
270 
271         float sample = frameData[offset];
272         mInfiniteRecording.write(sample);
273 
274         if (transformSample(sample, mOutputPhase)) {
275             resetAccumulator();
276             if (mMagnitude >= MIN_REQUIRED_MAGNITUDE) {
277                 // Analyze magnitude and phase on every period.
278                 double phaseError =
279                         Math.abs(calculatePhaseError(mPhaseOffset, mPreviousPhaseOffset));
280                 if (phaseError < mPhaseTolerance) {
281                     mMaxMagnitude = Math.max(mMagnitude, mMaxMagnitude);
282                 }
283                 mPreviousPhaseOffset = mPhaseOffset;
284 
285                 // Only look at the phase if we have a signal.
286                 if (mPhaseCount > 3) {
287                     // Accumulate phase error and average.
288                     mPhaseErrorSum += phaseError;
289                     mPhaseErrorCount++;
290                     mPhaseJitter = getAveragePhaseError();
291                 }
292 
293                 mPhaseCount++;
294             }
295         }
296         return result;
297     }
298 
updatePhaseIncrement()299     private void updatePhaseIncrement() {
300         mSinePeriod = getSampleRate() / TARGET_GLITCH_FREQUENCY;
301         mInverseSinePeriod = 1.0 / mSinePeriod;
302         mPhaseIncrement = 2.0 * Math.PI * mInverseSinePeriod;
303     }
304 
305     @Override
reset()306     public void reset() {
307         resetAccumulator();
308 
309         mOutputPhase = 0.0f;
310         mMagnitude = 0.0;
311         mMaxMagnitude = 0.0;
312         mPhaseOffset = 0.0;
313         mPreviousPhaseOffset = 0.0;
314         mPhaseJitter = INITIAL_JITTER;
315         mPhaseCount = 0;
316         mPhaseErrorSum = 0.0;
317         mPhaseErrorCount = 0.0;
318 
319         updatePhaseIncrement();
320     }
321 
322     @Override
analyzeBuffer(float[] audioData, int numChannels, int numFrames)323     public void analyzeBuffer(float[] audioData, int numChannels, int numFrames) {
324         int offset = mInputChannel;
325         for (int frameIndex = 0; frameIndex < numFrames; frameIndex++) {
326             // processOutputFrame(audioData, offset, numChannels);
327             processInputFrame(audioData, offset);
328             offset += numChannels;
329         }
330 
331 //        // Only look at the phase if we have a signal.
332 //        if (mMagnitude >= MIN_REQUIRED_MAGNITUDE) {
333 //            double phase = mPhaseOffset;
334 //            if (mPhaseCount > 3) {
335 //                double phaseError = calculatePhaseError(phase, mPhaseOffset);
336 //                // Accumulate phase error and average.
337 //                mPhaseErrorSum += phaseError;
338 //                mPhaseErrorCount++;
339 //                mPhaseJitter = getAveragePhaseError();
340 //            }
341 //
342 //            mPhaseOffset = phase;
343 //            mPhaseCount++;
344 //        }
345     }
346 }
347