AudioAnalyzePhase_F32.h

/*
 * AudioAnalyzePhase_F32.h
 *
 * 31 March 2020, Rev 8 April 2020  
 * Status Tested OK T3.6 and T4.0. 
 * Bob Larkin, in support of the library:
 * Chip Audette, OpenAudio, Apr 2017
 *     -------------------
 * There are two inputs, 0 and 1 (Left and Right)
 * There is one output, the phase angle between 0 & 1 expressed in
 * radians (180 degrees is Pi radians) or degrees.  This is a 180-degree
 * type of phase detector.  See RadioIQMixer_F32 for a 360 degree type.
 *
 * This block can be used to measure phase between two sinusoids, and the default IIR filter is suitable for this with a cut-off
 * frequency of 100 Hz. The only IIR configuration is 4-cascaded satages of BiQuad.  For this, 20 coefficients must be provided
 * in 4 times (b0, b1, b2, -a1, -a2) order (example below). This IIR filter inherently does not have very good
 * linearity in phase vs. frequency.  This can be a problem for communications systems.
 * As an alternative, a linear phase (as long as coefficients are symmetrical)
 * FIR filter can be set up with the begin method.  The built in FIR LP filter has a cutoff frequency of 4 kHz when used
 * at a 44.1 kHz sample rate.  This filter uses 53 coefficients (called taps).  Any FIR filter with 4 to 200 coefficients can be used
 * as set up by the begin method.
 *
 * DEFAULTS: 100 Hz IIR LP, output is in radians, and this does *NOT* need a call to begin().  This can be changed, including
 * using a FIR LP where linear phase is needed, or NO_LP_FILTER that leaves harmonics of the input frequency.  Method begin()
 * changes the options.  For instance, to use a 60 coefficient FIR the setup() in the .INO might do
 *    myAnalyzePhase.begin(FIR_LP_FILTER, &myFIRCoefficients[0], 60, DEGREES_PHASE);
 * If _pcoefficients is NULL, the coefficients will be left default.  For instance, to use the default 100 Hz IIR filter, with degree output
 *    myAnalyzePhase.begin(IIR_LP_FILTER, NULL, 20, DEGREES_PHASE);
 * To provide a new set of IIR coefficients (note strange coefficient order and negation for a() that CMSIS needs)
 *    myAnalyzePhase.begin(IIR_LP_FILTER, &myIIRCoefficients[0], 20, RADIANS_PHASE);
 * In begin() the pdConfig can be set (see #defines below).  The default is to use no limiter, but to measure the input levels over the 
 * block and use that to scale the multiplier output.  This will cause successive blocks to change slightly in output level due to
 * errors in level measurement, but is other wise fine.  If the limiter is used, the narrow band IIR filter should also be used to
 * prevent artifacts from "beats" between the sample rate and the input frequency.
 * 
 * Three different scaling routines are available following the LP filter.  These deal with the issue that the multiplier type
 * of phase detector produces an output proportional to the cosine of the phase angle between the two input sine waves.
 * If the inputs both have a magnitude ranging from -1.0 to 1.0, the output will be cos(phase difference).  Other values of
 * sine wave will multiply this by the product of the two maximum levels.  The selection of "fast" or "accurate" acos() will
 * make the output approximately the angle, as scaled by UNITS_MASK.  The ACOS_MASK bits in pdConfig, set by begin(), selects the
 * acos used.  Note that if acos function is used, the output range is 0 to pi radians, i.e., 0 to 180 degrees.  "Units" have no
 * effect when acos90 is not being used, as that would make little sense for the (-1,1) output.
 * 
 * Functions:
 *    setAnalyzePhaseConfig(const uint16_t LPType, float32_t *pCoeffs, uint16_t nCoeffs)
 *    setAnalyzePhaseConfig(const uint16_t LPType, float32_t *pCoeffs, uint16_t nCoeffs, uint16_t pdConfig)
 *      are used to chhange the output filter from the IIR default, where:
 *        LPType is  NO_LP_FILTER, IIR_LP_FILTER, FIR_LP_FILTER to select the output filter
 *        pCoeffs is a pointer to filter coefficients, either IIR or FIR 
 *        nCoeffs is the number of filter coefficients
 *        pdConfig is bitwise selection (default 0b1100) of
 *          Bit 0: 0=No Limiter (default)   1=Use limiter
 *          Bit 2 and 1: 00=Use no acos linearizer    01=undefined
 *                       10=Fast, math-continuous acos() (default)    11=Accurate acosf()
 *          Bit 3: 0=No scale of multiplier  1=scale to min-max (default)
 *          Bit 4: 0=Output in degrees    1=Output in radians (default)
 *    showError(uint16_t e) sets whether error printing comes from update (e=1) or not (e=0).
 * 
 * Examples:   AudioTestAnalyzePhase.ino and AudioTestSinCos.ino 
 * 
 * Some measured time data for a 128 size block, Teensy 3.6, parts of update():
 *   Default settings, total time 123 microseconds
 *   Overhead of update(), loading arrays, handling blocks, less than 2 microseconds
 *   Min-max calculation, 23 microseconds
 *   Multiplier DBMixer  8 microseconds
 *   IIR LPF (default filter)  57 microseconds
 *   53-term FIR filter  149 microseconds 
 *   Fast acos_32() linearizer  32 microseconds
 *   Accurate acosf(x) seems to vary (with x?), 150 to 350 microsecond range
 * 
 * Measured total update() time for the min-max scaling, fast acos(), and 53-term FIR filtering
 *      case is 214 microseconds for Teensy 3.6 and 45 microseconds for Teensy 4.0.
 *
 * Copyright (c) 2020 Bob Larkin
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#ifndef _analyze_phase_f32_h
#define _analyze_phase_f32_h

#define N_STAGES 4
#define NFIR_MAX 200
#define NO_LP_FILTER 0
#define IIR_LP_FILTER 1
#define FIR_LP_FILTER 2
#define RADIANS_PHASE 1.0
#define DEGREES_PHASE 57.295779

// Test the number of microseconds to execute update()
#define TEST_TIME 1

#define LIMITER_MASK 0b00001
#define ACOS_MASK    0b00110
#define SCALE_MASK   0b01000
#define UNITS_MASK   0b10000

#include "AudioStream_F32.h"
#include <math.h>

class AudioAnalyzePhase_F32 : public AudioStream_F32 {
//GUI: inputs:2, outputs:1  //this line used for automatic generation of GUI node
//GUI: shortName: AnalyzePhase
public:
    // Option of AudioSettings_F32 change to block size or sample rate:
    AudioAnalyzePhase_F32(void) :  AudioStream_F32(2, inputQueueArray_f32) { // default block_size and sampleRate_Hz
	// Initialize BiQuad IIR instance (ARM DSP Math Library)
        arm_biquad_cascade_df1_init_f32(&iir_inst, N_STAGES, &iir_coeffs[0],  &IIRStateF32[0]);
    }
    // Constructor including new block_size and/or sampleRate_Hz
    AudioAnalyzePhase_F32(const AudioSettings_F32 &settings) :  AudioStream_F32(2, inputQueueArray_f32) {
        block_size = settings.audio_block_samples;
	    sampleRate_Hz = settings.sample_rate_Hz;        
	// Initialize BiQuad IIR instance (ARM DSP Math Library)
        arm_biquad_cascade_df1_init_f32(&iir_inst, N_STAGES, &iir_coeffs[0],  &IIRStateF32[0]);
    }

    // Set AnalyzePhaseConfig while leaving pdConfig as is
    void setAnalyzePhaseConfig(const uint16_t _LPType, float32_t *_pCoeffs, uint16_t _nCoeffs) {
		 setAnalyzePhaseConfig(               _LPType,            _pCoeffs,          _nCoeffs,  pdConfig);
	}
    // Set AnalyzePhaseConfig in full generality
    void setAnalyzePhaseConfig(const uint16_t _LPType, float32_t *_pCoeffs, uint16_t _nCoeffs, uint16_t _pdConfig) {
        AudioNoInterrupts(); // No interrupts while changing parameters
        LPType = _LPType;
        if (LPType == NO_LP_FILTER) {
            //Serial.println("Advice: in AnalyzePhase, for NO_LP_FILTER the output contains 2nd harmonics");
            //Serial.println("    that need external filtering.");
        }
        else if (LPType == IIR_LP_FILTER) {
            if(_pCoeffs != NULL){
                pIirCoeffs = _pCoeffs;
                nIirCoeffs = _nCoeffs;
            }
            if (nIirCoeffs != 20){
                //Serial.println("Error, in AnalyzePhase, for IIR_LP_FILTER there must be 20 coefficients.");
                nIirCoeffs = 20;
            }
            arm_biquad_cascade_df1_init_f32(&iir_inst, N_STAGES, pIirCoeffs,  &IIRStateF32[0]);
         }   
         else if (LPType==FIR_LP_FILTER) {
            if(_pCoeffs != NULL){
                pFirCoeffs = _pCoeffs;
                nFirCoeffs = _nCoeffs;
            }
			if (nFirCoeffs<4 || nFirCoeffs>NFIR_MAX) {  // Too many or too few
                //Serial.print("Error, in AnalyzePhase, for FIR_LP_FILTER there must be >4 and <=");
                //Serial.print(NFIR_MAX);
                //Serial.println(" coefficients.");
                //Serial.println("    Restoring default IIR Filter.");
                LPType = IIR_LP_FILTER;
                pIirCoeffs = &iir_coeffs[0];
                nIirCoeffs = 20;                 // Number of coefficients 20
                pdConfig = 0b11100;
                LPType = IIR_LP_FILTER;    // Variables were set in setup() above
            }
            else {     //Acceptable number, so initialize it
                arm_fir_init_f32(&fir_inst, nFirCoeffs, pFirCoeffs,  &FIRStateF32[0], block_size);
		    }
        }
        pdConfig = _pdConfig;
        AudioInterrupts(); 
    }

    void showError(uint16_t e) {
        errorPrint = e;
    }

    void update(void);

private:
    float32_t sampleRate_Hz = AUDIO_SAMPLE_RATE_EXACT;
    uint16_t block_size = AUDIO_BLOCK_SAMPLES;
    // Two input data pointers
    audio_block_f32_t *inputQueueArray_f32[2];
    // Variables controlling the configuration
    uint16_t LPType = IIR_LP_FILTER;          // NO_LP_FILTER, IIR_LP_FILTER or FIR_LP_FILTER
    float32_t *pIirCoeffs = &iir_coeffs[0];   // Coefficients for IIR
    float32_t *pFirCoeffs = &fir_coeffs[0];   // Coefficients for FIR
    uint16_t nIirCoeffs = 20;                 // Number of coefficients 20
    uint16_t nFirCoeffs = 53;                 // Number of coefficients <=200
    uint16_t  pdConfig = 0b11100;             // No limiter, fast acos, scale multiplier, radians out;
    // Control error printing in update().  Should never be enabled
    // until all audio objects have been initialized.
    // Only used as 0 or 1 now, but 16 bits are available.
    uint16_t errorPrint = 0;
    
    //  *Temporary* - TEST_TIME allows measuring time in microseconds for each part of the update()
#if TEST_TIME
    elapsedMicros tElapse;
    int32_t iitt = 998000;     // count up to a million during startup
#endif

    /* FIR filter designed with http://t-filter.appspot.com
     * Sampling frequency: 44100 Hz
     * 0 Hz - 4000 Hz   gain = 1.0, ripple = 0.101 dB
     * 7000 - 22000 Hz  attenuation >= 81.8 dB
     * Suitable for measuring phase in communications systems with linear phase.
     */
    float32_t fir_coeffs[53] = {
    -0.000206064,-0.000525129,-0.00083518, -0.000774011,  2.5925E-05,
     0.001614912, 0.003431897, 0.004335125, 0.003127158, -0.000566047,
    -0.005566484,-0.009192163,-0.008417443,-0.001801824,  0.008839149,
     0.018273049, 0.019879265, 0.009349346,-0.011696836, -0.034389317,
    -0.045008839,-0.030706279, 0.013824834, 0.082060266,  0.156328996,
     0.213799940, 0.235420817, 0.213799940, 0.156328996,  0.082060266,
     0.013824834,-0.030706279,-0.045008839,-0.034389317, -0.011696836,
     0.009349346, 0.019879265, 0.018273049, 0.008839149, -0.001801824,
    -0.008417443,-0.009192163,-0.005566484,-0.000566047,  0.003127158,
     0.004335125, 0.003431897, 0.001614912, 2.5925E-05,  -0.000774011,
    -0.000835180,-0.000525129,-0.000206064 };

    // 8-pole Biquad fc=0.0025fs, -80 dB Iowa Hills
    // This is roughly the narrowest that doesn't have
    // artifacts from numerical errors more than about
    // 0.001 radians (0.06 deg), per experiments using F32.
    // b0,b1,b2,a1,a2 for each BiQuad.  Start with stage 0
    float32_t iir_coeffs[5 * N_STAGES]={
    0.08686551007982608,
    -0.1737214710369926,
    0.08686551007982608,
    1.9951804375779567,
    -0.9951899867006161,
    // and stage 1
    0.20909791845765324,
    -0.4181667739705088,
    0.20909791845765324,
    1.9965910753714984,
    -0.9966201383162961,
    // stage 2
    0.18360046797931723,
    -0.3671514768697197,
    0.18360046797931723,
    1.9981966389027592,
    -0.998246097991674,
    // stage 3
    0.03079484444321144,
    -0.061529427044071175,
    0.03079484444321144,
    1.999421284937329,
    -0.9994815467796806};

    // ARM DSP Math library IIR filter instance
    arm_biquad_casd_df1_inst_f32 iir_inst;

    // And a FIR type, as either can be used via begin()
    arm_fir_instance_f32 fir_inst;

    // Delay line space for the FIR
    float32_t FIRStateF32[AUDIO_BLOCK_SAMPLES + NFIR_MAX];

    // Delay line space for the Biquad, each arranged as {x[n-1], x[n-2], y[n-1], y[n-2]}
    float32_t IIRStateF32[4 * N_STAGES];
};
#endif