LoopbackAnalyzer.h revision ac9d979583c4a10c2ebcf9649d773edf7acce723
1/* 2 * Copyright (C) 2017 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 17/** 18 * Tools for measuring latency and for detecting glitches. 19 * These classes are pure math and can be used with any audio system. 20 */ 21 22#ifndef AAUDIO_EXAMPLES_LOOPBACK_ANALYSER_H 23#define AAUDIO_EXAMPLES_LOOPBACK_ANALYSER_H 24 25#include <algorithm> 26#include <assert.h> 27#include <cctype> 28#include <math.h> 29#include <stdio.h> 30#include <stdlib.h> 31#include <unistd.h> 32 33#include <audio_utils/sndfile.h> 34 35// Tag for machine readable results as property = value pairs 36#define LOOPBACK_RESULT_TAG "RESULT: " 37#define LOOPBACK_SAMPLE_RATE 48000 38 39#define MILLIS_PER_SECOND 1000 40 41#define MAX_ZEROTH_PARTIAL_BINS 40 42constexpr double MAX_ECHO_GAIN = 10.0; // based on experiments, otherwise autocorrelation too noisy 43 44static const float s_Impulse[] = { 45 0.0f, 0.0f, 0.0f, 0.0f, 0.2f, // silence on each side of the impulse 46 0.5f, 0.9999f, 0.0f, -0.9999, -0.5f, // bipolar 47 -0.2f, 0.0f, 0.0f, 0.0f, 0.0f 48}; 49 50class PseudoRandom { 51public: 52 PseudoRandom() {} 53 PseudoRandom(int64_t seed) 54 : mSeed(seed) 55 {} 56 57 /** 58 * Returns the next random double from -1.0 to 1.0 59 * 60 * @return value from -1.0 to 1.0 61 */ 62 double nextRandomDouble() { 63 return nextRandomInteger() * (0.5 / (((int32_t)1) << 30)); 64 } 65 66 /** Calculate random 32 bit number using linear-congruential method. */ 67 int32_t nextRandomInteger() { 68 // Use values for 64-bit sequence from MMIX by Donald Knuth. 69 mSeed = (mSeed * (int64_t)6364136223846793005) + (int64_t)1442695040888963407; 70 return (int32_t) (mSeed >> 32); // The higher bits have a longer sequence. 71 } 72 73private: 74 int64_t mSeed = 99887766; 75}; 76 77static double calculateCorrelation(const float *a, 78 const float *b, 79 int windowSize) 80{ 81 double correlation = 0.0; 82 double sumProducts = 0.0; 83 double sumSquares = 0.0; 84 85 // Correlate a against b. 86 for (int i = 0; i < windowSize; i++) { 87 float s1 = a[i]; 88 float s2 = b[i]; 89 // Use a normalized cross-correlation. 90 sumProducts += s1 * s2; 91 sumSquares += ((s1 * s1) + (s2 * s2)); 92 } 93 94 if (sumSquares >= 0.00000001) { 95 correlation = (float) (2.0 * sumProducts / sumSquares); 96 } 97 return correlation; 98} 99 100static int calculateCorrelations(const float *haystack, int haystackSize, 101 const float *needle, int needleSize, 102 float *results, int resultSize) 103{ 104 int maxCorrelations = haystackSize - needleSize; 105 int numCorrelations = std::min(maxCorrelations, resultSize); 106 107 for (int ic = 0; ic < numCorrelations; ic++) { 108 double correlation = calculateCorrelation(&haystack[ic], needle, needleSize); 109 results[ic] = correlation; 110 } 111 112 return numCorrelations; 113} 114 115/*==========================================================================================*/ 116/** 117 * Scan until we get a correlation of a single scan that goes over the tolerance level, 118 * peaks then drops back down. 119 */ 120static double findFirstMatch(const float *haystack, int haystackSize, 121 const float *needle, int needleSize, double threshold ) 122{ 123 int ic; 124 // How many correlations can we calculate? 125 int numCorrelations = haystackSize - needleSize; 126 double maxCorrelation = 0.0; 127 int peakIndex = -1; 128 double location = -1.0; 129 const double backThresholdScaler = 0.5; 130 131 for (ic = 0; ic < numCorrelations; ic++) { 132 double correlation = calculateCorrelation(&haystack[ic], needle, needleSize); 133 134 if( (correlation > maxCorrelation) ) { 135 maxCorrelation = correlation; 136 peakIndex = ic; 137 } 138 139 //printf("PaQa_FindFirstMatch: ic = %4d, correlation = %8f, maxSum = %8f\n", 140 // ic, correlation, maxSum ); 141 // Are we past what we were looking for? 142 if((maxCorrelation > threshold) && (correlation < backThresholdScaler * maxCorrelation)) { 143 location = peakIndex; 144 break; 145 } 146 } 147 148 return location; 149} 150 151typedef struct LatencyReport_s { 152 double latencyInFrames; 153 double confidence; 154} LatencyReport; 155 156// Apply a technique similar to Harmonic Product Spectrum Analysis to find echo fundamental. 157// Using first echo instead of the original impulse for a better match. 158static int measureLatencyFromEchos(const float *haystack, int haystackSize, 159 const float *needle, int needleSize, 160 LatencyReport *report) { 161 const double threshold = 0.1; 162 printf("measureLatencyFromEchos: haystackSize = %d, needleSize = %d\n", 163 haystackSize, needleSize); 164 165 // Find first peak 166 int first = (int) (findFirstMatch(haystack, 167 haystackSize, 168 needle, 169 needleSize, 170 threshold) + 0.5); 171 172 // Use first echo as the needle for the other echos because 173 // it will be more similar. 174 needle = &haystack[first]; 175 int again = (int) (findFirstMatch(haystack, 176 haystackSize, 177 needle, 178 needleSize, 179 threshold) + 0.5); 180 181 printf("measureLatencyFromEchos: first = %d, again at %d\n", first, again); 182 first = again; 183 184 // Allocate results array 185 int remaining = haystackSize - first; 186 const int maxReasonableLatencyFrames = 48000 * 2; // arbitrary but generous value 187 int numCorrelations = std::min(remaining, maxReasonableLatencyFrames); 188 float *correlations = new float[numCorrelations]; 189 float *harmonicSums = new float[numCorrelations](); // set to zero 190 191 // Generate correlation for every position. 192 numCorrelations = calculateCorrelations(&haystack[first], remaining, 193 needle, needleSize, 194 correlations, numCorrelations); 195 196 // Add higher harmonics mapped onto lower harmonics. 197 // This reinforces the "fundamental" echo. 198 const int numEchoes = 10; 199 for (int partial = 1; partial < numEchoes; partial++) { 200 for (int i = 0; i < numCorrelations; i++) { 201 harmonicSums[i / partial] += correlations[i] / partial; 202 } 203 } 204 205 // Find highest peak in correlation array. 206 float maxCorrelation = 0.0; 207 float sumOfPeaks = 0.0; 208 int peakIndex = 0; 209 const int skip = MAX_ZEROTH_PARTIAL_BINS; // skip low bins 210 for (int i = skip; i < numCorrelations; i++) { 211 if (harmonicSums[i] > maxCorrelation) { 212 maxCorrelation = harmonicSums[i]; 213 sumOfPeaks += maxCorrelation; 214 peakIndex = i; 215 printf("maxCorrelation = %f at %d\n", maxCorrelation, peakIndex); 216 } 217 } 218 219 report->latencyInFrames = peakIndex; 220 if (sumOfPeaks < 0.0001) { 221 report->confidence = 0.0; 222 } else { 223 report->confidence = maxCorrelation / sumOfPeaks; 224 } 225 226 delete[] correlations; 227 delete[] harmonicSums; 228 return 0; 229} 230 231class AudioRecording 232{ 233public: 234 AudioRecording() { 235 } 236 ~AudioRecording() { 237 delete[] mData; 238 } 239 240 void allocate(int maxFrames) { 241 delete[] mData; 242 mData = new float[maxFrames]; 243 mMaxFrames = maxFrames; 244 } 245 246 // Write SHORT data from the first channel. 247 int write(int16_t *inputData, int inputChannelCount, int numFrames) { 248 // stop at end of buffer 249 if ((mFrameCounter + numFrames) > mMaxFrames) { 250 numFrames = mMaxFrames - mFrameCounter; 251 } 252 for (int i = 0; i < numFrames; i++) { 253 mData[mFrameCounter++] = inputData[i * inputChannelCount] * (1.0f / 32768); 254 } 255 return numFrames; 256 } 257 258 // Write FLOAT data from the first channel. 259 int write(float *inputData, int inputChannelCount, int numFrames) { 260 // stop at end of buffer 261 if ((mFrameCounter + numFrames) > mMaxFrames) { 262 numFrames = mMaxFrames - mFrameCounter; 263 } 264 for (int i = 0; i < numFrames; i++) { 265 mData[mFrameCounter++] = inputData[i * inputChannelCount]; 266 } 267 return numFrames; 268 } 269 270 int size() { 271 return mFrameCounter; 272 } 273 274 float *getData() { 275 return mData; 276 } 277 278 void setSampleRate(int32_t sampleRate) { 279 mSampleRate = sampleRate; 280 } 281 282 int32_t getSampleRate() { 283 return mSampleRate; 284 } 285 286 int save(const char *fileName, bool writeShorts = true) { 287 SNDFILE *sndFile = nullptr; 288 int written = 0; 289 SF_INFO info = { 290 .frames = mFrameCounter, 291 .samplerate = mSampleRate, 292 .channels = 1, 293 .format = SF_FORMAT_WAV | (writeShorts ? SF_FORMAT_PCM_16 : SF_FORMAT_FLOAT) 294 }; 295 296 sndFile = sf_open(fileName, SFM_WRITE, &info); 297 if (sndFile == nullptr) { 298 printf("AudioRecording::save(%s) failed to open file\n", fileName); 299 return -errno; 300 } 301 302 written = sf_writef_float(sndFile, mData, mFrameCounter); 303 304 sf_close(sndFile); 305 return written; 306 } 307 308 int load(const char *fileName) { 309 SNDFILE *sndFile = nullptr; 310 SF_INFO info; 311 312 sndFile = sf_open(fileName, SFM_READ, &info); 313 if (sndFile == nullptr) { 314 printf("AudioRecording::load(%s) failed to open file\n", fileName); 315 return -errno; 316 } 317 318 assert(info.channels == 1); 319 320 allocate(info.frames); 321 mFrameCounter = sf_readf_float(sndFile, mData, info.frames); 322 323 sf_close(sndFile); 324 return mFrameCounter; 325 } 326 327private: 328 float *mData = nullptr; 329 int32_t mFrameCounter = 0; 330 int32_t mMaxFrames = 0; 331 int32_t mSampleRate = 48000; // common default 332}; 333 334// ==================================================================================== 335class LoopbackProcessor { 336public: 337 virtual ~LoopbackProcessor() = default; 338 339 340 virtual void reset() {} 341 342 virtual void process(float *inputData, int inputChannelCount, 343 float *outputData, int outputChannelCount, 344 int numFrames) = 0; 345 346 347 virtual void report() = 0; 348 349 virtual void printStatus() {}; 350 351 virtual int getResult() { 352 return -1; 353 } 354 355 virtual bool isDone() { 356 return false; 357 } 358 359 virtual int save(const char *fileName) { 360 (void) fileName; 361 return AAUDIO_ERROR_UNIMPLEMENTED; 362 } 363 364 virtual int load(const char *fileName) { 365 (void) fileName; 366 return AAUDIO_ERROR_UNIMPLEMENTED; 367 } 368 369 virtual void setSampleRate(int32_t sampleRate) { 370 mSampleRate = sampleRate; 371 } 372 373 int32_t getSampleRate() { 374 return mSampleRate; 375 } 376 377 // Measure peak amplitude of buffer. 378 static float measurePeakAmplitude(float *inputData, int inputChannelCount, int numFrames) { 379 float peak = 0.0f; 380 for (int i = 0; i < numFrames; i++) { 381 float pos = fabs(*inputData); 382 if (pos > peak) { 383 peak = pos; 384 } 385 inputData += inputChannelCount; 386 } 387 return peak; 388 } 389 390 391private: 392 int32_t mSampleRate = LOOPBACK_SAMPLE_RATE; 393}; 394 395class PeakDetector { 396public: 397 float process(float input) { 398 float output = mPrevious * mDecay; 399 if (input > output) { 400 output = input; 401 } 402 mPrevious = output; 403 return output; 404 } 405 406private: 407 float mDecay = 0.99f; 408 float mPrevious = 0.0f; 409}; 410 411 412static void printAudioScope(float sample) { 413 const int maxStars = 80; // arbitrary, fits on one line 414 char c = '*'; 415 if (sample < -1.0) { 416 sample = -1.0; 417 c = '$'; 418 } else if (sample > 1.0) { 419 sample = 1.0; 420 c = '$'; 421 } 422 int numSpaces = (int) (((sample + 1.0) * 0.5) * maxStars); 423 for (int i = 0; i < numSpaces; i++) { 424 putchar(' '); 425 } 426 printf("%c\n", c); 427} 428 429// ==================================================================================== 430/** 431 * Measure latency given a loopback stream data. 432 * Uses a state machine to cycle through various stages including: 433 * 434 */ 435class EchoAnalyzer : public LoopbackProcessor { 436public: 437 438 EchoAnalyzer() : LoopbackProcessor() { 439 mAudioRecording.allocate(2 * getSampleRate()); 440 mAudioRecording.setSampleRate(getSampleRate()); 441 } 442 443 void setSampleRate(int32_t sampleRate) override { 444 LoopbackProcessor::setSampleRate(sampleRate); 445 mAudioRecording.setSampleRate(sampleRate); 446 } 447 448 void reset() override { 449 mDownCounter = 200; 450 mLoopCounter = 0; 451 mMeasuredLoopGain = 0.0f; 452 mEchoGain = 1.0f; 453 mState = STATE_INITIAL_SILENCE; 454 } 455 456 virtual int getResult() { 457 return mState == STATE_DONE ? 0 : -1; 458 } 459 460 virtual bool isDone() { 461 return mState == STATE_DONE || mState == STATE_FAILED; 462 } 463 464 void setGain(float gain) { 465 mEchoGain = gain; 466 } 467 468 float getGain() { 469 return mEchoGain; 470 } 471 472 void report() override { 473 474 printf("EchoAnalyzer ---------------\n"); 475 printf(LOOPBACK_RESULT_TAG "measured.gain = %f\n", mMeasuredLoopGain); 476 printf(LOOPBACK_RESULT_TAG "echo.gain = %f\n", mEchoGain); 477 printf(LOOPBACK_RESULT_TAG "test.state = %d\n", mState); 478 if (mMeasuredLoopGain >= 0.9999) { 479 printf(" ERROR - clipping, turn down volume slightly\n"); 480 } else { 481 const float *needle = s_Impulse; 482 int needleSize = (int) (sizeof(s_Impulse) / sizeof(float)); 483 float *haystack = mAudioRecording.getData(); 484 int haystackSize = mAudioRecording.size(); 485 measureLatencyFromEchos(haystack, haystackSize, needle, needleSize, &mLatencyReport); 486 if (mLatencyReport.confidence < 0.01) { 487 printf(" ERROR - confidence too low = %f\n", mLatencyReport.confidence); 488 } else { 489 double latencyMillis = 1000.0 * mLatencyReport.latencyInFrames / getSampleRate(); 490 printf(LOOPBACK_RESULT_TAG "latency.frames = %8.2f\n", mLatencyReport.latencyInFrames); 491 printf(LOOPBACK_RESULT_TAG "latency.msec = %8.2f\n", latencyMillis); 492 printf(LOOPBACK_RESULT_TAG "latency.confidence = %8.6f\n", mLatencyReport.confidence); 493 } 494 } 495 } 496 497 void printStatus() override { 498 printf("state = %d, echo gain = %f ", mState, mEchoGain); 499 } 500 501 static void sendImpulse(float *outputData, int outputChannelCount) { 502 for (float sample : s_Impulse) { 503 *outputData = sample; 504 outputData += outputChannelCount; 505 } 506 } 507 508 void process(float *inputData, int inputChannelCount, 509 float *outputData, int outputChannelCount, 510 int numFrames) override { 511 int channelsValid = std::min(inputChannelCount, outputChannelCount); 512 float peak = 0.0f; 513 int numWritten; 514 int numSamples; 515 516 echo_state_t nextState = mState; 517 518 switch (mState) { 519 case STATE_INITIAL_SILENCE: 520 // Output silence at the beginning. 521 numSamples = numFrames * outputChannelCount; 522 for (int i = 0; i < numSamples; i++) { 523 outputData[i] = 0; 524 } 525 if (mDownCounter-- <= 0) { 526 nextState = STATE_MEASURING_GAIN; 527 //printf("%5d: switch to STATE_MEASURING_GAIN\n", mLoopCounter); 528 mDownCounter = 8; 529 } 530 break; 531 532 case STATE_MEASURING_GAIN: 533 sendImpulse(outputData, outputChannelCount); 534 peak = measurePeakAmplitude(inputData, inputChannelCount, numFrames); 535 // If we get several in a row then go to next state. 536 if (peak > mPulseThreshold) { 537 if (mDownCounter-- <= 0) { 538 //printf("%5d: switch to STATE_WAITING_FOR_SILENCE, measured peak = %f\n", 539 // mLoopCounter, peak); 540 mDownCounter = 8; 541 mMeasuredLoopGain = peak; // assumes original pulse amplitude is one 542 // Calculate gain that will give us a nice decaying echo. 543 mEchoGain = mDesiredEchoGain / mMeasuredLoopGain; 544 if (mEchoGain > MAX_ECHO_GAIN) { 545 printf("ERROR - loop gain too low. Increase the volume.\n"); 546 nextState = STATE_FAILED; 547 } else { 548 nextState = STATE_WAITING_FOR_SILENCE; 549 } 550 } 551 } else { 552 mDownCounter = 8; 553 } 554 break; 555 556 case STATE_WAITING_FOR_SILENCE: 557 // Output silence and wait for the echos to die down. 558 numSamples = numFrames * outputChannelCount; 559 for (int i = 0; i < numSamples; i++) { 560 outputData[i] = 0; 561 } 562 peak = measurePeakAmplitude(inputData, inputChannelCount, numFrames); 563 // If we get several in a row then go to next state. 564 if (peak < mSilenceThreshold) { 565 if (mDownCounter-- <= 0) { 566 nextState = STATE_SENDING_PULSE; 567 //printf("%5d: switch to STATE_SENDING_PULSE\n", mLoopCounter); 568 mDownCounter = 8; 569 } 570 } else { 571 mDownCounter = 8; 572 } 573 break; 574 575 case STATE_SENDING_PULSE: 576 mAudioRecording.write(inputData, inputChannelCount, numFrames); 577 sendImpulse(outputData, outputChannelCount); 578 nextState = STATE_GATHERING_ECHOS; 579 //printf("%5d: switch to STATE_GATHERING_ECHOS\n", mLoopCounter); 580 break; 581 582 case STATE_GATHERING_ECHOS: 583 numWritten = mAudioRecording.write(inputData, inputChannelCount, numFrames); 584 peak = measurePeakAmplitude(inputData, inputChannelCount, numFrames); 585 if (peak > mMeasuredLoopGain) { 586 mMeasuredLoopGain = peak; // AGC might be raising gain so adjust it on the fly. 587 // Recalculate gain that will give us a nice decaying echo. 588 mEchoGain = mDesiredEchoGain / mMeasuredLoopGain; 589 } 590 // Echo input to output. 591 for (int i = 0; i < numFrames; i++) { 592 int ic; 593 for (ic = 0; ic < channelsValid; ic++) { 594 outputData[ic] = inputData[ic] * mEchoGain; 595 } 596 for (; ic < outputChannelCount; ic++) { 597 outputData[ic] = 0; 598 } 599 inputData += inputChannelCount; 600 outputData += outputChannelCount; 601 } 602 if (numWritten < numFrames) { 603 nextState = STATE_DONE; 604 //printf("%5d: switch to STATE_DONE\n", mLoopCounter); 605 } 606 break; 607 608 case STATE_DONE: 609 default: 610 break; 611 } 612 613 mState = nextState; 614 mLoopCounter++; 615 } 616 617 int save(const char *fileName) override { 618 return mAudioRecording.save(fileName); 619 } 620 621 int load(const char *fileName) override { 622 return mAudioRecording.load(fileName); 623 } 624 625private: 626 627 enum echo_state_t { 628 STATE_INITIAL_SILENCE, 629 STATE_MEASURING_GAIN, 630 STATE_WAITING_FOR_SILENCE, 631 STATE_SENDING_PULSE, 632 STATE_GATHERING_ECHOS, 633 STATE_DONE, 634 STATE_FAILED 635 }; 636 637 int mDownCounter = 500; 638 int mLoopCounter = 0; 639 float mPulseThreshold = 0.02f; 640 float mSilenceThreshold = 0.002f; 641 float mMeasuredLoopGain = 0.0f; 642 float mDesiredEchoGain = 0.95f; 643 float mEchoGain = 1.0f; 644 echo_state_t mState = STATE_INITIAL_SILENCE; 645 646 AudioRecording mAudioRecording; // contains only the input after the gain detection burst 647 LatencyReport mLatencyReport; 648 // PeakDetector mPeakDetector; 649}; 650 651 652// ==================================================================================== 653/** 654 * Output a steady sinewave and analyze the return signal. 655 * 656 * Use a cosine transform to measure the predicted magnitude and relative phase of the 657 * looped back sine wave. Then generate a predicted signal and compare with the actual signal. 658 */ 659class SineAnalyzer : public LoopbackProcessor { 660public: 661 662 virtual int getResult() { 663 return mState == STATE_LOCKED ? 0 : -1; 664 } 665 666 void report() override { 667 printf("SineAnalyzer ------------------\n"); 668 printf(LOOPBACK_RESULT_TAG "peak.amplitude = %7.5f\n", mPeakAmplitude); 669 printf(LOOPBACK_RESULT_TAG "sine.magnitude = %7.5f\n", mMagnitude); 670 printf(LOOPBACK_RESULT_TAG "phase.offset = %7.5f\n", mPhaseOffset); 671 printf(LOOPBACK_RESULT_TAG "ref.phase = %7.5f\n", mPhase); 672 printf(LOOPBACK_RESULT_TAG "frames.accumulated = %6d\n", mFramesAccumulated); 673 printf(LOOPBACK_RESULT_TAG "sine.period = %6d\n", mPeriod); 674 printf(LOOPBACK_RESULT_TAG "test.state = %6d\n", mState); 675 printf(LOOPBACK_RESULT_TAG "frame.count = %6d\n", mFrameCounter); 676 // Did we ever get a lock? 677 bool gotLock = (mState == STATE_LOCKED) || (mGlitchCount > 0); 678 if (!gotLock) { 679 printf("ERROR - failed to lock on reference sine tone\n"); 680 } else { 681 // Only print if meaningful. 682 printf(LOOPBACK_RESULT_TAG "glitch.count = %6d\n", mGlitchCount); 683 } 684 } 685 686 void printStatus() override { 687 printf(" state = %d, glitches = %3d,", mState, mGlitchCount); 688 } 689 690 double calculateMagnitude(double *phasePtr = NULL) { 691 if (mFramesAccumulated == 0) { 692 return 0.0; 693 } 694 double sinMean = mSinAccumulator / mFramesAccumulated; 695 double cosMean = mCosAccumulator / mFramesAccumulated; 696 double magnitude = 2.0 * sqrt( (sinMean * sinMean) + (cosMean * cosMean )); 697 if( phasePtr != NULL ) 698 { 699 double phase = M_PI_2 - atan2( sinMean, cosMean ); 700 *phasePtr = phase; 701 } 702 return magnitude; 703 } 704 705 /** 706 * @param inputData contains microphone data with sine signal feedback 707 * @param outputData contains the reference sine wave 708 */ 709 void process(float *inputData, int inputChannelCount, 710 float *outputData, int outputChannelCount, 711 int numFrames) override { 712 float peak = measurePeakAmplitude(inputData, inputChannelCount, numFrames); 713 if (peak > mPeakAmplitude) { 714 mPeakAmplitude = peak; 715 } 716 717 for (int i = 0; i < numFrames; i++) { 718 float sample = inputData[i * inputChannelCount]; 719 720 float sinOut = sinf(mPhase); 721 722 switch (mState) { 723 case STATE_IMMUNE: 724 case STATE_WAITING_FOR_SIGNAL: 725 break; 726 case STATE_WAITING_FOR_LOCK: 727 mSinAccumulator += sample * sinOut; 728 mCosAccumulator += sample * cosf(mPhase); 729 mFramesAccumulated++; 730 // Must be a multiple of the period or the calculation will not be accurate. 731 if (mFramesAccumulated == mPeriod * 4) { 732 mPhaseOffset = 0.0; 733 mMagnitude = calculateMagnitude(&mPhaseOffset); 734 if (mMagnitude > mThreshold) { 735 if (fabs(mPreviousPhaseOffset - mPhaseOffset) < 0.001) { 736 mState = STATE_LOCKED; 737 //printf("%5d: switch to STATE_LOCKED\n", mFrameCounter); 738 } 739 mPreviousPhaseOffset = mPhaseOffset; 740 } 741 resetAccumulator(); 742 } 743 break; 744 745 case STATE_LOCKED: { 746 // Predict next sine value 747 float predicted = sinf(mPhase + mPhaseOffset) * mMagnitude; 748 // printf(" predicted = %f, actual = %f\n", predicted, sample); 749 750 float diff = predicted - sample; 751 if (fabs(diff) > mTolerance) { 752 mGlitchCount++; 753 //printf("%5d: Got a glitch # %d, predicted = %f, actual = %f\n", 754 // mFrameCounter, mGlitchCount, predicted, sample); 755 mState = STATE_IMMUNE; 756 //printf("%5d: switch to STATE_IMMUNE\n", mFrameCounter); 757 mDownCounter = mPeriod; // Set duration of IMMUNE state. 758 } 759 } break; 760 } 761 762 // Output sine wave so we can measure it. 763 outputData[i * outputChannelCount] = (sinOut * mOutputAmplitude) 764 + (mWhiteNoise.nextRandomDouble() * mNoiseAmplitude); 765 // printf("%5d: sin(%f) = %f, %f\n", i, mPhase, sinOut, mPhaseIncrement); 766 767 // advance and wrap phase 768 mPhase += mPhaseIncrement; 769 if (mPhase > M_PI) { 770 mPhase -= (2.0 * M_PI); 771 } 772 773 mFrameCounter++; 774 } 775 776 // Do these once per buffer. 777 switch (mState) { 778 case STATE_IMMUNE: 779 mDownCounter -= numFrames; 780 if (mDownCounter <= 0) { 781 mState = STATE_WAITING_FOR_SIGNAL; 782 //printf("%5d: switch to STATE_WAITING_FOR_SIGNAL\n", mFrameCounter); 783 } 784 break; 785 case STATE_WAITING_FOR_SIGNAL: 786 if (peak > mThreshold) { 787 mState = STATE_WAITING_FOR_LOCK; 788 //printf("%5d: switch to STATE_WAITING_FOR_LOCK\n", mFrameCounter); 789 resetAccumulator(); 790 } 791 break; 792 case STATE_WAITING_FOR_LOCK: 793 case STATE_LOCKED: 794 break; 795 } 796 797 } 798 799 void resetAccumulator() { 800 mFramesAccumulated = 0; 801 mSinAccumulator = 0.0; 802 mCosAccumulator = 0.0; 803 } 804 805 void reset() override { 806 mGlitchCount = 0; 807 mState = STATE_IMMUNE; 808 mPhaseIncrement = 2.0 * M_PI / mPeriod; 809 printf("phaseInc = %f for period %d\n", mPhaseIncrement, mPeriod); 810 resetAccumulator(); 811 } 812 813private: 814 815 enum sine_state_t { 816 STATE_IMMUNE, 817 STATE_WAITING_FOR_SIGNAL, 818 STATE_WAITING_FOR_LOCK, 819 STATE_LOCKED 820 }; 821 822 int mPeriod = 79; 823 double mPhaseIncrement = 0.0; 824 double mPhase = 0.0; 825 double mPhaseOffset = 0.0; 826 double mPreviousPhaseOffset = 0.0; 827 double mMagnitude = 0.0; 828 double mThreshold = 0.005; 829 double mTolerance = 0.01; 830 int32_t mFramesAccumulated = 0; 831 double mSinAccumulator = 0.0; 832 double mCosAccumulator = 0.0; 833 int32_t mGlitchCount = 0; 834 double mPeakAmplitude = 0.0; 835 int mDownCounter = 4000; 836 int32_t mFrameCounter = 0; 837 float mOutputAmplitude = 0.75; 838 839 PseudoRandom mWhiteNoise; 840 float mNoiseAmplitude = 0.00; // Used to experiment with warbling caused by DRC. 841 842 sine_state_t mState = STATE_IMMUNE; 843}; 844 845 846#undef LOOPBACK_SAMPLE_RATE 847#undef LOOPBACK_RESULT_TAG 848 849#endif /* AAUDIO_EXAMPLES_LOOPBACK_ANALYSER_H */ 850