diff --git a/AudioAnalyzePhase_F32.cpp b/AudioAnalyzePhase_F32.cpp new file mode 100644 index 0000000..da1b69d --- /dev/null +++ b/AudioAnalyzePhase_F32.cpp @@ -0,0 +1,160 @@ +/* + * AudioAnalyzePhase_F32.cpp + * + * 31 March 2020 Rev 8 April 2020 + * Bob Larkin, in support of the library: + * Chip Audette, OpenAudio, Apr 2017 + * + * * 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. + */ +/* There are two inputs, I and Q (Left and Right or 0 and 1) + * There is one output, the phase angle between the two inputs expressed in + * radians (180 degrees is Pi radians). See AudioAnalyzePhase_F32.h + * for details. + */ + +#include "AudioAnalyzePhase_F32.h" +#include +#include "mathDSP_F32.h" + +void AudioAnalyzePhase_F32::update(void) { + audio_block_f32_t *block0, *block1; + uint16_t i; + float32_t mult0, mult1, min0, max0, minmax0, min1, max1, minmax1; + float32_t min_max = 1.0f; // Filled in correctly, later + mathDSP_F32 mathDSP1; + +#if TEST_TIME + if (iitt++ >1000000) iitt = -10; + uint32_t t1, t2; + t1 = tElapse; +#endif + + // Get first input, 0 + block0 = AudioStream_F32::receiveWritable_f32(0); + if (!block0) { + if(errorPrint) Serial.println("AN_PHASE ERR: No input memory"); + return; + } + + // Get second input, 1 + block1 = AudioStream_F32::receiveWritable_f32(1); + if (!block1){ + if(errorPrint) Serial.println("AN_PHASE ERR: No input memory"); + AudioStream_F32::release(block0); + return; + } + + // Limiter on both inputs; can be used with narrow IIR filter + if (pdConfig & LIMITER_MASK) { + for (uint16_t j = 0; jlength; j++) { + if (block0->data[j]>=0.0f) + block0->data[j] = 1.0f; + else + block0->data[j] = -1.0f; + + if (block1->data[j]>=0.0f) + block1->data[j] = 1.0f; + else + block1->data[j] = -1.0f; + } + } + + // If needed, find a scaling factor for the multiplier type phase det + // This would not normally be used with a limiter---they are + // both trying to solve the same problem. + if(pdConfig & SCALE_MASK) { + min0 = 1.0E6f; max0 = -1.0E6f; + min1 = 1.0E6f; max1 = -1.0E6f; + for (i=0; i < block0->length; i++){ + if(block0->data[i] < min0) min0 = block0->data[i]; + if(block0->data[i] > max0) max0 = block0->data[i]; + if(block1->data[i] < min1) min1 = block1->data[i]; + if(block1->data[i] > max1) max1 = block1->data[i]; + } + minmax0 = max0 - min0; // 0 to 2 + minmax1 = max1 - min1; // 0 to 2 + min_max = minmax0 * minmax1; // 0 to 4 + } + + // multiply the two inputs to get phase plus second harmonic. Low pass + // filter and then apply ArcCos of the result to return the phase difference + // in radians. Multiply and leave result in 1 + arm_mult_f32(block0->data, block1->data, block1->data, block0->length); + + if (LPType == NO_LP_FILTER) { + for (i=0; i < block0->length; i++) + block0->data[i] = block1->data[i]; // Move back to block0 + } + else if(LPType == IIR_LP_FILTER) { + // Now filter 1, leaving result in 0. This is 4 BiQuads + arm_biquad_cascade_df1_f32(&iir_inst, block1->data, block0->data, block0->length); + } + else { // Alternate FIR filter for FIR_LP_FILTER + arm_fir_f32(&fir_inst, block1->data, block0->data, block0->length); + } + AudioStream_F32::release(block1); // Not needed further + + /* For variable, pdConfig: + * LIMITER_MASK: 0=No Limiter 1=Use limiter + * ACOS_MASK: 00=Use no acos linearizer 01=undefined + * 10=Fast, math-continuous acos() (default) 11=Accurate acos() + * SCALE_MASK: 0=No scale of multiplier 1=scale to min-max (default) + * UNITS_MASK: 0=Output in degrees 1=Output in radians (default) + */ + if(pdConfig & SCALE_MASK) + mult0 = 8.0f / min_max; + else + mult0 = 2.0f; + + if(pdConfig & UNITS_MASK) + mult1 = 1.0; + else + mult1 = 57.295779f; + + // Optionally, apply ArcCos() to linearize block0 output in radians * units to scale to degrees + if (pdConfig & ACOS_MASK) { + if((pdConfig & ACOS_MASK) == 0b00110) { // Accurate acosf from C-library + for (uint i = 0; idata[i] = mult1 * acosf(mult0 * block0->data[i]); + } + } + else if ((pdConfig & ACOS_MASK) == 0b00100) { // Fast acos from polynomial + for (uint i = 0; idata[i] = mult1 * mathDSP1.acos_f32(mult0 * block0->data[i]); + } + } + } + // Not applying linearization or scaling, at all. Works for multiplier outputs near 0 (90 deg difference) + else { + float32_t mult = 8.0f * mult0 * mult1; + for (uint i = 0; idata[i] = mult * block0->data[i]; + } + } + AudioStream_F32::transmit(block0, 0); + AudioStream_F32::release(block0); + #if TEST_TIME + t2 = tElapse; + if(iitt++ < 0) {Serial.print("At AnalyzePhase end, microseconds = "); Serial.println (t2 - t1); } + t1 = tElapse; + #endif +} // End update() diff --git a/AudioAnalyzePhase_F32.h b/AudioAnalyzePhase_F32.h new file mode 100644 index 0000000..fbff1fc --- /dev/null +++ b/AudioAnalyzePhase_F32.h @@ -0,0 +1,271 @@ +/* + * 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 + +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 diff --git a/AudioFilterEqualizer_F32.cpp b/AudioFilterEqualizer_F32.cpp new file mode 100644 index 0000000..3ef85c8 --- /dev/null +++ b/AudioFilterEqualizer_F32.cpp @@ -0,0 +1,192 @@ +/* AudioFilterEqualizer_F32.cpp + * + * Bob Larkin, W7PUA 8 May 2020 + * + * 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. + */ + +#include "AudioFilterEqualizer_F32.h" + +void AudioFilterEqualizer_F32::update(void) { + audio_block_f32_t *block, *block_new; + +#if TEST_TIME_EQ + if (iitt++ >1000000) iitt = -10; + uint32_t t1, t2; + t1 = tElapse; +#endif + + block = AudioStream_F32::receiveReadOnly_f32(); + if (!block) return; + + // If there's no coefficient table, give up. + if (cf32f == NULL) { + AudioStream_F32::release(block); + return; + } + + block_new = AudioStream_F32::allocate_f32(); // get a block for the FIR output + if (block_new) { + //apply the FIR + arm_fir_f32(&fir_inst, block->data, block_new->data, block->length); + AudioStream_F32::transmit(block_new); // send the FIR output + AudioStream_F32::release(block_new); + } + AudioStream_F32::release(block); + +#if TEST_TIME_EQ + t2 = tElapse; + if(iitt++ < 0) {Serial.print("At AnalyzePhase end, microseconds = "); Serial.println (t2 - t1); } + t1 = tElapse; +#endif +} + +/* equalizerNew() calculates the Equalizer FIR filter coefficients. Works from: + * uint16_t equalizerNew(uint16_t _nBands, float32_t *feq, float32_t *adb, + uint16_t _nFIR, float32_t *_cf32f, float32_t kdb) + * nBands Number of equalizer bands + * feq Pointer to array feq[] of nBands breakpoint frequencies, fractions of sample rate, Hz + * adb Pointer to array aeq[] of nBands levels, in dB, for the feq[] defined frequency bands + * nFIR The number of FIR coefficients (taps) used in the equalzer + * cf32f Pointer to an array of float to hold FIR coefficients + * kdb A parameter that trades off sidelobe levels for sharpness of band transition. + * kdb=30 sharp cutoff, poor sidelobes + * kdb=60 slow cutoff, low sidelobes + * + * The arrays, feq[], aeq[] and cf32f[] are supplied by the calling .INO + * + * Returns: 0 if successful, or an error code if not. + * Errors: 1 = Too many bands, 50 max + * 2 = sidelobe level out of range, must be > 0 + * 3 = nFIR out of range + * + * Note - This function runs at setup time, and there is no need to fret about + * processor speed. Likewise, local arrays are created on the stack and are + * available for other use when this function closes. + */ +uint16_t AudioFilterEqualizer_F32::equalizerNew(uint16_t _nBands, float32_t *feq, float32_t *adb, + uint16_t _nFIR, float32_t *_cf32f, float32_t kdb) { + uint16_t i, j; + uint16_t nHalfFIR; + float32_t beta, kbes; + float32_t q, xj2, scaleXj2, WindowWt; + float32_t fNorm[50]; // Normalized to the sampling frequency + float32_t aVolts[50]; // Convert from dB to "quasi-Volts" + mathDSP_F32 mathEqualizer; // For Bessel function + + // Make private copies + cf32f = _cf32f; + nFIR = _nFIR; + nBands = _nBands; + + // Check range of nFIR + if (nFIR<5 || nFIR>EQUALIZER_MAX_COEFFS) + return ERR_EQ_NFIR; + + // The number of FIR coefficients needs to be odd + if (2*(nFIR/2) == nFIR) + nFIR -= 1; // We just won't use the last element of the array + nHalfFIR = (nFIR - 1)/2; // If nFIR=199, nHalfFIR=99 + + for (int kk = 0; kk50) return ERR_EQ_BANDS; + for (i=0; i50) + beta = 0.1102*(kdb-8.7); + else if (kdb>20.96 && kdb<=50.0) + beta = 0.58417*powf((kdb-20.96), 0.4) + 0.07886*(kdb-20.96); + else + beta=0.0; + // Note: i0f is the floating point in & out zero'th order Bessel function (see mathDSP_F32.h) + kbes = 1.0f / mathEqualizer.i0f(beta); // An additional derived parameter used in loop + + // Apply the Kaiser window + scaleXj2 = 2.0f/(float32_t)nFIR; + scaleXj2 *= scaleXj2; + for (j=0; j<=nHalfFIR; j++) { // For 199 Taps, this is 0 to 99 + xj2 = (int16_t)(0.5f+(float32_t)j); + xj2 = scaleXj2*xj2*xj2; + WindowWt=kbes*(mathEqualizer.i0f(beta*sqrt(1.0-xj2))); + cf32f[nHalfFIR + j] *= WindowWt; // Apply the Kaiser window to upper half + cf32f[nHalfFIR - j] = cf32f[nHalfFIR +j]; // and create the lower half + } + // And fill in the members of fir_inst + arm_fir_init_f32(&fir_inst, nFIR, (float32_t *)cf32f, &StateF32[0], (uint32_t)block_size); + return 0; +} + +/* Calculate response in dB. Leave nFreq point result in array rdb[] supplied + * by the calling .INO See Parks and Burris, "Digital Filter Design," p27 (Type 1). + */ +void AudioFilterEqualizer_F32::getResponse(uint16_t nFreq, float32_t *rdb) { + uint16_t i, j; + float32_t bt; + float32_t piOnNfreq; + uint16_t nHalfFIR; + + nHalfFIR = (nFIR - 1)/2; + piOnNfreq = MF_PI / (float32_t)nFreq; + for (i=0; i1000000) iitt = -10; + uint32_t t1, t2; + t1 = tElapse; +#endif + + blockIn = AudioStream_F32::receiveReadOnly_f32(); + if (!blockIn) return; + + // If there's no coefficient table, give up. + if (cf32f == NULL) { + AudioStream_F32::release(blockIn); + return; + } + + blockOut = AudioStream_F32::allocate_f32(); // get a block for the FIR output + if (blockOut) { + // The FIR update + arm_fir_f32(&fir_inst, blockIn->data, blockOut->data, blockIn->length); + AudioStream_F32::transmit(blockOut); // send the FIR output + AudioStream_F32::release(blockOut); + } + AudioStream_F32::release(blockIn); + +#if TEST_TIME_FIRG + t2 = tElapse; + if(iitt++ < 0) {Serial.print("At FIRGeneral end, microseconds = "); Serial.println (t2 - t1); // } + Serial.print("numtaps = "); Serial.println (fir_inst.numTaps); + } + t1 = tElapse; +#endif +} + +/* FIRGeneralNew() calculates the generalFIR filter coefficients. Works from: + * adb Pointer to nFIR/2 array adb[] of levels, in dB + * nFIR The number of FIR coefficients (taps) used + * cf32f Pointer to an array of float to hold nFIR FIR coefficients + * kdb A parameter that trades off sidelobe levels for sharpness of band transition. + * kdb=30 sharp cutoff, poor sidelobes + * kdb=60 slow cutoff, low sidelobes + * pStateArray Pointer to 128+nFIR array of floats, used here briefly and then + * passed to the FIR routine as working storage + * + * The arrays, adb[], cf32f[] and pStateArray[] are supplied by the calling .INO + * + * Returns: 0 if successful, or an error code if not. + * Errors: 1 = NU + * 2 = sidelobe level out of range, must be > 0 + * 3 = nFIR out of range + * + * Note - This function runs at setup time, so slowness is not an issue + */ +uint16_t AudioFilterFIRGeneral_F32::FIRGeneralNew( + float32_t *adb, uint16_t _nFIR, float32_t *_cf32f, float32_t kdb, + float32_t *pStateArray) { + + uint16_t i, k, n; + uint16_t nHalfFIR; + float32_t beta, kbes; + float32_t num, xn2, scaleXn2, WindowWt; + float32_t M; // Burris and Parks, (2.16) + mathDSP_F32 mathEqualizer; // For Bessel function + bool even; + + cf32f = _cf32f; // Make pivate copies + nFIR = _nFIR; + + // Check range of nFIR + if (nFIR<4) + return ERR_FIRGEN_NFIR; + + M = 0.5f*(float32_t)(nFIR - 1); + // The number of FIR coefficients is even or odd + if (2*(nFIR/2) == nFIR) { + even = true; + nHalfFIR = nFIR/2; + } + else { + even = false; + nHalfFIR = (nFIR - 1)/2; + } + + for (i=0; i 0.1f) // kdb==0.0 means no window + cf32f[nHalfFIR - n - 1] *= WindowWt; // Apply window, reverse subscripts + cf32f[nHalfFIR + n] = cf32f[nHalfFIR - n - 1]; // and create the upper half + } + } + else { // nFIR is odd, nHalfFIR = (nFIR - 1)/2 + for (n=0; n<=nHalfFIR; n++) { // For 21 Taps, this is 0 to 10, including center + xn2 = (int16_t)(0.5f+(float32_t)n); + xn2 = scaleXn2*xn2*xn2; + WindowWt=kbes*(mathEqualizer.i0f(beta*sqrt(1.0-xn2))); + if(kdb > 0.1f) + cf32f[nHalfFIR - n] *= WindowWt; + // 21 taps, for n=0, nHalfFIR-n = 10, for n=1 nHalfFIR-n=9, for n=nHalfFIR, nHalfFIR-n=0 + cf32f[nHalfFIR + n] = cf32f[nHalfFIR - n]; // and create upper half (rewrite center is OK) + } + } + // And finally, fill in the members of fir_inst given in update() to the ARM FIR routine. + AudioNoInterrupts(); + arm_fir_init_f32(&fir_inst, nFIR, (float32_t *)cf32f, &pStateArray[0], (uint32_t)block_size); + AudioInterrupts(); + return 0; +} + +// FIRGeneralLoad() allows an array of nFIR FIR coefficients to be loaded. They come from an .INO +// supplied array. Also, pStateArray[] is .INO supplied and must be (block_size + nFIR) in size. +uint16_t AudioFilterFIRGeneral_F32::LoadCoeffs(uint16_t _nFIR, float32_t *_cf32f, float32_t *pStateArray) { + nFIR = _nFIR; + cf32f = _cf32f; + if (nFIR<4) // Check range of nFIR + return ERR_FIRGEN_NFIR; + for(int i=0; i<(nFIR+AUDIO_BLOCK_SAMPLES); i++) // Zero, to be sure + pStateArray[i] = 0.0f; + AudioNoInterrupts(); + arm_fir_init_f32(&fir_inst, nFIR, &cf32f[0], &pStateArray[0], (uint32_t)block_size); + AudioInterrupts(); + return 0; +} + +/* Calculate frequency response in dB. Leave nFreq point result in array rdb[] supplied + * by the calling .INO See B&P p27 (Type 1 and 2). Be aware that if nFIR*nFreq is big, + * like 100,000 or more, this will take a while to calcuulate all the cosf(). Normally, + * this is not an issue as this is an infrequent calculation. + * This function assumes that the phase of the FIR is linear with frequency, i.e., + * the coefficients are symmetrical about the middle. Otherwise, doubling of values, + * as is done here, is not valid. + */ +void AudioFilterFIRGeneral_F32::getResponse(uint16_t nFreq, float32_t *rdb) { + uint16_t i, n; + float32_t bt; + float32_t piOnNfreq; + uint16_t nHalfFIR; + float32_t M; + + piOnNfreq = MF_PI / (float32_t)nFreq; + // The number of FIR coefficients, even or odd? + if (2*(nFIR/2) == nFIR) { // it is even + nHalfFIR = nFIR/2; + M = 0.5f*(float32_t)(nFIR - 1); + for (i=0; i1000000) iitt = -10; + uint32_t t1, t2; + t1 = tElapse; +#endif + + // Get input to FM Detector block + blockIn = AudioStream_F32::receiveWritable_f32(0); + if (!blockIn) { + if(errorPrintFM) Serial.println("FMDET-ERR: No input memory"); + return; + } + + // If there's no coefficient table, give up. + if (fir_IQ_Coeffs == NULL) { + if(errorPrintFM) Serial.println("FMDET-ERR: No IQ FIR Coefficients"); + AudioStream_F32::release(blockIn); + return; + } + + if (fir_Out_Coeffs == NULL) { + if(errorPrintFM) Serial.println("FMDET-ERR: No Out FIR Coefficients"); + AudioStream_F32::release(blockIn); + return; + } + + // Try to get a block for the FM output + blockOut = AudioStream_F32::allocate_f32(); + if (!blockOut){ // Didn't have any + if(errorPrintFM) Serial.println("FMDET-ERR: No Output Memory"); + AudioStream_F32::release(blockIn); + return; + } + + // Generate sine and cosine of center frequency and double-balance mix + // these with the input signal to produce an intermediate result + // saved as v_i[] and v_q[] + for (i=0; i < block_size; i++) { + phaseS += phaseIncrement; + if (phaseS > 512.0f) + phaseS -= 512.0f; + index_sine = (uint16_t) phaseS; + deltaPhase = phaseS -(float32_t) index_sine; + /* Read two nearest values of input value from the sin table */ + a = sinTable512_f32[index_sine]; + b = sinTable512_f32[index_sine+1]; + // Linear interpolation and multiplying (DBMixer) with input + v_i[i] = blockIn->data[i] * (a + 0.001953125*(b-a)*deltaPhase); + + /* Repeat for cosine by adding 90 degrees phase */ + index_sine = (index_sine + 128) & 0x01ff; + /* Read two nearest values of input value from the sin table */ + a = sinTable512_f32[index_sine]; + b = sinTable512_f32[index_sine+1]; + /* deltaPhase will be the same as used for sin */ + v_q[i] = blockIn->data[i] * (a + 0.001953125*(b-a)*deltaPhase); + } + + // Do I FIR and Q FIR. We can borrow blockIn and blockOut at this point + //void arm_fir_f32( const arm_fir_instance_f32* S, float32_t* pSrc, float32_t* pDst, uint32_t blockSize) + arm_fir_f32(&FMDet_I_inst, v_i, blockIn->data, (uint32_t)blockIn->length); + arm_fir_f32(&FMDet_Q_inst, v_q, blockOut->data, (uint32_t)blockOut->length); + // Do ATAN2, differentiation and de-emphasis in single loop + for(i=0; idata[i], (float)blockIn->data[i]); + // Apply differentiator by subtracting last value of atan2 + if(dtemp1>MF_PI_2 && diffLast<-MF_PI_2) // Probably a wrap around + dtemp2 = dtemp1 - diffLast - MF_TWOPI; + else if(dtemp1<-MF_PI_2 && diffLast >MF_PI_2) // Probably a reverse wrap around + dtemp2 = dtemp1 - diffLast + MF_TWOPI; + else + dtemp2 = dtemp1 - diffLast; // Differentiate + diffLast = dtemp1; // Ready for next time through loop + // Data point is now dtemp2. Apply single pole de-emphasis LPF, in place + dLast = Kdem * dtemp2 + OneMinusKdem * dLast; + blockIn->data[i] = dLast; // and save to an array + } + + // Do output FIR filter. Data now in blockIn. + arm_fir_f32(&FMDet_Out_inst, blockIn->data, blockOut->data, (uint32_t)blockIn->length); + AudioStream_F32::release(blockIn); + + // Transmit the data + AudioStream_F32::transmit(blockOut, 0); + AudioStream_F32::release(blockOut); +#if TEST_TIME_FM + t2 = tElapse; + if(iitt++ < 0) {Serial.print("At end of FM Det "); Serial.println (t2 - t1); } + t1 = tElapse; +#endif +} diff --git a/RadioFMDetector_F32.h b/RadioFMDetector_F32.h new file mode 100644 index 0000000..a860f65 --- /dev/null +++ b/RadioFMDetector_F32.h @@ -0,0 +1,277 @@ +/* + * RadioFMDetector_F32 + * 22 March 2020 Bob Larkin + * With much credit to: + * Chip Audette (OpenAudio) Feb 2017 + * Building from AudioFilterFIR from Teensy Audio Library + * (AudioFilterFIR credited to Pete (El Supremo)) + * and of course, to PJRC for the Teensy and Teensy Audio Library + * + * 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. + */ + +/* This consists of a single input at some frequency, such as 10 to 20 kHz and + * an output, such as 0 to 5 kHz. The output level is linearly dependent on the + * frequency of the input sine wave frequency, i.e., an it is an FM detector. + * The input needs to be band limited below the lower frequency side of the + * input, typically 10 kHz. This is not part of this block. + * + * NOTE: Due to the sample frequencies we are working with, like 44.1 kHz, this + * detector cannot handle full FM broadcast bandwidths. It is suitable for + * NBFM as used in communications, marine radio, ham radio, etc. + * + * The output can be FIR filtered using default parameters, + * or using coefficients from an array. A separate single pole de-emphasis filer + * is included that again can be programmed. + * + * Internally, the detector uses a pair of mixers (multipliers) that generate the + * in-phase and quadrature inputs to a atan2 type of phase detector. These + * mixers have two output signals at the difference (desired) and sum (undesired) + * frequencies. The high frequency sum signal can be filtered (For a 15 kHz center, + * with an input band of 10 to 20 kHz the sum signal will be from 25 to 35 kHz that + * wraps around the 22 kHz half-sample point to produce 19 to 9 kHz. This needs to + * be removed before the atan2. A pair of FIR filters, using FIR_IQ_Coeffs + * are used. These are again programmable and default to a 29-tap LPF with + * a 5 kHz cutoff. + * + * Status: Tested static, tested with FM modulated Fluke 6061B. + * An input of about 60 microvolts to the SGTL5000 gave 12 dB SINAD. + * The output sounded good. Tested T3.6 and T4.0. No known bugs + * + * Output: Float, sensitivity is 2*pi*(f - fCenter)*sample_rate_Hz + * For 44117Hz samplerate, this is 0.000142421 per Hz + * + * Accuracy: The function used is precise. However, the approximations, such + * fastAtan2, slightly limit the accuracy. A 200 point sample of a + * 14 kHz input had an average error of 0.03 Hz + * and a standard deviation of 0.81 Hz. + * The largest errors in this sample were about +/- 1.7 Hz. This is + * with the default filters. + * + * Functions: + * frequency(float fCenter ) sets the center frequency in Hz, default 15000. + * + * filterOut(float *firCoeffs, uint nFIR, float Kdem) sets output filtering where: + * float32_t* firCoeffs is an array of coefficients + * uint nFIR is the number of coefficients + * float32_t Kdem is the de-emphasis frequency factor, where + * Kdem = 1/(0.5+(tau*fsample)) and tau is the de-emphasis + * time constant, typically 0.0005 second and fsample is + * the sample frequency, typically 44117. + * + * filterIQ(float *fir_IQ_Coeffs, uint nFIR_IQ) sets output filtering where: + * float32_t* fir_IQ_Coeffs is an array of coefficients + * uint nFIR_IQ is the number of coefficients, max 60 + * + * setSampleRate_Hz(float32_t _sampleRate_Hz) allows dynamic changing of + * the sample rate (experimental as of May 2020). + * + * returnInitializeFMError() Returns the initialization errors. + * B0001 (value 1) is an error in the IQ FIR Coefficients or quantity. + * B0010 (value 2) is an error in the Output FIR Coefficients or quantity. + * B0100 (value 4) is an error in the de-emphasis constant + * B1000 (value 8) is center frequency above half-sample frequency. + * All for debug. + * + * showError(uint16_t e) Turns error printing in the update function on (e=1) + * or off (e=0). For debug only. + * + * Time: For T3.6, an update of a 128 sample block, 370 microseconds, or + * 2.9 microseconds per data point. + * For T4.0, 87 microseconds, or 0.68 microseconds per data point. + * + * Error checking: See functions setSampleRate_Hz() and returnInitializeFMError() + * above. + */ + +#ifndef _radioFMDetector_f32_h +#define _radioFMDetector_f32_h + +#include "mathDSP_F32.h" +#include "AudioStream_F32.h" +#include "arm_math.h" + +#define MAX_FIR_IQ_COEFFS 100 +#define MAX_FIR_OUT_COEFFS 120 + +#define TEST_TIME_FM 0 + +class RadioFMDetector_F32 : public AudioStream_F32 { +//GUI: inputs:1, outputs:1 //this line used for automatic generation of GUI node +//GUI: shortName: FMDetector +public: + // Default block size and sample rate: + RadioFMDetector_F32(void) : AudioStream_F32(1, inputQueueArray_f32) { + initializeFM(); + } + // Option of AudioSettings_F32 change to block size and/or sample rate: + RadioFMDetector_F32(const AudioSettings_F32 &settings) : AudioStream_F32(1, inputQueueArray_f32) { + sampleRate_Hz = settings.sample_rate_Hz; + block_size = settings.audio_block_samples; + initializeFM(); + } + + // Provide for changing input center frequency, in Hz + void frequency(float32_t _fCenter) { + fCenter = _fCenter; + phaseIncrement = 512.0f * fCenter / sampleRate_Hz; + } + + // Provide for user FIR for I and Q signals to user supplied array + void filterIQ(float32_t* _fir_IQ_Coeffs, int _nFIR_IQ) { + if( fir_IQ_Coeffs==NULL || nFIR_IQ<4 || nFIR_IQ>MAX_FIR_IQ_COEFFS ) { + initializeFMErrors |= 1; + return; + } + fir_IQ_Coeffs = _fir_IQ_Coeffs; + nFIR_IQ = _nFIR_IQ; + initializeFM(); + } + + // Provide for changing to user FIR for detector output, (and user de-emphasis) + void filterOut(float32_t *_fir_Out_Coeffs, int _nFIR_Out, float32_t _Kdem) { + if( _fir_Out_Coeffs==NULL || _nFIR_Out<4 || _nFIR_Out>MAX_FIR_OUT_COEFFS) { + initializeFMErrors |= 2; + return; + } + if( _Kdem<0.0001 || _Kdem>1.0 ) { + initializeFMErrors |= 4; + return; + } + fir_Out_Coeffs = _fir_Out_Coeffs; + nFIR_Out = _nFIR_Out; + Kdem = _Kdem; + OneMinusKdem = 1.0f - Kdem; + initializeFM(); + } + + void setSampleRate_Hz(float32_t _sampleRate_Hz) { + if (fCenter > _sampleRate_Hz/2.0f) { // Check freq range + initializeFMErrors |= 8; + return; + } + sampleRate_Hz = _sampleRate_Hz; + // update phase increment for new frequency + phaseIncrement = 512.0f * fCenter / sampleRate_Hz; + } + + void showError(uint16_t e) { + errorPrintFM = e; + } + + uint16_t returnInitializeFMError(void) { + return initializeFMErrors; + } + + void update(void); + +private: + // One input data pointer + audio_block_f32_t *inputQueueArray_f32[1]; + float32_t fCenter = 15000.0f; + float32_t phaseS = 0.0f; + float32_t phaseS_C = 128.00f; + float32_t phaseIncrement = 512.0f*15000.0f/AUDIO_SAMPLE_RATE_EXACT; + float32_t sampleRate_Hz = AUDIO_SAMPLE_RATE_EXACT; + uint16_t block_size = AUDIO_BLOCK_SAMPLES; + // De-emphasis constant + float32_t Kdem = 0.045334f; + float32_t OneMinusKdem = 0.954666f; + // Save last data point of atan2 for differentiator + float32_t diffLast = 0.0f; + // Save last data point for next update of de-emphasis filter + float32_t dLast = 0.0f; + // Control error printing in update(), normally off + uint16_t errorPrintFM = 0; + // Monitor constructor errors + uint16_t initializeFMErrors = 0; + uint16_t nFIR_IQ = 29; + float32_t* fir_IQ_Coeffs = fir_IQ29; + uint16_t nFIR_Out = 39; + float32_t* fir_Out_Coeffs = fir_Out39; +#if TEST_TIME_FM +elapsedMicros tElapse; +int32_t iitt = 999000; // count up to a million during startup +#endif + // ARM CMSIS FIR filter instances and State vectors, sized for max, max + arm_fir_instance_f32 FMDet_I_inst; + float32_t State_I_F32[AUDIO_BLOCK_SAMPLES + MAX_FIR_IQ_COEFFS]; // 228 + + arm_fir_instance_f32 FMDet_Q_inst; + float32_t State_Q_F32[AUDIO_BLOCK_SAMPLES + MAX_FIR_IQ_COEFFS]; // 248 + + arm_fir_instance_f32 FMDet_Out_inst; + float32_t State_Out_F32[AUDIO_BLOCK_SAMPLES + MAX_FIR_OUT_COEFFS]; + + // Initialize the FM Detector, part of setting up and changing parameters + void initializeFM(void) { + if (fir_IQ_Coeffs && nFIR_IQ <= MAX_FIR_IQ_COEFFS) { +/* the instance setup call + * void arm_fir_init_f32( + * arm_fir_instance_f32* S, points to instance of floating-point FIR filter structure. + * uint16_t numTaps, Number of filter coefficients in the filter. + * float32_t* pCoeffs, points to the filter coefficients buffer. + * float32_t* pState, points to the state buffer. + * uint32_t blockSize) Number of samples that are processed per call. + */ + arm_fir_init_f32(&FMDet_I_inst, nFIR_IQ, (float32_t*)fir_IQ_Coeffs, &State_I_F32[0], (uint32_t)block_size); + arm_fir_init_f32(&FMDet_Q_inst, nFIR_IQ, (float32_t*)fir_IQ_Coeffs, &State_Q_F32[0], (uint32_t)block_size); + } + else initializeFMErrors |= B0001; + + if (fir_Out_Coeffs && nFIR_Out <= MAX_FIR_OUT_COEFFS) { + arm_fir_init_f32(&FMDet_Out_inst, nFIR_Out, (float32_t*)fir_Out_Coeffs, &State_Out_F32[0], (uint32_t)block_size); + } + else initializeFMErrors |= B0010; + dLast = 0.0; + } + + /* FIR filter designed with http://t-filter.appspot.com + * fs = 44100 Hz, < 5kHz ripple 0.29 dB, >9 kHz, -62 dB, 29 taps + */ + float32_t fir_IQ29[29] = { + -0.000970689f, -0.004690292f, -0.008256345f, -0.007565650f, + 0.001524420f, 0.015435011f, 0.021920240f, 0.008211937f, + -0.024286413f, -0.052184700f, -0.040532507f, 0.031248107f, + 0.146902412f, 0.255179564f, 0.299445269f, 0.255179564f, + 0.146902412f, 0.031248107f, -0.040532507f, -0.052184700f, + -0.024286413f, 0.008211937f, 0.021920240f, 0.015435011f, + 0.001524420f, -0.007565650f, -0.008256345f, -0.004690292f, + -0.000970689f}; + + /* FIR filter designed with http://t-filter.appspot.com + * fs = 44100 Hz, < 3kHz ripple 0.36 dB, >6 kHz, -60 dB, 39 taps + * Corrected to give DC gain = 1.00 + */ + float32_t fir_Out39[39] = { + -0.0008908477f, -0.0008401274f, -0.0001837353f, 0.0017556005f, + 0.0049353322f, 0.0084952916f, 0.0107668722f, 0.0097441685f, + 0.0039877576f, -0.0063455016f, -0.0188069300f, -0.0287453055f, + -0.0303831521f, -0.0186809770f, 0.0085931270f, 0.0493875744f, + 0.0971742012f, 0.1423015880f, 0.1745838382f, 0.1863024485f, + 0.1745838382f, 0.1423015880f, 0.0971742012f, 0.0493875744f, + 0.0085931270f, -0.0186809770f, -0.0303831521f, -0.0287453055f, + -0.0188069300f, -0.0063455016f, 0.0039877576f, 0.0097441685f, + 0.0107668722f, 0.0084952916f, 0.0049353322f, 0.0017556005f, + -0.0001837353f, -0.0008401274f, -0.0008908477f }; + +}; +#endif diff --git a/examples/AudioTestAnalyzePhase_F32/AudioTestAnalyzePhase_F32.ino b/examples/AudioTestAnalyzePhase_F32/AudioTestAnalyzePhase_F32.ino new file mode 100644 index 0000000..a46b7e0 --- /dev/null +++ b/examples/AudioTestAnalyzePhase_F32/AudioTestAnalyzePhase_F32.ino @@ -0,0 +1,114 @@ +/* Test AudioAnalyzePhase_F32.cpp RSL 7 April 2020 + * Generates 2 sine waves of different phase, but same frequency + * and measures the phase difference. + */ + +//Include files +#include +//#include +//#include +//#include +#include +#include +//#include + +//Create audio objects +// Input object creates stream, even though not used. I16 object to allow T4.x. +AudioInputI2S audioInI2S1; +AudioConvert_I16toF32 cnvt1; +// And the objects for the phase measurement: +AudioSynthWaveformSine_F32 sine1; +AudioSynthWaveformSine_F32 sine2; +AudioAnalyzePhase_F32 phase1; +AudioRecordQueue_F32 queue1_F; +AudioConnection patchCord1(audioInI2S1, 0, cnvt1, 0); +AudioConnection_F32 patchCord3(sine1, 0, phase1, 0); +AudioConnection_F32 patchCord4(sine2, 0, phase1, 1); +AudioConnection_F32 patchCord5(phase1, 0, queue1_F, 0); +AudioControlSGTL5000 sgtl5000_1; + +#define NBLOCKS 8 + float dt1[128*NBLOCKS]; + float *pq, *pd; + int k, i; + +// ==================== SETUP() ============================================ +void setup(void) { + //Start the USB serial link (to enable debugging) + Serial.begin(300); delay(500); + + AudioMemory(2); // Allocate Int16 audio data blocks + AudioMemory_F32(15); // Allocate Float32 audio data blocks + sgtl5000_1.enable(); + + AudioNoInterrupts(); + sine1.amplitude(1.0); + sine1.frequency(12345); + sine1.phase(0.0); + + sine2.amplitude(1.0); + sine2.frequency(12345); + // The next call sets the phase difference + sine2.phase(60.0); // This phase reationship will measure +60 degrees + + AudioInterrupts(); + + // The next argument can be set to 1 to show update() errors. But, it will show false + // errors before the audiostream is up and running. + phase1.showError(0); // For diagnostics + +/* For variable pdConfig (default 0b1100): + * 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 acos() + * Bit 3: 0=No scale of multiplier 1=scale to min-max (default) + * Bit 4: 0=Output in degrees 1=Output in radians (default) + * + * Uncomment one of the next 4 examples, or leave all 4 commented + * out and get the default settings for begin(), using about + * 123 microseconds for 128 block size. + * Times tu are time spent in update() on T3.6 for full 128 point block. + */ + //#1 - This uses min-max scaling, fast acos(), and FIR filtering + phase1.setAnalyzePhaseConfig(FIR_LP_FILTER, NULL, 53, 0b01100); // tu = 213 microseconds + + //#2 - This uses min-max scaling and the Accurate acos() +// phase1.setAnalyzePhaseConfig(FIR_LP_FILTER, NULL, 53, 0b01110); // tu <= 531 microseconds + + //#3 - This uses min-max scaling, IIR Filter and the no acos() linearization +// phase1.setAnalyzePhaseConfig(IIR_LP_FILTER, NULL, 53, 0b01000); // tu = 96 microseconds + + //#4 - This uses no scaling (use two magnitude 1.0 sine waves), + // and the no acos() linearization. No LP filtering +// phase1.setAnalyzePhaseConfig(NO_LP_FILTER, NULL, 20, 0b10000); // tu = 28 uSec + + i = 0; k=0; + queue1_F.begin(); + Serial.println("Setup complete."); +}; + +void loop(void) { + // Collect 128xNBLOCKS samples and output to Serial + // This "if" will be active for i on (0, NBLOCKS-1) + if (queue1_F.available() >= 1 && i>=0 && i +#include "analyze_peak_f32.h" +#include "analyze_rms_f32.h" + +AudioInputI2S i2s1; +AudioSynthWaveformSine_F32 sine1; +AudioAnalyzeRMS_F32 rms1; +AudioAnalyzePeak_F32 peak1; +AudioConnection_F32 patchCord1(sine1, 0, rms1, 0); +AudioConnection_F32 patchCord2(sine1, 0, peak1, 0); + +uint16_t n = 0; + +void setup(void) { + AudioMemory(5); //allocate Int16 audio data blocks + AudioMemory_F32(5); //allocate Float32 audio data blocks + Serial.begin(300); delay(1000); + + // Default amlitude +/- 1.0 + sine1.frequency(1000.0); + + // Set next to 0 to suppress print errors in update() + rms1.showError(1); + peak1.showError(1); +} + +void loop(void) { + if (n & 1) { + while(!rms1.available() ) ; //Wait + Serial.print("RMS value = "); + Serial.println(rms1.read(), 7 ); + while(!peak1.available() ) ; //Wait + Serial.print("Peak value = "); + Serial.println(peak1.read(), 7 ); + } + else { + while(!rms1.available() ) ; //Wait + Serial.print("RMS value = "); + Serial.println(rms1.read(), 7 ); + while(!peak1.available() ) ; //Wait + Serial.print("Peak to peak value = "); + Serial.println(peak1.readPeakToPeak(), 7 ); + } + n++; + // The RMS and Peak data collection runs during + // delay() because of hardware interrupts + delay(1000); +} diff --git a/examples/AudioTestSinCos/AudioTestSinCos.ino b/examples/AudioTestSinCos/AudioTestSinCos.ino new file mode 100644 index 0000000..f88b50b --- /dev/null +++ b/examples/AudioTestSinCos/AudioTestSinCos.ino @@ -0,0 +1,96 @@ +/* AudioTestSinCos.ino Bob Larkin 19 April 2020 + * Generates 1024 samples of Sin and Cos from the + * AudioSynthSinCos_F32. Samples sent to the USB + * Serial.print. Also, the sine and cosine signals are + * connected to the AudioAnalyzePhase_32 and a third + * Queue to Serial.print the, close to 90 degree, + * phase difference. + */ + +#include "Audio.h" +#include +#include "DSP_TeensyAudio_F32.h" + +#define NBLOCKS 8 +// Necessary foraudio stream +AudioInputI2S i2s1; +AudioConvert_I16toF32 convert; +AudioConnection pcI16(i2s1, 0, convert, 0); +// And the experiment +AudioSynthSineCosine_F32 sincos1; +AudioRecordQueue_F32 queue1; +AudioRecordQueue_F32 queue2; +AudioRecordQueue_F32 queue3; +AudioAnalyzePhase_F32 phase1; +AudioConnection_F32 patchCord1(sincos1, 0, queue1, 0); +AudioConnection_F32 patchCord2(sincos1, 1, queue2, 0); +AudioConnection_F32 patchCord3(sincos1, 0, phase1, 0); +AudioConnection_F32 patchCord4(sincos1, 1, phase1, 1); +AudioConnection_F32 patchCord5(phase1, 0, queue3, 0); + +float dt1[128*NBLOCKS]; // Place to save sin +float dt2[128*NBLOCKS]; // and cos +float dt3[128*NBLOCKS]; // and phase angle +float *pq1, *pd1, *pq2, *pd2, *pd3, *pq3; +int i, k; + +void setup(void) { + AudioMemory(5); //allocate Int16 audio data blocks + AudioMemory_F32(20); //allocate Float32 audio data blocks + Serial.begin(300); delay(1000); + + // simple() is not needed here, as it is default unless amplitude or + // phaseS_C_r is changed. But, it can be used to restore simple after + // involking changes. Here it is just a sample of usage. + sincos1.simple(true); + // Default amlitude +/- 1.0 + sincos1.frequency(1212.345); + // If either or both of the following two are uncommented, the + // detailed, slower, sincos will be used: + // sincos1.amplitude(0.9); + // sincos1.phaseS_C_r(120.00*M_PI/180.0); // Test non-90 degree, like 120 deg + + // Set next to 1 to print errors in update(), or 0 for no print. For debug only + phase1.showError(0); + + queue1.begin(); + queue2.begin(); + queue3.begin(); + i = 0; k=0; +} + +void loop(void) { + // Collect 128xNBLOCKS samples and output to Serial + // This "if" will be active for i on (0, NBLOCKS-1) + if (queue1.available() >= 1 && queue2.available() + && queue3.available() && i>=0 && i + +// Uncomment the lines "SINE" for internally generated sine wave. +// Uncomment the lines "ADC" to use the SGTL5000 Teensy audio adaptor +// AudioInputI2S_F32 i2sIn; // ADC +AudioSynthWaveformSine_F32 sine1; // SINE +RadioFMDetector_F32 fmDet1; +AudioRecordQueue_F32 queue1; +AudioOutputI2S_F32 i2sOut; // Leave in for timing +// AudioControlSGTL5000 sgtl5000_1; // ADC + +// AudioConnection_F32 connect0(i2sIn, 0, fmDet1, 0); // ADC +AudioConnection_F32 connect0(sine1, 0, fmDet1, 0); // SINE +// AudioConnection_F32 connect1(sine1, 0, queue1, 0); +AudioConnection_F32 connect3(fmDet1, 0, i2sOut, 0); +AudioConnection_F32 connect5(fmDet1, 0, queue1, 0); + +float dt1[512]; // Place to save output +float *pq1, *pd1; +uint16_t k; +int i; + +void setup(void) { + AudioMemory(5); + AudioMemory_F32(5); + Serial.begin(300); delay(1000); // Any rate is OK + Serial.println("Serial Started"); + + // sgtl5000_1.enable(); // ADC + // sgtl5000_1.inputSelect(AUDIO_INPUT_LINEIN); //ADC + + sine1.frequency(14000.0); // SINE + + // The FM detector has error checking during object construction + // when Serial.print is not available. See RadioFMDetector_F32.h: + Serial.print("FM Initialization errors: "); + Serial.println( fmDet1.returnInitializeFMError() ); + + // The following enables error checking inside of the "ubdate()" + // Output goes to the Serial (USB) Monitor. Normally, this is quiet. + fmDet1.showError(1); + + queue1.begin(); + i = 0; k=0; +} + +void loop(void) { + // Collect 512 samples and output to Serial + // This "if" will be active for i = 0,1,2,3 + if (queue1.available() >= 1) { + if( i>=0 && i<4) { + pq1 = queue1.readBuffer(); + pd1 = &dt1[i*128]; + for(k=0; k<128; k++) { + *pd1++ = *pq1++; + } + queue1.freeBuffer(); + if(i++==3) { + i=4; // Only collect 4 blocks + queue1.end(); // No more data to queue1 + } + } + else { + queue1.freeBuffer(); + } + } + // We have 512 data samples. Serial.print them + if(i == 4) { + Serial.println("512 Time in seconds and FM Output samples:"); + for (k=0; k<512; k++) { + Serial.print (0.000022667*(float32_t)k, 6); Serial.print (","); + Serial.println (dt1[k],7); + } + i = 5; + } + if(i==5) { + i = 6; + Serial.print("CPU: Percent Usage, Max: "); + Serial.print(AudioProcessorUsage()); + Serial.print(", "); + Serial.println(AudioProcessorUsageMax()); + + Serial.print("Int16 Memory Usage, Max: "); + Serial.print(AudioMemoryUsage()); + Serial.print(", "); + Serial.println(AudioMemoryUsageMax()); + + Serial.print("Float Memory Usage, Max: "); + Serial.print(AudioMemoryUsage_F32()); + Serial.print(", "); + Serial.println(AudioMemoryUsageMax_F32()); + Serial.println(); + } +} diff --git a/examples/ReceiverFM/hilbert121A.h b/examples/ReceiverFM/hilbert121A.h new file mode 100644 index 0000000..d624d04 --- /dev/null +++ b/examples/ReceiverFM/hilbert121A.h @@ -0,0 +1,123 @@ +// Following is 121 term Hilbert FIR filter +float32_t hilbert121A[121] = { + 0.000000000000000000, + 0.000773378567767513, + 0.000000000000000000, + 0.001046207887980644, + 0.000000000000000000, + 0.001368896533613985, + 0.000000000000000000, + 0.001746769975247667, + 0.000000000000000000, + 0.002185555845922462, + 0.000000000000000000, + 0.002691457154069645, + 0.000000000000000000, + 0.003271251311125927, + 0.000000000000000000, + 0.003932423233774751, + 0.000000000000000000, + 0.004683343721596901, + 0.000000000000000000, + 0.005533508538632429, + 0.000000000000000000, + 0.006493859804516438, + 0.000000000000000000, + 0.007577220484233372, + 0.000000000000000000, + 0.008798886675905997, + 0.000000000000000000, + 0.010177443901536392, + 0.000000000000000000, + 0.011735907609641917, + 0.000000000000000000, + 0.013503343224246872, + 0.000000000000000000, + 0.015517212970554440, + 0.000000000000000000, + 0.017826854793349920, + 0.000000000000000000, + 0.020498780519188083, + 0.000000000000000000, + 0.023625003856774591, + 0.000000000000000000, + 0.027336628208641155, + 0.000000000000000000, + 0.031827023036304102, + 0.000000000000000000, + 0.037393534868609392, + 0.000000000000000000, + 0.044517689704988733, + 0.000000000000000000, + 0.054032871748808158, + 0.000000000000000000, + 0.067515548043274365, + 0.000000000000000000, + 0.088347125250410385, + 0.000000000000000000, + 0.125324201622410869, + 0.000000000000000000, + 0.210709715079613419, + 0.000000000000000000, + 0.634897508268964295, + 0.000000000000000000, +-0.634897508268964295, + 0.000000000000000000, +-0.210709715079613419, + 0.000000000000000000, +-0.125324201622410869, + 0.000000000000000000, +-0.088347125250410385, + 0.000000000000000000, +-0.067515548043274365, + 0.000000000000000000, +-0.054032871748808158, + 0.000000000000000000, +-0.044517689704988733, + 0.000000000000000000, +-0.037393534868609392, + 0.000000000000000000, +-0.031827023036304102, + 0.000000000000000000, +-0.027336628208641155, + 0.000000000000000000, +-0.023625003856774591, + 0.000000000000000000, +-0.020498780519188083, + 0.000000000000000000, +-0.017826854793349920, + 0.000000000000000000, +-0.015517212970554440, + 0.000000000000000000, +-0.013503343224246872, + 0.000000000000000000, +-0.011735907609641917, + 0.000000000000000000, +-0.010177443901536392, + 0.000000000000000000, +-0.008798886675905997, + 0.000000000000000000, +-0.007577220484233372, + 0.000000000000000000, +-0.006493859804516438, + 0.000000000000000000, +-0.005533508538632429, + 0.000000000000000000, +-0.004683343721596901, + 0.000000000000000000, +-0.003932423233774751, + 0.000000000000000000, +-0.003271251311125927, + 0.000000000000000000, +-0.002691457154069645, + 0.000000000000000000, +-0.002185555845922462, + 0.000000000000000000, +-0.001746769975247667, + 0.000000000000000000, +-0.001368896533613985, + 0.000000000000000000, +-0.001046207887980644, + 0.000000000000000000, +-0.000773378567767513, + 0.000000000000000000}; diff --git a/examples/ReceiverFM/hilbert19A.h b/examples/ReceiverFM/hilbert19A.h new file mode 100644 index 0000000..155b132 --- /dev/null +++ b/examples/ReceiverFM/hilbert19A.h @@ -0,0 +1,21 @@ +// Following is a 19 term Hilbert FIR filter +float hilbert19A[]={ + 0.003006666677728199, + 0.000000000000000000, + 0.017439390086760032, + 0.000000000000000000, + 0.058235158985840196, + 0.000000000000000000, + 0.161905323559397157, + 0.000000000000000000, + 0.617838316334015092, + 0.000000000000000000, +-0.617838316334015092, + 0.000000000000000000, +-0.161905323559397157, + 0.000000000000000000, +-0.058235158985840196, + 0.000000000000000000, +-0.017439390086760032, + 0.000000000000000000, +-0.003006666677728199}; diff --git a/examples/ReceiverFM/hilbert251A.h b/examples/ReceiverFM/hilbert251A.h new file mode 100644 index 0000000..745db86 --- /dev/null +++ b/examples/ReceiverFM/hilbert251A.h @@ -0,0 +1,253 @@ +// Following is 251 term Hilbert FIR filter +float32_t hilbert251A[]={ +0.0000003255, +0.0000000000, +0.0000030702, +0.0000000000, +0.0000089286, +0.0000000000, +0.0000183061, +0.0000000000, +0.0000316287, +0.0000000000, +0.0000493436, +0.0000000000, +0.0000719193, +0.0000000000, +0.0000998451, +0.0000000000, +0.0001336320, +0.0000000000, +0.0001738120, +0.0000000000, +0.0002209393, +0.0000000000, +0.0002755899, +0.0000000000, +0.0003383625, +0.0000000000, +0.0004098790, +0.0000000000, +0.0004907853, +0.0000000000, +0.0005817525, +0.0000000000, +0.0006834782, +0.0000000000, +0.0007966881, +0.0000000000, +0.0009221383, +0.0000000000, +0.0010606178, +0.0000000000, +0.0012129515, +0.0000000000, +0.0013800041, +0.0000000000, +0.0015626848, +0.0000000000, +0.0017619529, +0.0000000000, +0.0019788241, +0.0000000000, +0.0022143787, +0.0000000000, +0.0024697715, +0.0000000000, +0.0027462425, +0.0000000000, +0.0030451312, +0.0000000000, +0.0033678928, +0.0000000000, +0.0037161183, +0.0000000000, +0.0040915578, +0.0000000000, +0.0044961498, +0.0000000000, +0.0049320558, +0.0000000000, +0.0054017033, +0.0000000000, +0.0059078375, +0.0000000000, +0.0064535860, +0.0000000000, +0.0070425380, +0.0000000000, +0.0076788436, +0.0000000000, +0.0083673390, +0.0000000000, +0.0091137048, +0.0000000000, +0.0099246683, +0.0000000000, +0.0108082660, +0.0000000000, +0.0117741868, +0.0000000000, +0.0128342256, +0.0000000000, +0.0140028938, +0.0000000000, +0.0152982506, +0.0000000000, +0.0167430570, +0.0000000000, +0.0183664064, +0.0000000000, +0.0202060801, +0.0000000000, +0.0223120327, +0.0000000000, +0.0247516963, +0.0000000000, +0.0276183140, +0.0000000000, +0.0310445375, +0.0000000000, +0.0352256211, +0.0000000000, +0.0404611696, +0.0000000000, +0.0472354231, +0.0000000000, +0.0563851215, +0.0000000000, +0.0694911881, +0.0000000000, +0.0899418673, +0.0000000000, +0.1265473875, +0.0000000000, +0.2116132716, +0.0000000000, +0.6358933477, +0.0000000000, +-0.6358933478, +0.0000000000, +-0.2116132717, +0.0000000000, +-0.1265473876, +0.0000000000, +-0.0899418674, +0.0000000000, +-0.0694911882, +0.0000000000, +-0.0563851216, +0.0000000000, +-0.0472354232, +0.0000000000, +-0.0404611697, +0.0000000000, +-0.0352256212, +0.0000000000, +-0.0310445376, +0.0000000000, +-0.0276183141, +0.0000000000, +-0.0247516964, +0.0000000000, +-0.0223120328, +0.0000000000, +-0.0202060802, +0.0000000000, +-0.0183664065, +0.0000000000, +-0.0167430571, +0.0000000000, +-0.0152982507, +0.0000000000, +-0.0140028939, +0.0000000000, +-0.0128342257, +0.0000000000, +-0.0117741869, +0.0000000000, +-0.0108082661, +0.0000000000, +-0.0099246684, +0.0000000000, +-0.0091137049, +0.0000000000, +-0.0083673391, +0.0000000000, +-0.0076788437, +0.0000000000, +-0.0070425381, +0.0000000000, +-0.0064535861, +0.0000000000, +-0.0059078376, +0.0000000000, +-0.0054017034, +0.0000000000, +-0.0049320559, +0.0000000000, +-0.0044961499, +0.0000000000, +-0.0040915579, +0.0000000000, +-0.0037161184, +0.0000000000, +-0.0033678929, +0.0000000000, +-0.0030451313, +0.0000000000, +-0.0027462426, +0.0000000000, +-0.0024697716, +0.0000000000, +-0.0022143788, +0.0000000000, +-0.0019788242, +0.0000000000, +-0.0017619530, +0.0000000000, +-0.0015626849, +0.0000000000, +-0.0013800042, +0.0000000000, +-0.0012129516, +0.0000000000, +-0.0010606179, +0.0000000000, +-0.0009221384, +0.0000000000, +-0.0007966882, +0.0000000000, +-0.0006834783, +0.0000000000, +-0.0005817526, +0.0000000000, +-0.0004907854, +0.0000000000, +-0.0004098791, +0.0000000000, +-0.0003383626, +0.0000000000, +-0.0002755900, +0.0000000000, +-0.0002209394, +0.0000000000, +-0.0001738121, +0.0000000000, +-0.0001336321, +0.0000000000, +-0.0000998452, +0.0000000000, +-0.0000719194, +0.0000000000, +-0.0000493437, +0.0000000000, +-0.0000316288, +0.0000000000, +-0.0000183062, +0.0000000000, +-0.0000089287, +0.0000000000, +-0.0000030703, +0.0000000000, +-0.0000003256}; diff --git a/examples/ReceiverPart1/Rcvr1Outputs.gnumeric b/examples/ReceiverPart1/Rcvr1Outputs.gnumeric new file mode 100644 index 0000000..d8b30c9 Binary files /dev/null and b/examples/ReceiverPart1/Rcvr1Outputs.gnumeric differ diff --git a/examples/ReceiverPart1/ReceiverPart1.ino b/examples/ReceiverPart1/ReceiverPart1.ino new file mode 100644 index 0000000..855a489 --- /dev/null +++ b/examples/ReceiverPart1/ReceiverPart1.ino @@ -0,0 +1,90 @@ +/* ReceiverPart1.ino Bob Larkin 20 April 2020 + * This is "Hello World" of Radio design. It can + * receive a Lower Sideband signal, but that is it. + * To test the floating point radio receiver blocks: + * AudioFilter90Deg_F32 + * RadioIQMixer_F32 + */ + +#include "Audio.h" +#include +#include "DSP_TeensyAudio_F32.h" + +AudioInputI2S i2s1; +AudioSynthWaveformSine_F32 sine1; // Test signal +RadioIQMixer_F32 iqmixer1; +AudioFilter90Deg_F32 hilbert1; +// NOTE: Use AudioMixer4_F32 *from Tympan* +AudioMixer4_F32 sum1; +AudioRecordQueue_F32 queue1; // The LSB output +AudioRecordQueue_F32 queue2; // The test sine wave + +AudioConnection_F32 patchCord0(sine1, 0, iqmixer1, 0); +AudioConnection_F32 patchCord2(sine1, 0, iqmixer1, 1); +AudioConnection_F32 patchCord4(iqmixer1, 0, hilbert1, 0); +AudioConnection_F32 patchCord5(iqmixer1, 1, hilbert1, 1); +AudioConnection_F32 patchCord6(hilbert1, 0, sum1, 0); +AudioConnection_F32 patchCord7(hilbert1, 1, sum1, 1); +AudioConnection_F32 patchCord8(sum1, 0, queue1, 0); +AudioConnection_F32 patchCord9(sine1, 0, queue2, 0); + +// Pick one of these +#include "hilbert19A.h" +#include "hilbert121A.h" +#include "hilbert251A.h" + +float dt1[128]; +float dt2[128]; +float *pq1, *pd1, *pq2, *pd2; +int i, k; + +void setup(void) { + AudioMemory(5); + AudioMemory_F32(10); + Serial.begin(300); delay(1000); + sine1.frequency(14000.0); + iqmixer1.frequency(16000.0); + // Pick one of the three, the 251 does not have + // enough samples here to show the full build-up. + // hilbert1.begin(hilbert19A, 19); + // hilbert1.begin(hilbert121A, 121); + hilbert1.begin(hilbert121A, 121); + sum1.gain(0, 1.0); // Leave at 1.0 + sum1.gain(1, -1.0); // -1 for LSB out + iqmixer1.showError(1); // Prints update() errors + hilbert1.showError(1); + queue1.begin(); + queue2.begin(); + i = 0; k=0; +} + +void loop(void) { + // Collect 128 samples and output to Serial + // This "if" will be active for i == 0 + if (queue1.available() >= 1 && i==0) { + pq1 = queue1.readBuffer(); + pd1 = &dt1[0]; + pq2 = queue2.readBuffer(); + pd2 = &dt2[0]; + for(k=0; k<128; k++) { + *pd1++ = *pq1++; + *pd2++ = *pq2++; + } + i=1; // Only collect 1 block + Serial.print("Maximum F32 memory usage at loop:"); + Serial.println( AudioMemoryUsageMax_F32() ); + queue1.freeBuffer(); + queue2.freeBuffer(); + queue1.end(); // No more data to queue1 + queue2.end(); // No more data to queue2 + } + if(i == 1) { + Serial.println("128 Samples: Time (sec), LSB Out, Sine Wave In"); + for (k=0; k<128; k++) { + Serial.print (0.000022667*(float32_t)k, 6); Serial.print (","); + Serial.print (dt1[k],7); Serial.print (","); + Serial.println (dt2[k],7); + } + i = 2; + } +} diff --git a/examples/ReceiverPart1/hilbert121A.h b/examples/ReceiverPart1/hilbert121A.h new file mode 100644 index 0000000..d624d04 --- /dev/null +++ b/examples/ReceiverPart1/hilbert121A.h @@ -0,0 +1,123 @@ +// Following is 121 term Hilbert FIR filter +float32_t hilbert121A[121] = { + 0.000000000000000000, + 0.000773378567767513, + 0.000000000000000000, + 0.001046207887980644, + 0.000000000000000000, + 0.001368896533613985, + 0.000000000000000000, + 0.001746769975247667, + 0.000000000000000000, + 0.002185555845922462, + 0.000000000000000000, + 0.002691457154069645, + 0.000000000000000000, + 0.003271251311125927, + 0.000000000000000000, + 0.003932423233774751, + 0.000000000000000000, + 0.004683343721596901, + 0.000000000000000000, + 0.005533508538632429, + 0.000000000000000000, + 0.006493859804516438, + 0.000000000000000000, + 0.007577220484233372, + 0.000000000000000000, + 0.008798886675905997, + 0.000000000000000000, + 0.010177443901536392, + 0.000000000000000000, + 0.011735907609641917, + 0.000000000000000000, + 0.013503343224246872, + 0.000000000000000000, + 0.015517212970554440, + 0.000000000000000000, + 0.017826854793349920, + 0.000000000000000000, + 0.020498780519188083, + 0.000000000000000000, + 0.023625003856774591, + 0.000000000000000000, + 0.027336628208641155, + 0.000000000000000000, + 0.031827023036304102, + 0.000000000000000000, + 0.037393534868609392, + 0.000000000000000000, + 0.044517689704988733, + 0.000000000000000000, + 0.054032871748808158, + 0.000000000000000000, + 0.067515548043274365, + 0.000000000000000000, + 0.088347125250410385, + 0.000000000000000000, + 0.125324201622410869, + 0.000000000000000000, + 0.210709715079613419, + 0.000000000000000000, + 0.634897508268964295, + 0.000000000000000000, +-0.634897508268964295, + 0.000000000000000000, +-0.210709715079613419, + 0.000000000000000000, +-0.125324201622410869, + 0.000000000000000000, +-0.088347125250410385, + 0.000000000000000000, +-0.067515548043274365, + 0.000000000000000000, +-0.054032871748808158, + 0.000000000000000000, +-0.044517689704988733, + 0.000000000000000000, +-0.037393534868609392, + 0.000000000000000000, +-0.031827023036304102, + 0.000000000000000000, +-0.027336628208641155, + 0.000000000000000000, +-0.023625003856774591, + 0.000000000000000000, +-0.020498780519188083, + 0.000000000000000000, +-0.017826854793349920, + 0.000000000000000000, +-0.015517212970554440, + 0.000000000000000000, +-0.013503343224246872, + 0.000000000000000000, +-0.011735907609641917, + 0.000000000000000000, +-0.010177443901536392, + 0.000000000000000000, +-0.008798886675905997, + 0.000000000000000000, +-0.007577220484233372, + 0.000000000000000000, +-0.006493859804516438, + 0.000000000000000000, +-0.005533508538632429, + 0.000000000000000000, +-0.004683343721596901, + 0.000000000000000000, +-0.003932423233774751, + 0.000000000000000000, +-0.003271251311125927, + 0.000000000000000000, +-0.002691457154069645, + 0.000000000000000000, +-0.002185555845922462, + 0.000000000000000000, +-0.001746769975247667, + 0.000000000000000000, +-0.001368896533613985, + 0.000000000000000000, +-0.001046207887980644, + 0.000000000000000000, +-0.000773378567767513, + 0.000000000000000000}; diff --git a/examples/ReceiverPart1/hilbert19A.h b/examples/ReceiverPart1/hilbert19A.h new file mode 100644 index 0000000..155b132 --- /dev/null +++ b/examples/ReceiverPart1/hilbert19A.h @@ -0,0 +1,21 @@ +// Following is a 19 term Hilbert FIR filter +float hilbert19A[]={ + 0.003006666677728199, + 0.000000000000000000, + 0.017439390086760032, + 0.000000000000000000, + 0.058235158985840196, + 0.000000000000000000, + 0.161905323559397157, + 0.000000000000000000, + 0.617838316334015092, + 0.000000000000000000, +-0.617838316334015092, + 0.000000000000000000, +-0.161905323559397157, + 0.000000000000000000, +-0.058235158985840196, + 0.000000000000000000, +-0.017439390086760032, + 0.000000000000000000, +-0.003006666677728199}; diff --git a/examples/ReceiverPart1/hilbert251A.h b/examples/ReceiverPart1/hilbert251A.h new file mode 100644 index 0000000..745db86 --- /dev/null +++ b/examples/ReceiverPart1/hilbert251A.h @@ -0,0 +1,253 @@ +// Following is 251 term Hilbert FIR filter +float32_t hilbert251A[]={ +0.0000003255, +0.0000000000, +0.0000030702, +0.0000000000, +0.0000089286, +0.0000000000, +0.0000183061, +0.0000000000, +0.0000316287, +0.0000000000, +0.0000493436, +0.0000000000, +0.0000719193, +0.0000000000, +0.0000998451, +0.0000000000, +0.0001336320, +0.0000000000, +0.0001738120, +0.0000000000, +0.0002209393, +0.0000000000, +0.0002755899, +0.0000000000, +0.0003383625, +0.0000000000, +0.0004098790, +0.0000000000, +0.0004907853, +0.0000000000, +0.0005817525, +0.0000000000, +0.0006834782, +0.0000000000, +0.0007966881, +0.0000000000, +0.0009221383, +0.0000000000, +0.0010606178, +0.0000000000, +0.0012129515, +0.0000000000, +0.0013800041, +0.0000000000, +0.0015626848, +0.0000000000, +0.0017619529, +0.0000000000, +0.0019788241, +0.0000000000, +0.0022143787, +0.0000000000, +0.0024697715, +0.0000000000, +0.0027462425, +0.0000000000, +0.0030451312, +0.0000000000, +0.0033678928, +0.0000000000, +0.0037161183, +0.0000000000, +0.0040915578, +0.0000000000, +0.0044961498, +0.0000000000, +0.0049320558, +0.0000000000, +0.0054017033, +0.0000000000, +0.0059078375, +0.0000000000, +0.0064535860, +0.0000000000, +0.0070425380, +0.0000000000, +0.0076788436, +0.0000000000, +0.0083673390, +0.0000000000, +0.0091137048, +0.0000000000, +0.0099246683, +0.0000000000, +0.0108082660, +0.0000000000, +0.0117741868, +0.0000000000, +0.0128342256, +0.0000000000, +0.0140028938, +0.0000000000, +0.0152982506, +0.0000000000, +0.0167430570, +0.0000000000, +0.0183664064, +0.0000000000, +0.0202060801, +0.0000000000, +0.0223120327, +0.0000000000, +0.0247516963, +0.0000000000, +0.0276183140, +0.0000000000, +0.0310445375, +0.0000000000, +0.0352256211, +0.0000000000, +0.0404611696, +0.0000000000, +0.0472354231, +0.0000000000, +0.0563851215, +0.0000000000, +0.0694911881, +0.0000000000, +0.0899418673, +0.0000000000, +0.1265473875, +0.0000000000, +0.2116132716, +0.0000000000, +0.6358933477, +0.0000000000, +-0.6358933478, +0.0000000000, +-0.2116132717, +0.0000000000, +-0.1265473876, +0.0000000000, +-0.0899418674, +0.0000000000, +-0.0694911882, +0.0000000000, +-0.0563851216, +0.0000000000, +-0.0472354232, +0.0000000000, +-0.0404611697, +0.0000000000, +-0.0352256212, +0.0000000000, +-0.0310445376, +0.0000000000, +-0.0276183141, +0.0000000000, +-0.0247516964, +0.0000000000, +-0.0223120328, +0.0000000000, +-0.0202060802, +0.0000000000, +-0.0183664065, +0.0000000000, +-0.0167430571, +0.0000000000, +-0.0152982507, +0.0000000000, +-0.0140028939, +0.0000000000, +-0.0128342257, +0.0000000000, +-0.0117741869, +0.0000000000, +-0.0108082661, +0.0000000000, +-0.0099246684, +0.0000000000, +-0.0091137049, +0.0000000000, +-0.0083673391, +0.0000000000, +-0.0076788437, +0.0000000000, +-0.0070425381, +0.0000000000, +-0.0064535861, +0.0000000000, +-0.0059078376, +0.0000000000, +-0.0054017034, +0.0000000000, +-0.0049320559, +0.0000000000, +-0.0044961499, +0.0000000000, +-0.0040915579, +0.0000000000, +-0.0037161184, +0.0000000000, +-0.0033678929, +0.0000000000, +-0.0030451313, +0.0000000000, +-0.0027462426, +0.0000000000, +-0.0024697716, +0.0000000000, +-0.0022143788, +0.0000000000, +-0.0019788242, +0.0000000000, +-0.0017619530, +0.0000000000, +-0.0015626849, +0.0000000000, +-0.0013800042, +0.0000000000, +-0.0012129516, +0.0000000000, +-0.0010606179, +0.0000000000, +-0.0009221384, +0.0000000000, +-0.0007966882, +0.0000000000, +-0.0006834783, +0.0000000000, +-0.0005817526, +0.0000000000, +-0.0004907854, +0.0000000000, +-0.0004098791, +0.0000000000, +-0.0003383626, +0.0000000000, +-0.0002755900, +0.0000000000, +-0.0002209394, +0.0000000000, +-0.0001738121, +0.0000000000, +-0.0001336321, +0.0000000000, +-0.0000998452, +0.0000000000, +-0.0000719194, +0.0000000000, +-0.0000493437, +0.0000000000, +-0.0000316288, +0.0000000000, +-0.0000183062, +0.0000000000, +-0.0000089287, +0.0000000000, +-0.0000030703, +0.0000000000, +-0.0000003256}; diff --git a/examples/ReceiverPart2/ReceiverPart2.ino b/examples/ReceiverPart2/ReceiverPart2.ino new file mode 100644 index 0000000..83125ce --- /dev/null +++ b/examples/ReceiverPart2/ReceiverPart2.ino @@ -0,0 +1,177 @@ +/* ReceiverPart2.ino Bob Larkin 29 April 2020 + * This is a simple SP radio design. It can receive 2 modes, + * Single Sideband (SSB) and Narrow Band FM (NBFM). SSB + * breaks into Lower Sidband (LSB) and Upper Sideband (USB). + * It gets even better in that AM can be received on either + * LSB or USB by tuning to the AM carrier frequency. + * + * The signal path is switched between SSB and FM by a Chip + * Audette AudioSwitch4_F32 switch block. This keeps resources + * from being used in blocks that are not needed. + * + * We are restricted to receiving in the 8 to 22 kHz range, + * so to use this on the air would require frequency conversion + * hardware. + * + * Details on the blocks are part of the include .h files for the + * blocks and in the INO files + * TestFM.ino + * ReceiverPart1.ino + * + * Input and Output is via the Teensy Audio Adaptor that uses + * a SGTL5000 CODEC, + * + * See the individual .h files for each block for more details. + * + * Measured peak-to-peak levels, using AudioAnalyzePeak_F32, + * all done at overload point for ADC: + * At ADC (i2sIn), max, 2.0 + * At IQ Mixer out, max, 2.0 + * At 90 deg Phase out, 2.0 + * At Sum out 2.0 + * At LPF FIR Out 2.0 + * + * FM Det gives 0.50 out for about 5.6 kHz p-p deviation + * + * T3.6 Processor load, measured: 16% for NBFM + * 30% for LSB or USB 29 tap LPF + * T4.0 Processor load, measured: 4.3% for NBFM + * 6.5% for LSB or USB 29 tap LPF + */ + +#include "Audio.h" +#include +#include "DSP_TeensyAudio_F32.h" + +// ********* Mini Control Panel ********* +// Set mode and gain here and re-compile + +// Here is the mode switch +#define LSB 1 +#define USB 2 +#define NBFM 3 +uint16_t mode = NBFM; // <--Select mode + +int gainControlDB = 0; // Set SSB gain in dB. 0 dB is a gain of 1.0 + +// ***************************************** +// To work with T4.0 the I2S routine outputs 16-bit integer (I16). Then +// use Audette I16 to F32 convert. Same below for output, in reverse. +AudioInputI2S i2sIn; +AudioConvert_I16toF32 cnvrt1; +AudioSwitch4_F32 switch1; // Select SSB or FM +RadioIQMixer_F32 iqmixer1; +AudioFilter90Deg_F32 hilbert1; +AudioMixer4_F32 sum1; // Summing node for the SSB receiver +AudioFilterFIR_F32 fir1; // Low Pass Filter to frequency limit the SSB +RadioFMDetector_F32 fmdet1; // NBFM from 10 to 20 kHz +AudioMixer4_F32 sum2; // SSB and NBFM rejoin here +AudioConvert_F32toI16 cnvrt2; // Left +AudioConvert_F32toI16 cnvrt3; // Right +AudioOutputI2S i2sOut; +AudioAnalyzePeak_F32 peak1; +AudioControlSGTL5000 sgtl5000_1; + +AudioConnection conI16_1(i2sIn, 0, cnvrt1, 0); // ADC +AudioConnection_F32 connect0(cnvrt1, 0, switch1, 0); // Analog to Digital +AudioConnection_F32 connect1(switch1, 0, iqmixer1, 0); // SSB Input +AudioConnection_F32 connect2(switch1, 0, iqmixer1, 1); // SSB Input +AudioConnection_F32 connect3(switch1, 1, fmdet1, 0); // FM input +AudioConnection_F32 connect4(iqmixer1, 0, hilbert1, 0); // Broadband 90 deg phase +AudioConnection_F32 connect5(iqmixer1, 1, hilbert1, 1); +AudioConnection_F32 connect6(hilbert1, 0, sum1, 0); // Sideband select +AudioConnection_F32 connect7(hilbert1, 1, sum1, 1); +AudioConnection_F32 connect8(sum1, 0, fir1, 0); // Limit audio BW +AudioConnection_F32 connect9(fir1, 0, sum2, 0); // Output of SSB +AudioConnection_F32 connectA(fmdet1, 0, sum2, 1); // Output of FM +AudioConnection_F32 connectC(sum2, 0, cnvrt2, 0); // Out to the CODEC left +AudioConnection_F32 connectD(sum2, 0, cnvrt3, 0); // and right +AudioConnection_F32 connectE(sum2, 0, peak1, 0); +AudioConnection conI16_2(cnvrt2, 0, i2sOut, 0); // DAC +AudioConnection conI16_3(cnvrt3, 0, i2sOut, 1); // DAC + +// Filter for AudioFilter90Deg_F32 hilbert1 +#include "hilbert251A.h" + +/* FIR filter designed with http://t-filter.appspot.com + * fs = 44100 Hz, < 5kHz ripple 0.29 dB, >9 kHz, -62 dB, 29 taps + */ +float32_t fir_IQ29[29] = { +-0.000970689f, -0.004690292f, -0.008256345f, -0.007565650f, + 0.001524420f, 0.015435011f, 0.021920240f, 0.008211937f, +-0.024286413f, -0.052184700f, -0.040532507f, 0.031248107f, + 0.146902412f, 0.255179564f, 0.299445269f, 0.255179564f, + 0.146902412f, 0.031248107f, -0.040532507f, -0.052184700f, +-0.024286413f, 0.008211937f, 0.021920240f, 0.015435011f, + 0.001524420f, -0.007565650f, -0.008256345f, -0.004690292f, +-0.000970689f}; + +void setup(void) { + float32_t vGain; + AudioMemory(5); + AudioMemory_F32(5); + Serial.begin(300); delay(1000); + + // Enable the audio shield, select input, and enable output + sgtl5000_1.enable(); // ADC + sgtl5000_1.inputSelect(AUDIO_INPUT_LINEIN); + + // The switch is single pole 4 position, numbered (0, 3) 0=SSB, 1 = FM + if(mode==LSB || mode==USB){ switch1.setChannel(0); Serial.println("SSB"); } + else if(mode==NBFM) { switch1.setChannel(1); Serial.println("NBFM"); } + + iqmixer1.frequency(15000.0); // LO Frequency, typically 10 to 15 kHz + + hilbert1.begin(hilbert251A, 251); // Set the Hilbert transform FIR filter + + sum1.gain(0, 1.0f); // Leave set at 1.0 + if(mode==LSB) sum1.gain(1, -1.0f); // -1 for LSB out + else if(mode==USB) sum1.gain(1, 1.0f); // +1 for USB and for NBFM, we don't care + + // The FM detector has error checking during object construction + // when Serial.print is not available. See RadioFMDetector_F32.h: + Serial.print("FM Initialization errors: "); + Serial.println( fmdet1.returnInitializeFMError() ); + + // See RadioFMDetector_F32.h for information on functions for modifying the + // FM Detector. Default values are used here, starting with a 15 kHz center frequency. + + fir1.begin(fir_IQ29, 29); // 5 kHz bandwidth set for radio in SSB + + // The gainControlDB goes in 1 dB steps. Convert here to a voltage ratio + vGain = powf(10.0f, ((float32_t)gainControlDB)/20.0 ); + // And apply that ratio to the output summing block. Gain here for SSB only + sum2.gain(0, vGain); Serial.print("vGain = "); Serial.println(vGain, 4); + sum2.gain(1, 1.0f); // FM gain + + // The following enable error checking inside of the blocks indicated. + // Output goes to the Serial (USB) Monitor. Use for debug. + //hilbert1.showError(1); // Should show input error when in FM + //fmdet1.showError(1); // Should show input error when in LSB or USB +} + +void loop(void) { +// Here is where the adjustment of the volume control could go. +// And anything else that needs regular attention, other +// than the audio stream. + + if (peak1.available() ) {Serial.print("P-P ="); Serial.println(peak1.readPeakToPeak(), 6);} + else Serial.println("Peak-Peak not available"); + Serial.print("CPU: Percent Usage, Max: "); + Serial.print(AudioProcessorUsage()); + Serial.print(", "); + Serial.print(AudioProcessorUsageMax()); + Serial.print(" "); + Serial.print("Int16 Memory: "); + Serial.print(AudioMemoryUsage()); + Serial.print(", "); + Serial.print(AudioMemoryUsageMax()); + Serial.print(" "); + Serial.print("Float Memory: "); + Serial.print(AudioMemoryUsage_F32()); + Serial.print(", "); + Serial.println(AudioMemoryUsageMax_F32()); + Serial.println(); + + delay(1000); +} diff --git a/examples/ReceiverPart2/hilbert121A.h b/examples/ReceiverPart2/hilbert121A.h new file mode 100644 index 0000000..d624d04 --- /dev/null +++ b/examples/ReceiverPart2/hilbert121A.h @@ -0,0 +1,123 @@ +// Following is 121 term Hilbert FIR filter +float32_t hilbert121A[121] = { + 0.000000000000000000, + 0.000773378567767513, + 0.000000000000000000, + 0.001046207887980644, + 0.000000000000000000, + 0.001368896533613985, + 0.000000000000000000, + 0.001746769975247667, + 0.000000000000000000, + 0.002185555845922462, + 0.000000000000000000, + 0.002691457154069645, + 0.000000000000000000, + 0.003271251311125927, + 0.000000000000000000, + 0.003932423233774751, + 0.000000000000000000, + 0.004683343721596901, + 0.000000000000000000, + 0.005533508538632429, + 0.000000000000000000, + 0.006493859804516438, + 0.000000000000000000, + 0.007577220484233372, + 0.000000000000000000, + 0.008798886675905997, + 0.000000000000000000, + 0.010177443901536392, + 0.000000000000000000, + 0.011735907609641917, + 0.000000000000000000, + 0.013503343224246872, + 0.000000000000000000, + 0.015517212970554440, + 0.000000000000000000, + 0.017826854793349920, + 0.000000000000000000, + 0.020498780519188083, + 0.000000000000000000, + 0.023625003856774591, + 0.000000000000000000, + 0.027336628208641155, + 0.000000000000000000, + 0.031827023036304102, + 0.000000000000000000, + 0.037393534868609392, + 0.000000000000000000, + 0.044517689704988733, + 0.000000000000000000, + 0.054032871748808158, + 0.000000000000000000, + 0.067515548043274365, + 0.000000000000000000, + 0.088347125250410385, + 0.000000000000000000, + 0.125324201622410869, + 0.000000000000000000, + 0.210709715079613419, + 0.000000000000000000, + 0.634897508268964295, + 0.000000000000000000, +-0.634897508268964295, + 0.000000000000000000, +-0.210709715079613419, + 0.000000000000000000, +-0.125324201622410869, + 0.000000000000000000, +-0.088347125250410385, + 0.000000000000000000, +-0.067515548043274365, + 0.000000000000000000, +-0.054032871748808158, + 0.000000000000000000, +-0.044517689704988733, + 0.000000000000000000, +-0.037393534868609392, + 0.000000000000000000, +-0.031827023036304102, + 0.000000000000000000, +-0.027336628208641155, + 0.000000000000000000, +-0.023625003856774591, + 0.000000000000000000, +-0.020498780519188083, + 0.000000000000000000, +-0.017826854793349920, + 0.000000000000000000, +-0.015517212970554440, + 0.000000000000000000, +-0.013503343224246872, + 0.000000000000000000, +-0.011735907609641917, + 0.000000000000000000, +-0.010177443901536392, + 0.000000000000000000, +-0.008798886675905997, + 0.000000000000000000, +-0.007577220484233372, + 0.000000000000000000, +-0.006493859804516438, + 0.000000000000000000, +-0.005533508538632429, + 0.000000000000000000, +-0.004683343721596901, + 0.000000000000000000, +-0.003932423233774751, + 0.000000000000000000, +-0.003271251311125927, + 0.000000000000000000, +-0.002691457154069645, + 0.000000000000000000, +-0.002185555845922462, + 0.000000000000000000, +-0.001746769975247667, + 0.000000000000000000, +-0.001368896533613985, + 0.000000000000000000, +-0.001046207887980644, + 0.000000000000000000, +-0.000773378567767513, + 0.000000000000000000}; diff --git a/examples/ReceiverPart2/hilbert19A.h b/examples/ReceiverPart2/hilbert19A.h new file mode 100644 index 0000000..155b132 --- /dev/null +++ b/examples/ReceiverPart2/hilbert19A.h @@ -0,0 +1,21 @@ +// Following is a 19 term Hilbert FIR filter +float hilbert19A[]={ + 0.003006666677728199, + 0.000000000000000000, + 0.017439390086760032, + 0.000000000000000000, + 0.058235158985840196, + 0.000000000000000000, + 0.161905323559397157, + 0.000000000000000000, + 0.617838316334015092, + 0.000000000000000000, +-0.617838316334015092, + 0.000000000000000000, +-0.161905323559397157, + 0.000000000000000000, +-0.058235158985840196, + 0.000000000000000000, +-0.017439390086760032, + 0.000000000000000000, +-0.003006666677728199}; diff --git a/examples/ReceiverPart2/hilbert251A.h b/examples/ReceiverPart2/hilbert251A.h new file mode 100644 index 0000000..745db86 --- /dev/null +++ b/examples/ReceiverPart2/hilbert251A.h @@ -0,0 +1,253 @@ +// Following is 251 term Hilbert FIR filter +float32_t hilbert251A[]={ +0.0000003255, +0.0000000000, +0.0000030702, +0.0000000000, +0.0000089286, +0.0000000000, +0.0000183061, +0.0000000000, +0.0000316287, +0.0000000000, +0.0000493436, +0.0000000000, +0.0000719193, +0.0000000000, +0.0000998451, +0.0000000000, +0.0001336320, +0.0000000000, +0.0001738120, +0.0000000000, +0.0002209393, +0.0000000000, +0.0002755899, +0.0000000000, +0.0003383625, +0.0000000000, +0.0004098790, +0.0000000000, +0.0004907853, +0.0000000000, +0.0005817525, +0.0000000000, +0.0006834782, +0.0000000000, +0.0007966881, +0.0000000000, +0.0009221383, +0.0000000000, +0.0010606178, +0.0000000000, +0.0012129515, +0.0000000000, +0.0013800041, +0.0000000000, +0.0015626848, +0.0000000000, +0.0017619529, +0.0000000000, +0.0019788241, +0.0000000000, +0.0022143787, +0.0000000000, +0.0024697715, +0.0000000000, +0.0027462425, +0.0000000000, +0.0030451312, +0.0000000000, +0.0033678928, +0.0000000000, +0.0037161183, +0.0000000000, +0.0040915578, +0.0000000000, +0.0044961498, +0.0000000000, +0.0049320558, +0.0000000000, +0.0054017033, +0.0000000000, +0.0059078375, +0.0000000000, +0.0064535860, +0.0000000000, +0.0070425380, +0.0000000000, +0.0076788436, +0.0000000000, +0.0083673390, +0.0000000000, +0.0091137048, +0.0000000000, +0.0099246683, +0.0000000000, +0.0108082660, +0.0000000000, +0.0117741868, +0.0000000000, +0.0128342256, +0.0000000000, +0.0140028938, +0.0000000000, +0.0152982506, +0.0000000000, +0.0167430570, +0.0000000000, +0.0183664064, +0.0000000000, +0.0202060801, +0.0000000000, +0.0223120327, +0.0000000000, +0.0247516963, +0.0000000000, +0.0276183140, +0.0000000000, +0.0310445375, +0.0000000000, +0.0352256211, +0.0000000000, +0.0404611696, +0.0000000000, +0.0472354231, +0.0000000000, +0.0563851215, +0.0000000000, +0.0694911881, +0.0000000000, +0.0899418673, +0.0000000000, +0.1265473875, +0.0000000000, +0.2116132716, +0.0000000000, +0.6358933477, +0.0000000000, +-0.6358933478, +0.0000000000, +-0.2116132717, +0.0000000000, +-0.1265473876, +0.0000000000, +-0.0899418674, +0.0000000000, +-0.0694911882, +0.0000000000, +-0.0563851216, +0.0000000000, +-0.0472354232, +0.0000000000, +-0.0404611697, +0.0000000000, +-0.0352256212, +0.0000000000, +-0.0310445376, +0.0000000000, +-0.0276183141, +0.0000000000, +-0.0247516964, +0.0000000000, +-0.0223120328, +0.0000000000, +-0.0202060802, +0.0000000000, +-0.0183664065, +0.0000000000, +-0.0167430571, +0.0000000000, +-0.0152982507, +0.0000000000, +-0.0140028939, +0.0000000000, +-0.0128342257, +0.0000000000, +-0.0117741869, +0.0000000000, +-0.0108082661, +0.0000000000, +-0.0099246684, +0.0000000000, +-0.0091137049, +0.0000000000, +-0.0083673391, +0.0000000000, +-0.0076788437, +0.0000000000, +-0.0070425381, +0.0000000000, +-0.0064535861, +0.0000000000, +-0.0059078376, +0.0000000000, +-0.0054017034, +0.0000000000, +-0.0049320559, +0.0000000000, +-0.0044961499, +0.0000000000, +-0.0040915579, +0.0000000000, +-0.0037161184, +0.0000000000, +-0.0033678929, +0.0000000000, +-0.0030451313, +0.0000000000, +-0.0027462426, +0.0000000000, +-0.0024697716, +0.0000000000, +-0.0022143788, +0.0000000000, +-0.0019788242, +0.0000000000, +-0.0017619530, +0.0000000000, +-0.0015626849, +0.0000000000, +-0.0013800042, +0.0000000000, +-0.0012129516, +0.0000000000, +-0.0010606179, +0.0000000000, +-0.0009221384, +0.0000000000, +-0.0007966882, +0.0000000000, +-0.0006834783, +0.0000000000, +-0.0005817526, +0.0000000000, +-0.0004907854, +0.0000000000, +-0.0004098791, +0.0000000000, +-0.0003383626, +0.0000000000, +-0.0002755900, +0.0000000000, +-0.0002209394, +0.0000000000, +-0.0001738121, +0.0000000000, +-0.0001336321, +0.0000000000, +-0.0000998452, +0.0000000000, +-0.0000719194, +0.0000000000, +-0.0000493437, +0.0000000000, +-0.0000316288, +0.0000000000, +-0.0000183062, +0.0000000000, +-0.0000089287, +0.0000000000, +-0.0000030703, +0.0000000000, +-0.0000003256}; diff --git a/examples/TestEqualizer1/TestEqualizer1.ino b/examples/TestEqualizer1/TestEqualizer1.ino new file mode 100644 index 0000000..a08f2f0 --- /dev/null +++ b/examples/TestEqualizer1/TestEqualizer1.ino @@ -0,0 +1,81 @@ +/* TestEqualizer1.ino Bob Larkin 10 May 2020 + * This is a test of the Filter Equalizer for Teensy Audio. + * This version is for the Chip Audette _F32 Library. + */ + +#include "Audio.h" +#include "OpenAudio_ArduinoLibrary.h" +#include "DSP_TeensyAudio_F32.h" + +AudioInputI2S i2s1; // Start interrupts +AudioSynthWaveformSine_F32 sine1; // Test signal +AudioFilterEqualizer_F32 equalize1; +AudioRecordQueue_F32 queue1; // The LSB output +AudioConnection_F32 patchCord1(sine1, 0, equalize1, 0); +AudioConnection_F32 patchCord2(equalize1, 0, queue1, 0); + +//nBands = 10 This octave band equalizer is set strangely to demonstrate the Equalizer +float fBand[] = { 40.0, 80.0, 160.0, 320.0, 640.0, 1280.0, 2560.0, 5120.0, 10240.0, 22058.5}; +float dbBand[] = {10.0, 2.0, -2.0, -5.0, -2.0, -4.0, -20.0, 6.0, 10.0, -100}; +float equalizeCoeffs[249]; +float dt1[128]; +float *pq1, *pd1; +int i, k; +float32_t dBResponse[500]; // Show lots of detail + +void setup(void) { + AudioMemory(5); + AudioMemory_F32(10); + Serial.begin(300); delay(1000); + Serial.println("*** Test Audio Equalizer ***"); + sine1.frequency(1000.0f); + + // Initialize the equalizer with 10 bands, 199 FIR coefficients and -65 dB sidelobes + uint16_t eq = equalize1.equalizerNew(10, &fBand[0], &dbBand[0], 199, &equalizeCoeffs[0], 65); + if (eq == ERR_EQ_BANDS) + Serial.println("Equalizer failed: Invalid number of frequency bands."); + else if (eq == ERR_EQ_SIDELOBES) + Serial.println("Equalizer failed: Invalid sidelobe specification."); + else if (eq == ERR_EQ_NFIR) + Serial.println("Equalizer failed: Invalid number of FIR coefficients."); + else + Serial.println("Equalizer initialized successfully."); + + // Get frequency response in dB for 500 points, uniformly spaced from 0 to 21058 Hz + equalize1.getResponse(500, dBResponse); + Serial.println("Freq Hz, Response dB"); + for(int m=0; m<500; m++) { // Print the response in Hz, dB, suitable for a spread sheet + Serial.print((float32_t)m * 22058.5f / 500.0f); + Serial.print(","); + Serial.println(dBResponse[m]); + } + i = -10; k=0; +} + +void loop(void) { + if(i<0) i++; // Get past startup + if(i==0) queue1.begin(); + + // Collect 128 samples and output to Serial + // This "if" will be active for i == 0 + if (queue1.available() >= 1 && i >= 0) { + pq1 = queue1.readBuffer(); + pd1 = &dt1[0]; + for(k=0; k<128; k++) { + *pd1++ = *pq1++; + } + i=1; // Only collect 1 block + Serial.print("Maximum F32 memory usage at loop:"); + Serial.println( AudioMemoryUsageMax_F32() ); + queue1.freeBuffer(); + queue1.end(); // No more data to queue1 + } + if(i == 1) { + // Uncomment the next 3 lines to printout a sample of the sine wave. + /*Serial.println("128 Samples: "); + for (k=0; k<128; k++) { + Serial.println (dt1[k],7); */ + } + i = 2; + delay(2); +} diff --git a/examples/TestEqualizer1/hilbert121A.h b/examples/TestEqualizer1/hilbert121A.h new file mode 100644 index 0000000..d624d04 --- /dev/null +++ b/examples/TestEqualizer1/hilbert121A.h @@ -0,0 +1,123 @@ +// Following is 121 term Hilbert FIR filter +float32_t hilbert121A[121] = { + 0.000000000000000000, + 0.000773378567767513, + 0.000000000000000000, + 0.001046207887980644, + 0.000000000000000000, + 0.001368896533613985, + 0.000000000000000000, + 0.001746769975247667, + 0.000000000000000000, + 0.002185555845922462, + 0.000000000000000000, + 0.002691457154069645, + 0.000000000000000000, + 0.003271251311125927, + 0.000000000000000000, + 0.003932423233774751, + 0.000000000000000000, + 0.004683343721596901, + 0.000000000000000000, + 0.005533508538632429, + 0.000000000000000000, + 0.006493859804516438, + 0.000000000000000000, + 0.007577220484233372, + 0.000000000000000000, + 0.008798886675905997, + 0.000000000000000000, + 0.010177443901536392, + 0.000000000000000000, + 0.011735907609641917, + 0.000000000000000000, + 0.013503343224246872, + 0.000000000000000000, + 0.015517212970554440, + 0.000000000000000000, + 0.017826854793349920, + 0.000000000000000000, + 0.020498780519188083, + 0.000000000000000000, + 0.023625003856774591, + 0.000000000000000000, + 0.027336628208641155, + 0.000000000000000000, + 0.031827023036304102, + 0.000000000000000000, + 0.037393534868609392, + 0.000000000000000000, + 0.044517689704988733, + 0.000000000000000000, + 0.054032871748808158, + 0.000000000000000000, + 0.067515548043274365, + 0.000000000000000000, + 0.088347125250410385, + 0.000000000000000000, + 0.125324201622410869, + 0.000000000000000000, + 0.210709715079613419, + 0.000000000000000000, + 0.634897508268964295, + 0.000000000000000000, +-0.634897508268964295, + 0.000000000000000000, +-0.210709715079613419, + 0.000000000000000000, +-0.125324201622410869, + 0.000000000000000000, +-0.088347125250410385, + 0.000000000000000000, +-0.067515548043274365, + 0.000000000000000000, +-0.054032871748808158, + 0.000000000000000000, +-0.044517689704988733, + 0.000000000000000000, +-0.037393534868609392, + 0.000000000000000000, +-0.031827023036304102, + 0.000000000000000000, +-0.027336628208641155, + 0.000000000000000000, +-0.023625003856774591, + 0.000000000000000000, +-0.020498780519188083, + 0.000000000000000000, +-0.017826854793349920, + 0.000000000000000000, +-0.015517212970554440, + 0.000000000000000000, +-0.013503343224246872, + 0.000000000000000000, +-0.011735907609641917, + 0.000000000000000000, +-0.010177443901536392, + 0.000000000000000000, +-0.008798886675905997, + 0.000000000000000000, +-0.007577220484233372, + 0.000000000000000000, +-0.006493859804516438, + 0.000000000000000000, +-0.005533508538632429, + 0.000000000000000000, +-0.004683343721596901, + 0.000000000000000000, +-0.003932423233774751, + 0.000000000000000000, +-0.003271251311125927, + 0.000000000000000000, +-0.002691457154069645, + 0.000000000000000000, +-0.002185555845922462, + 0.000000000000000000, +-0.001746769975247667, + 0.000000000000000000, +-0.001368896533613985, + 0.000000000000000000, +-0.001046207887980644, + 0.000000000000000000, +-0.000773378567767513, + 0.000000000000000000}; diff --git a/examples/TestEqualizer1/hilbert19A.h b/examples/TestEqualizer1/hilbert19A.h new file mode 100644 index 0000000..155b132 --- /dev/null +++ b/examples/TestEqualizer1/hilbert19A.h @@ -0,0 +1,21 @@ +// Following is a 19 term Hilbert FIR filter +float hilbert19A[]={ + 0.003006666677728199, + 0.000000000000000000, + 0.017439390086760032, + 0.000000000000000000, + 0.058235158985840196, + 0.000000000000000000, + 0.161905323559397157, + 0.000000000000000000, + 0.617838316334015092, + 0.000000000000000000, +-0.617838316334015092, + 0.000000000000000000, +-0.161905323559397157, + 0.000000000000000000, +-0.058235158985840196, + 0.000000000000000000, +-0.017439390086760032, + 0.000000000000000000, +-0.003006666677728199}; diff --git a/examples/TestEqualizer1/hilbert251A.h b/examples/TestEqualizer1/hilbert251A.h new file mode 100644 index 0000000..745db86 --- /dev/null +++ b/examples/TestEqualizer1/hilbert251A.h @@ -0,0 +1,253 @@ +// Following is 251 term Hilbert FIR filter +float32_t hilbert251A[]={ +0.0000003255, +0.0000000000, +0.0000030702, +0.0000000000, +0.0000089286, +0.0000000000, +0.0000183061, +0.0000000000, +0.0000316287, +0.0000000000, +0.0000493436, +0.0000000000, +0.0000719193, +0.0000000000, +0.0000998451, +0.0000000000, +0.0001336320, +0.0000000000, +0.0001738120, +0.0000000000, +0.0002209393, +0.0000000000, +0.0002755899, +0.0000000000, +0.0003383625, +0.0000000000, +0.0004098790, +0.0000000000, +0.0004907853, +0.0000000000, +0.0005817525, +0.0000000000, +0.0006834782, +0.0000000000, +0.0007966881, +0.0000000000, +0.0009221383, +0.0000000000, +0.0010606178, +0.0000000000, +0.0012129515, +0.0000000000, +0.0013800041, +0.0000000000, +0.0015626848, +0.0000000000, +0.0017619529, +0.0000000000, +0.0019788241, +0.0000000000, +0.0022143787, +0.0000000000, +0.0024697715, +0.0000000000, +0.0027462425, +0.0000000000, +0.0030451312, +0.0000000000, +0.0033678928, +0.0000000000, +0.0037161183, +0.0000000000, +0.0040915578, +0.0000000000, +0.0044961498, +0.0000000000, +0.0049320558, +0.0000000000, +0.0054017033, +0.0000000000, +0.0059078375, +0.0000000000, +0.0064535860, +0.0000000000, +0.0070425380, +0.0000000000, +0.0076788436, +0.0000000000, +0.0083673390, +0.0000000000, +0.0091137048, +0.0000000000, +0.0099246683, +0.0000000000, +0.0108082660, +0.0000000000, +0.0117741868, +0.0000000000, +0.0128342256, +0.0000000000, +0.0140028938, +0.0000000000, +0.0152982506, +0.0000000000, +0.0167430570, +0.0000000000, +0.0183664064, +0.0000000000, +0.0202060801, +0.0000000000, +0.0223120327, +0.0000000000, +0.0247516963, +0.0000000000, +0.0276183140, +0.0000000000, +0.0310445375, +0.0000000000, +0.0352256211, +0.0000000000, +0.0404611696, +0.0000000000, +0.0472354231, +0.0000000000, +0.0563851215, +0.0000000000, +0.0694911881, +0.0000000000, +0.0899418673, +0.0000000000, +0.1265473875, +0.0000000000, +0.2116132716, +0.0000000000, +0.6358933477, +0.0000000000, +-0.6358933478, +0.0000000000, +-0.2116132717, +0.0000000000, +-0.1265473876, +0.0000000000, +-0.0899418674, +0.0000000000, +-0.0694911882, +0.0000000000, +-0.0563851216, +0.0000000000, +-0.0472354232, +0.0000000000, +-0.0404611697, +0.0000000000, +-0.0352256212, +0.0000000000, +-0.0310445376, +0.0000000000, +-0.0276183141, +0.0000000000, +-0.0247516964, +0.0000000000, +-0.0223120328, +0.0000000000, +-0.0202060802, +0.0000000000, +-0.0183664065, +0.0000000000, +-0.0167430571, +0.0000000000, +-0.0152982507, +0.0000000000, +-0.0140028939, +0.0000000000, +-0.0128342257, +0.0000000000, +-0.0117741869, +0.0000000000, +-0.0108082661, +0.0000000000, +-0.0099246684, +0.0000000000, +-0.0091137049, +0.0000000000, +-0.0083673391, +0.0000000000, +-0.0076788437, +0.0000000000, +-0.0070425381, +0.0000000000, +-0.0064535861, +0.0000000000, +-0.0059078376, +0.0000000000, +-0.0054017034, +0.0000000000, +-0.0049320559, +0.0000000000, +-0.0044961499, +0.0000000000, +-0.0040915579, +0.0000000000, +-0.0037161184, +0.0000000000, +-0.0033678929, +0.0000000000, +-0.0030451313, +0.0000000000, +-0.0027462426, +0.0000000000, +-0.0024697716, +0.0000000000, +-0.0022143788, +0.0000000000, +-0.0019788242, +0.0000000000, +-0.0017619530, +0.0000000000, +-0.0015626849, +0.0000000000, +-0.0013800042, +0.0000000000, +-0.0012129516, +0.0000000000, +-0.0010606179, +0.0000000000, +-0.0009221384, +0.0000000000, +-0.0007966882, +0.0000000000, +-0.0006834783, +0.0000000000, +-0.0005817526, +0.0000000000, +-0.0004907854, +0.0000000000, +-0.0004098791, +0.0000000000, +-0.0003383626, +0.0000000000, +-0.0002755900, +0.0000000000, +-0.0002209394, +0.0000000000, +-0.0001738121, +0.0000000000, +-0.0001336321, +0.0000000000, +-0.0000998452, +0.0000000000, +-0.0000719194, +0.0000000000, +-0.0000493437, +0.0000000000, +-0.0000316288, +0.0000000000, +-0.0000183062, +0.0000000000, +-0.0000089287, +0.0000000000, +-0.0000030703, +0.0000000000, +-0.0000003256}; diff --git a/examples/TestEqualizer1Audio/TestEqualizer1Audio.ino b/examples/TestEqualizer1Audio/TestEqualizer1Audio.ino new file mode 100644 index 0000000..f01a285 --- /dev/null +++ b/examples/TestEqualizer1Audio/TestEqualizer1Audio.ino @@ -0,0 +1,105 @@ +/* TestEqualizer1Audio.ino Bob Larkin 18 May 2020 + * This is a test of the Filter Equalizer using the Audette OpenAudio_ArduinoLibrary. + * Runs two different equalizers, switching every 5 seconds to demonstrate + * dynamic filter changing. + */ + +#include "Audio.h" +#include "OpenAudio_ArduinoLibrary.h" +#include "DSP_TeensyAudio_F32.h" +#include + +// Signals from ADC to Equalizer to DAC using Teensy Audio Adaptor +AudioInputI2S i2sIn; +AudioConvert_I16toF32 convertItoF1; +AudioConvert_I16toF32 convertItoF2; +AudioFilterEqualizer_F32 equalize1; +AudioFilterEqualizer_F32 equalize2; +AudioConvert_F32toI16 convertFtoI1; +AudioConvert_F32toI16 convertFtoI2; +AudioOutputI2S i2sOut; +AudioConnection patchCord1(i2sIn, 0, convertItoF1, 0); +AudioConnection patchCord3(i2sIn, 1, convertItoF2, 0); +AudioConnection_F32 patchCord5(convertItoF1, 0, equalize1, 0); +AudioConnection_F32 patchCord7(convertItoF2, 0, equalize2, 0); +AudioConnection_F32 patchCord9(equalize1, 0, convertFtoI1, 0); +AudioConnection_F32 patchCordB(equalize2, 0, convertFtoI2, 0); +AudioConnection patchCordD(convertFtoI1, 0, i2sOut, 0); +AudioConnection patchCordF(convertFtoI2, 0, i2sOut, 1); +AudioControlSGTL5000 codec; + +// Some sort of octave band equalizer as a one alternative, 10 bands +float32_t fBand1[] = { 40.0, 80.0, 160.0, 320.0, 640.0, 1280.0, 2560.0, 5120.0, 10240.0, 22058.5}; +float32_t dbBand1[] = {10.0, 2.0, -2.0, -5.0, -2.0, -4.0, -10.0, 6.0, 10.0, -100}; + +// To contrast, put a strange bandpass filter as an alternative, 5 bands +float32_t fBand2[] = { 300.0, 500.0, 800.0, 1000.0, 22058.5}; +float32_t dbBand2[] = {-60.0, 0.0, -20.0, 0.0, -60.0}; + +float32_t equalizeCoeffs[200]; +int16_t k = 0; + +void setup(void) { + AudioMemory(5); + AudioMemory_F32(10); + Serial.begin(300); delay(1000); + Serial.println("*** Test Audio Equalizer ***"); + codec.enable(); + codec.inputSelect(AUDIO_INPUT_LINEIN); + codec.adcHighPassFilterDisable(); // necessary to suppress noise + codec.lineInLevel(2); + codec.volume(0.8); + + // Initialize the equalizer with 10 bands, fBand1[] 199 FIR coefficients + // -65 dB sidelobes, 16384 Max coefficient value + uint16_t eq = equalize1.equalizerNew(10, &fBand1[0], &dbBand1[0], 30, &equalizeCoeffs[0], 60.0); + if (eq == ERR_EQ_BANDS) + Serial.println("Equalizer failed: Invalid number of frequency bands."); + else if (eq == ERR_EQ_SIDELOBES) + Serial.println("Equalizer failed: Invalid sidelobe specification."); + else if (eq == ERR_EQ_NFIR) + Serial.println("Equalizer failed: Invalid number of FIR coefficients."); + else + Serial.println("Equalizer initialized successfully."); + + eq = equalize2.equalizerNew(10, &fBand1[0], &dbBand1[0], 30, &equalizeCoeffs[0], 60.0); + if (eq == ERR_EQ_BANDS) + Serial.println("Equalizer failed: Invalid number of frequency bands."); + else if (eq == ERR_EQ_SIDELOBES) + Serial.println("Equalizer failed: Invalid sidelobe specification."); + else if (eq == ERR_EQ_NFIR) + Serial.println("Equalizer failed: Invalid number of FIR coefficients."); + else + Serial.println("Equalizer initialized successfully."); + // for (k=0; k<200; k++) Serial.println(equalizeCoeffs[k]); +} + +void loop(void) { +#if 0 // Change between two filters every 10 seconds. + + // To run with just the 10-band equalizer, comment out the entire loop with "/* ... */" + + if(k == 0) { + k = 1; + equalize1.equalizerNew(10, &fBand1[0], &dbBand1[0], 200, &equalizeCoeffs[0], 65.0); + equalize2.equalizerNew(10, &fBand1[0], &dbBand1[0], 200, &equalizeCoeffs[0], 65.0); + } + else { // Swap back and forth + k = 0; + equalize1.equalizerNew(5, &fBand2[0], &dbBand2[0], 100, &equalizeCoeffs[0], 40.0); + equalize2.equalizerNew(5, &fBand2[0], &dbBand2[0], 100, &equalizeCoeffs[0], 40.0); + } +#endif + delay(10000); // Interrupts will keep the audio going + Serial.print("Proc = "); + Serial.print(AudioProcessorUsage()); + Serial.print(" ("); + Serial.print(AudioProcessorUsageMax()); + Serial.print("), Mem = "); + Serial.print(AudioMemoryUsage()); + Serial.print(" ("); + Serial.print(AudioMemoryUsageMax()); + Serial.println(")"); + AudioProcessorUsageMaxReset(); + AudioMemoryUsageMaxReset(); +} diff --git a/examples/TestFIRGeneral3/TestFIRGeneral3.ino b/examples/TestFIRGeneral3/TestFIRGeneral3.ino new file mode 100644 index 0000000..1389dfb --- /dev/null +++ b/examples/TestFIRGeneral3/TestFIRGeneral3.ino @@ -0,0 +1,88 @@ +/* TestEqualizer3General.ino Bob Larkin 22 May 2020 + * This is a test of the Filter Equalizer for Teensy Audio. + * This version is for the Chip Audette _F32 Library. + * See TestEqualizer2.ino for the int16 version. + */ + +#include "Audio.h" +#include "OpenAudio_ArduinoLibrary.h" +#include "DSP_TeensyAudio_F32.h" + +AudioInputI2S i2s1; +AudioSynthWaveformSine_F32 sine1; // Test signal +AudioFilterFIRGeneral_F32 firg1; +AudioRecordQueue_F32 queue1; // The LSB output +AudioConnection_F32 patchCord1(sine1, 0, firg1, 0); +AudioConnection_F32 patchCord2(firg1, 0, queue1, 0); + +float dbA[51]; +float equalizeCoeffs[51]; +float dt1[128]; +float *pq1, *pd1; +int i, k; +float32_t dBResponse[500]; // Show lots of detail + +void setup(void) { + AudioMemory(5); + AudioMemory_F32(10); + Serial.begin(300); delay(1000); + Serial.println("*** Test FIR General routine for F32 ***"); + sine1.frequency(1000.0f); + + // This simple case corresponds to Burris and Parks, "Digital Filter Design." + // Example 3.1when sidelobes are set to 0.0. + for (i=0; i<6; i++) dbA[i] = 0.0f; + for (i=6; i<11; i++) dbA[i] = -120.0f; + + // Initialize the FIR with 21 points, 21 FIR coefficients and -0 dB sidelobes + uint16_t eq = firg1.FIRGeneralNew(&dbA[0], 21, &equalizeCoeffs[0], 0.0); + if (eq == ERR_EQ_BANDS) + Serial.println("FIR General failed: Invalid number of frequency points."); + else if (eq == ERR_EQ_SIDELOBES) + Serial.println("FIR General failed: Invalid sidelobe specification."); + else if (eq == ERR_EQ_NFIR) + Serial.println("FIR General failed: Invalid number of FIR coefficients."); + else + Serial.println("FIR General initialized successfully."); + + Serial.println("FIR Coeffs"); + for(i=0; i<21; i++) Serial.println(equalizeCoeffs[i], 7); + Serial.println(""); + + // Get frequency response in dB for 200 points, uniformly spaced from 0 to 21058 Hz + firg1.getResponse(200, dBResponse); + Serial.println("Freq Hz, Response dB"); + for(int m=0; m<200; m++) { // Print the response in Hz, dB, suitable for a spread sheet + Serial.print((float32_t)m * 22058.5f / 200.0f); + Serial.print(","); + Serial.println(dBResponse[m]); + } + + i = -10; k=0; +} + +void loop(void) { + if(i<0) i++; // Get past startup + if(i==0) queue1.begin(); + + // Collect 128 samples and output to Serial + // This "if" will be active for i == 0 + if (queue1.available() >= 1 && i >= 0) { + pq1 = queue1.readBuffer(); + pd1 = &dt1[0]; + for(k=0; k<128; k++) { + *pd1++ = *pq1++; + } + i=1; // Only collect 1 block + queue1.freeBuffer(); + queue1.end(); // No more data to queue1 + } + if(i == 1) { + // Uncomment the next 3 lines to printout a sample of the sine wave. + Serial.println("128 Samples: "); + for (k=0; k<128; k++) + Serial.println (dt1[k],7); + i = 2; + } + delay(2); +} diff --git a/examples/TestFIRGeneralLarge4/HP55ieeeRSL.gif b/examples/TestFIRGeneralLarge4/HP55ieeeRSL.gif new file mode 100644 index 0000000..c91e64d Binary files /dev/null and b/examples/TestFIRGeneralLarge4/HP55ieeeRSL.gif differ diff --git a/examples/TestFIRGeneralLarge4/KaiserBetaSidelobes.gif b/examples/TestFIRGeneralLarge4/KaiserBetaSidelobes.gif new file mode 100644 index 0000000..8ff338b Binary files /dev/null and b/examples/TestFIRGeneralLarge4/KaiserBetaSidelobes.gif differ diff --git a/examples/TestFIRGeneralLarge4/Test4.gnumeric b/examples/TestFIRGeneralLarge4/Test4.gnumeric new file mode 100644 index 0000000..f4da4cf Binary files /dev/null and b/examples/TestFIRGeneralLarge4/Test4.gnumeric differ diff --git a/examples/TestFIRGeneralLarge4/TestFIRGeneralLarge4.ino b/examples/TestFIRGeneralLarge4/TestFIRGeneralLarge4.ino new file mode 100644 index 0000000..cffbbac --- /dev/null +++ b/examples/TestFIRGeneralLarge4/TestFIRGeneralLarge4.ino @@ -0,0 +1,127 @@ +/* TestFIRGeneralLarge4.ino Bob Larkin 24 May 2020 + * Test the generation of FIR filters and obtaining their + * responses. See TestFIRGeneralLarge5.ino to run the filters. + * This is a test of the AudioFilterFIR_F32 for Teensy Audio. + */ + +#include "Audio.h" +#include "OpenAudio_ArduinoLibrary.h" +#include "DSP_TeensyAudio_F32.h" + +AudioInputI2S i2s1; // Creates timing +AudioFilterFIRGeneral_F32 firg1; + +float32_t dbA[2000]; +float32_t workSpace[4128]; +float32_t equalizeCoeffs[4000]; +int i, k; +float32_t dBResponse[5000]; // Extreme size, can show lots of detail! + +// Following is based on Rabiner, McGonegal and Paul, FWFIR test case, but has +// had the windowing effect divided out to show their basic filter response. +float32_t hpCoeff[55]={ +-0.003643388f, -0.007195844f, 0.012732424f, -0.007796006f, -0.004276372f, + 0.013760401f, -0.012262941f, 0.000000304f, 0.013553423f, -0.016818445f, + 0.005786355f, 0.011693321f, -0.021220577f, 0.013364335f, 0.007566100f, +-0.025227439f, 0.023410856f, -0.000000305f, -0.028612852f, 0.037841329f, +-0.014052131f, -0.031182658f, 0.063661835f, -0.046774701f, -0.032787389f, + 0.151364865f, -0.257518316f, 0.300000700f, -0.257518316f, 0.151364865f, +-0.032787389f, -0.046774701f, 0.063661835f, -0.031182658f, -0.014052131f, + 0.037841329f, -0.028612852f, -0.000000305f, 0.023410856f, -0.025227439f, + 0.007566100f, 0.013364335f, -0.021220577f, 0.011693321f, 0.005786355f, +-0.016818445f, 0.013553423f, 0.000000304f, -0.012262941f, 0.013760401f, +-0.004276372f, -0.007796006f, 0.012732424f, -0.007195844f, -0.003643388}; + +void setup(void) { + uint16_t nFIR = 4; + float32_t sidelobes = 45; + AudioMemory(5); + AudioMemory_F32(10); + Serial.begin(300); delay(1000); + Serial.println("*** Test Response of FIR General routine for F32 ***"); + +// Un-comment one of the following, or use the direct Coefficient load that follows + +// Test case for odd nFIR from Burris and Parks Example 3.1 + nFIR = 21; sidelobes = 0.0; // No window + dbA[0]=0.0; dbA[1]=0.0; dbA[2]=0.0; dbA[3]=0.0; dbA[4]=0.0; dbA[5]=0.0; + dbA[6]=-140.0; dbA[7]=-140.0; dbA[8]=-140.0; dbA[9]=-140.0; dbA[10]=-140.0; + +// Test case for even nFIR from Burris and Parks Example 3.2 +// nFIR = 20; sidelobes = 0.0; +// dbA[0]=0.0; dbA[1]=0.0; dbA[2]=0.0; dbA[3]=0.0; dbA[4]=-120.0; dbA[5]=-120.0; +// dbA[6]=-120.0; dbA[7]=-120.0; dbA[8]=-120.0; dbA[9]=-120.0; + +// Test case for odd nFIR from Burris and Parks Example 3.1 with Kaiser Window +// nFIR = 21; sidelobes = 45.0f; +// dbA[0]=0.0; dbA[1]=0.0; dbA[2]=0.0; dbA[3]=0.0; dbA[4]=0.0; dbA[5]=0.0; +// dbA[6]=-140.0; dbA[7]=-140.0; dbA[8]=-140.0; dbA[9]=-140.0; dbA[10]=-140.0f; + +// Test case for even nFIR from Burris and Parks Example 3.2 +// nFIR = 20; sidelobes = 45.0f; +// dbA[0]=0.0; dbA[1]=0.0; dbA[2]=0.0; dbA[3]=0.0; dbA[4]=0.0; dbA[5]=-120.0; +// dbA[6]=-120.0; dbA[7]=-120.0; dbA[8]=-120.0; dbA[9]=-120.0; + +// Rabiner, McGonegal and Paul, FWFIR test case, HP Kaiser SL=60, fco=0.35 +// nFIR = 55; sidelobes = 60.0; +// for(i=0; i<28; i++) { +// if (i<20) dbA[i] = -120.0f; +// else dbA[i] = 0.0f; +// } + +// The next 3 examples are tests of extremely large FIR arrays and only for Teensy 4 +// Generate a response curve for our "octave" equalizer, tried with 30 dB sidelobes +// nFIR = 3999; sidelobes = 30.0f; +// dbA[0] = 12.0f; // Avoid log(0) +// for(i=1; i<2000; i++) +// dbA[i] = 12.0f*cosf(10.0f*log10f(11.02925f*(float32_t)i)); + +// Generate a slowly changing response, up 12 dB at bass and treble and -12 dB +// for mid-frequencies. Scaled to "octaves" by the form cos(log(f)) +// nFIR = 3999; sidelobes = 30.0f; +// dbA[0] = 12.0f; +// dbA[1] = 12.0f; +// for(i=2; i<2000; i++) +// dbA[i] = -12.0f*cosf(2.2f*log10f(11.02925f*(float32_t)i)); + +// Narrow band "CW" filter, tried with 60 dB sidelobes (60 to 80 or 70 to 72 +// nFIR = 3999; sidelobes = 60; +// for(i=0; i<2000; i++) { +// if (i<70 || i>72) dbA[i] = -120.0f; +// else dbA[i] = 0.0f; +// } + + // Initialize the FIR filter design and run + uint16_t eq = firg1.FIRGeneralNew(&dbA[0], nFIR, &equalizeCoeffs[0], sidelobes, &workSpace[0]); + if (eq == ERR_EQ_BANDS) + Serial.println("FIR General failed: Invalid number of frequency points."); + else if (eq == ERR_EQ_SIDELOBES) + Serial.println("FIR General failed: Invalid sidelobe specification."); + else if (eq == ERR_EQ_NFIR) + Serial.println("FIR General failed: Invalid number of FIR coefficients."); + else + Serial.println("FIR General initialized successfully."); + +// Load (if not commented) non-windowed HP filter coefficients, test case: +// firg1.LoadCoeffs(55, hpCoeff, workSpace); // Direct load of FIR coefficients, already calculated + + Serial.println("FIR Coeffs"); + for(i=0; i72) dbA[i] = -120.0f; + else dbA[i] = 0.0f; + } + + // Initialize the FIR filter design and run + uint16_t eq = firg1.FIRGeneralNew(&dbA[0], nFIR, &equalizeCoeffs[0], sidelobes, &workSpace[0]); + if (eq == ERR_EQ_BANDS) + Serial.println("FIR General failed: Invalid number of frequency points."); + else if (eq == ERR_EQ_SIDELOBES) + Serial.println("FIR General failed: Invalid sidelobe specification."); + else if (eq == ERR_EQ_NFIR) + Serial.println("FIR General failed: Invalid number of FIR coefficients."); + else + Serial.println("FIR General initialized successfully."); +} + +void loop(void) { +} diff --git a/radioNoiseBlanker_F32.cpp b/radioNoiseBlanker_F32.cpp new file mode 100644 index 0000000..f415e64 --- /dev/null +++ b/radioNoiseBlanker_F32.cpp @@ -0,0 +1,90 @@ +/* + * radioNoiseBlanker_F32.cpp + * + * 22 March 2020 + * Bob Larkin, in support of the library: + * Chip Audette, OpenAudio, Apr 2017 + * ------------------- + * + * MIT License, Use at your own risk. +*/ + +#include "radioNoiseBlanker_F32.h" + +void radioNoiseBlanker_F32::update(void) { + audio_block_f32_t *blockIn, *blockOut=NULL; + uint16_t i; + float32_t absSignal; + + // Get input block // <data + // based on whether gate is open or not. + for(i=0; idata[i]; // ith data + delayData[(i+in_index) & delayBufferMask] = datai; // Put ith data to circular delay buffer + + absSignal = fabsf(datai); // Rectified I-F + runningSum += fabsf(datai); // Update by adding one rectified point + runningSum -= delayData[(i + in_index - RUNNING_SUM_SIZE) & delayBufferMask]; // & subtract one + + pulseTime++; // This keeps track of leading and trailing delays of the gate pulse + if (absSignal > (threshold * runningSum)) { // A noise pulse event + if (gateOn == true) { // Signals are flowing, this is beginning of noise pulse event + gateOn = false; + pulseTime = -nAnticipation; + } + else { // gateOn==false, we are already in a noise pulse event + if (pulseTime > 0) { // Waiting for pulse event to end + pulseTime = 0; // Keep waiting + } + } + } + else { // Noise pulse is below threshold + if (gateOn == true) { // Signals are flowing normally + pulseTime = -9999; + } + else { // gateOn==false, we are already in a noise pulse event + if(pulseTime >= nDecay) { // End of a pulse event, turn gate on + gateOn = true; + pulseTime = -9999; + } + } + } + // Ready to enter I-F data to output, offset in time by "nAnticipation" + if (pulseTime == -9999) + blockOut->data[i] = delayData[(256 + i - nAnticipation) & delayBufferMask]; // Need 256?? + else // -nAnticipation < pulseTime < nDecay i.e., blanked out + blockOut->data[i] = 0.0f; + } // End of loop point by point over input 128 data points + + AudioStream_F32::release (blockIn); + AudioStream_F32::transmit(blockOut, 0); // send the delayed or blanked data + AudioStream_F32::release (blockOut); + + // Update pointer in_index to delay line for next 128 update + in_index = (in_index + block_size) & delayBufferMask; +} + diff --git a/radioNoiseBlanker_F32.h b/radioNoiseBlanker_F32.h new file mode 100644 index 0000000..b0a0fbf --- /dev/null +++ b/radioNoiseBlanker_F32.h @@ -0,0 +1,126 @@ +/* radioNoiseBlanker_F32.h + * 15 May 2020 Bob Larkin + * + * Thanks to PJRC and Pau Stoffregen for the Teensy processor, Teensyduino + * and Teensy Audio. Thanks to Chip Audette for the F32 extension + * work. Alll of that makes this possible. Bob + * + * 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, development funding 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. + */ + +/* + * An I-F signal comes in with occassional noise pulses. This block watches for + * the pulses and turns off the I-F while the pulse exists. In order to + * be as smart as possible this looks ahead by a number of samples, nAnticipation. + * Likewise, the I-F is left off for nDecay samples after the pulse ends. + * Various methods could be to be used to "turn off" the i-f, + * including replacement with waveforms. + * As of this initial write, zeros are used in the waveform. + * + * A threshold can be adjusted via setNB(). This is compared with the + * average rectified voltage being received. If this is too small, like 1.5 + * or 2.0, we will be loosing data by blanking. If we set it too high, like 20.0, + * we will not blank noise pulses. Experiments will find a god setting. + * With a sine wave input and no impulse noise, the average rectified signal + * should be about 0.637. To catch the top of that requires a threshold of + * 1/0.637 = 1.57. That would seem to be a very minimal threshold setting. + * + * Status: Tested on computer data + * Functions: + * setNB( ... ); // See below + * showError(e); // e=1 means to print update() error messages + * enable(runNB); // When runNB=false the signals bypasses the update() processing + * Examples: + * TestNoiseBlanker1.ino 128 data for plotting and examination + * Time: Update() of 128 samples 32 microseconds + */ + +#ifndef _radio_noise_blanker_f32_h +#define _radio_noise_blanker_f32_h + +#include "AudioStream_F32.h" +#include "arm_math.h" + +#define RUNNING_SUM_SIZE 125 + +class radioNoiseBlanker_F32 : public AudioStream_F32 { +//GUI: inputs:1, outputs:1 //this line used for automatic generation of GUI node +//GUI: shortName: NoiseBlanker +public: + // Option of AudioSettings_F32 change to block size (no sample rate dependent variables here): + radioNoiseBlanker_F32(void) : AudioStream_F32(1, inputQueueArray_f32) { + block_size = AUDIO_BLOCK_SAMPLES; + } + radioNoiseBlanker_F32(const AudioSettings_F32 &settings) : AudioStream_F32(1, inputQueueArray_f32) { + block_size = settings.audio_block_samples; + } + + void enable(bool _runNB) { + runNB = _runNB; + } + + void setNoiseBlanker(float32_t _threshold, uint16_t _nAnticipation, uint16_t _nDecay) { + if (_threshold < 0.0) threshold = 0.0; + else threshold = _threshold/(float32_t)RUNNING_SUM_SIZE; + + if(_nAnticipation < 1) nAnticipation = 1; + else if (_nAnticipation >125) nAnticipation = 125; + else nAnticipation = _nAnticipation; + + if (_nDecay < 1) nDecay = 1; + else if (_nDecay > 10) nDecay = 10; // Is this handled OK at end of block? + else nDecay = _nDecay; + } + + void showError(uint16_t e) { // This is for debug + errorPrint = e; + } + + void update(void); + +private: + uint16_t block_size = AUDIO_BLOCK_SAMPLES; + // Input data pointers + audio_block_f32_t *inputQueueArray_f32[1]; + + // Control error printing in update() 0=No print + uint16_t errorPrint = 0; + + // To look ahead we need a delay of up to 128 samples. Too much removes good data. + // Too little enters noise pulse data to the output. This can be a simple circular + // buffer if we make the buffer a power of 2 in length and binary-truncate the index. + // Choose 2^8 = 256. + float32_t delayData[256]; // The circular delay line + uint16_t in_index = 0; // Pointer to next block update entry + // And a mask to make the circular buffer limit to a power of 2 + uint16_t delayBufferMask = 0X00FF; + + // Three variables to allow .INO control of Noise Blanker + float32_t threshold = 1.0E6f; // Start disabled + uint16_t nAnticipation = 5; + uint16_t nDecay = 8; + + int16_t pulseTime = -9999; // Tracks nAnticipation and nDecay + bool gateOn = true; // Signals going through NB + float32_t runningSum = 0.0; + bool runNB = false; +}; +#endif diff --git a/readme.md b/readme.md index 1daf09b..55cd819 100644 --- a/readme.md +++ b/readme.md @@ -25,6 +25,12 @@ OpenAudio Library for Teensy 20-Moved in radio function RadioIQMixer_F32. 21-Added Example INO, FineFreqShift_OA.ino with Hilbert shifter and delay stereo. 22-Brought in support stuff mathDSP_F32 and sinTable512_F32.h. +23-Brought in AudioAnalyzePhase_F32 and Example +24-Brought in AudioFilterEqualyzer_F32 and 2 Example +25-Brought in AudioFilterFIRGeneral_F32 and 2 Example INOs +26-Brought in RadioFMDetector_F32 and Example +27-Brought in synth_sin_cos_F32 and test example +28-Brought in RadioNoiseBlanker_F32 and Example **Purpose**: The purpose of this library is to build upon the [Teensy Audio Library](http://www.pjrc.com/teensy/td_libs_Audio.html) to enable new functionality for real-time audio processing. diff --git a/synth_sin_cos_f32.cpp b/synth_sin_cos_f32.cpp new file mode 100644 index 0000000..4d18ffa --- /dev/null +++ b/synth_sin_cos_f32.cpp @@ -0,0 +1,103 @@ +/* synth_sin_cos_f32.cpp + * + * SynthSinCos_F32 Bob Larkin April 17, 2020 + * + * Based on Chip Audette's OpenAudio sine(), that was + * Modeled on: AudioSynthWaveformSine from Teensy Audio Library + * + * Purpose: Create sine and cosine wave of given amplitude, frequency + * and phase. Outputs in float32_t floating point. + * Routines are from the arm CMSIS library and use a 512 point lookup + * table with linear interpolation to achieve float accuracy limits. + * + * 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. + */ + +#include "synth_sin_cos_f32.h" + +// 513 values of the sine wave in a float array: +#include "sinTable512_f32.h" + +void AudioSynthSineCosine_F32::update(void) { + audio_block_f32_t *blockS, *blockC; + uint16_t index, i; + float32_t a, b, deltaPhase, phaseC; + + blockS = AudioStream_F32::allocate_f32(); // Output blocks + if (!blockS) return; + + blockC = AudioStream_F32::allocate_f32(); + if (!blockC) { + AudioStream_F32::release(blockS); + return; + } + + // doSimple has amplitude (-1, 1) and sin/cos differ by 90.00 degrees. + if (doSimple) { + for (i=0; i < block_length; i++) { + phaseS += phaseIncrement; + if (phaseS > 512.0f) + phaseS -= 512.0f; + index = (uint16_t) phaseS; + deltaPhase = phaseS -(float32_t) index; + /* Read two nearest values of input value from the sin table */ + a = sinTable512_f32[index]; + b = sinTable512_f32[index+1]; + blockS->data[i] = a + 0.001953125*(b-a)*deltaPhase; /* Linear interpolation process */ + + /* Repeat for cosine by adding 90 degrees phase */ + index = (index + 128) & 0x01ff; + /* Read two nearest values of input value from the sin table */ + a = sinTable512_f32[index]; + b = sinTable512_f32[index+1]; + /* deltaPhase will be the same as used for sin */ + blockC->data[i] = a + 0.001953125*(b-a)*deltaPhase; /* Linear interpolation process */ + } + } + else { // Do a more flexible update, i.e., not doSimple + for (i=0; i < block_length; i++) { + phaseS += phaseIncrement; + if (phaseS > 512.0f) phaseS -= 512.0f; + index = (uint16_t) phaseS; + deltaPhase = phaseS -(float32_t) index; + /* Read two nearest values of input value from the sin table */ + a = sinTable512_f32[index]; + b = sinTable512_f32[index+1]; + blockS->data[i] = amplitude_pk * (a + 0.001953125*(b-a)*deltaPhase); /* Linear interpolation process */ + + /* Shift forward phaseS_C and get cos. First, the calculation of index of the table */ + phaseC = phaseS + phaseS_C; + if (phaseC > 512.0f) phaseC -= 512.0f; + index = (uint16_t) phaseC; + deltaPhase = phaseC -(float32_t) index; + /* Read two nearest values of input value from the sin table */ + a = sinTable512_f32[index]; + b = sinTable512_f32[index+1]; + blockC->data[i] = amplitude_pk * (a + 0.001953125*(b-a)*deltaPhase); /* Linear interpolation process */ + } + } + AudioStream_F32::transmit(blockS, 0); + AudioStream_F32::release (blockS); + AudioStream_F32::transmit(blockC, 1); + AudioStream_F32::release (blockC); + return; +} + diff --git a/synth_sin_cos_f32.h b/synth_sin_cos_f32.h new file mode 100644 index 0000000..77b50aa --- /dev/null +++ b/synth_sin_cos_f32.h @@ -0,0 +1,162 @@ +/* synth_sin_cos_f32.h + * AudioSynthSineCosine_F32 + * + * Status: Checked for function and accuracy. 19 April 2020 + * + * Created: Bob Larkin 15 April 2020 + * + * Based on Chip Audette's OpenAudio sine(), that was + * Modeled on: AudioSynthWaveformSine from Teensy Audio Library + * + * Purpose: Create sine and cosine wave of given amplitude, frequency + * and phase. Outputs are audio_block_f32_t blocks of float32_t. + * Routines are from the arm CMSIS library and use a 512 point lookup + * table with linear interpolation to achieve float accuracy limits. + * + * This provides for setting the phase of the sine, setting the difference + * in phase between the sine and cosine and setting the + * (-amplitude, amplitude) range. If these are at the default values, + * called doSimple, the caluclation is faster. + * For doSimple either true or false, the frequency can be changed + * at will using the frequency() method. + * + * Defaults: + * Frequency: 1000.0 Hz + * Phase of Sine: 0.0 radians (0.0 deg) + * Phase of Cosine: pi/2 radians (90.0 deg) ahead of Sine + * Amplitude: -1.0 to 1.0 + * + * Time: T3.6 update() block of 128 with doSimple is 36 microseconds + * Same using flexible doSimple=false is 49 microseconds + * T4.0 update() block of 128 with doSimple is 16 microseconds + * Same using flexible doSimple=false is 24 microseconds + * + * 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 synth_sin_cos_f32_h_ +#define synth_sin_cos_f32_h_ + +#include "AudioStream_F32.h" +#include "arm_math.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#ifndef M_PI_2 +#define M_PI_2 1.57079632679489661923 +#endif + +#ifndef M_TWOPI +#define M_TWOPI (M_PI * 2.0) +#endif + +#define MF2_PI 6.2831853f + +class AudioSynthSineCosine_F32 : public AudioStream_F32 { +//GUI: inputs:0, outputs:2 //this line used for automatic generation of GUI node +//GUI: shortName:SineCosine //this line used for automatic generation of GUI node +public: + AudioSynthSineCosine_F32(void) : AudioStream_F32(0, NULL) { } //uses default AUDIO_SAMPLE_RATE from AudioStream.h + + AudioSynthSineCosine_F32(const AudioSettings_F32 &settings) : AudioStream_F32(0, NULL) { + setSampleRate_Hz(settings.sample_rate_Hz); + setBlockLength(settings.audio_block_samples); + } + + void frequency(float32_t fr) { // Frequency in Hz + freq = fr; + if (freq < 0.0f) freq = 0.0f; + else if (freq > sample_rate_Hz/2.0f) freq = sample_rate_Hz/2.0f; + phaseIncrement = 512.0f * freq / sample_rate_Hz; + } + + /* Externally, phase comes in the range (0,2*M_PI) keeping with C math functions + * Internally, the full circle is represented as (0.0, 512.0). This is + * convenient for finding the entry to the sine table. + */ + void phase_r(float32_t a) { + while (a < 0.0f) a += MF2_PI; + while (a > MF2_PI) a -= MF2_PI; + phaseS = 512.0f * a / MF2_PI; + doSimple = false; + return; + } + + // phaseS_C_r is the number of radians that the cosine output leads the + // sine output. The default is M_PI_2 = pi/2 = 1.57079633 radians, + // corresponding to 90.00 degrees cosine leading sine. + void phaseS_C_r(float32_t a) { + while (a < 0.0f) a += MF2_PI; + while (a > MF2_PI) a -= MF2_PI; + // Internally a full circle is 512.00 of phase + phaseS_C = 512.0f * a / MF2_PI; + doSimple = false; + return; + } + + // The amplitude, a, is the peak, as in zero-to-peak. This produces outputs + // ranging from -a to +a. Both outputs are the same amplitude. + void amplitude(float32_t a) { + amplitude_pk = a; + doSimple = false; + return; + } + + // Speed up calculations by setting phaseS_C=90deg, amplitude=1 + // Note, s=true will override any setting of phaseS_C_r or amplitude. + void simple(bool s) { + doSimple = s; + if(doSimple) { + phaseS_C = 128.0f; + amplitude_pk = 1.0f; + } + return; + } + + void setSampleRate_Hz(float32_t fs_Hz) { + // Check freq range + if (freq > sample_rate_Hz/2.0f) freq = sample_rate_Hz/2.f; + // update phase increment for new frequency + phaseIncrement = 512.0f * freq / fs_Hz; + } + + void setBlockLength(uint16_t bl) { + if(bl > 128) bl = 128; + block_length = bl; + } + + virtual void update(void); + +private: + float32_t freq = 1000.0f; + float32_t phaseS = 0.0f; + float32_t phaseS_C = 128.00; + float32_t amplitude_pk = 1.0f; + float32_t sample_rate_Hz = AUDIO_SAMPLE_RATE; + float32_t phaseIncrement = 512.0f * freq /sample_rate_Hz; + uint16_t block_length = 128; + // if only freq() is used, the complexities of phase, phaseS_C, + // and amplitude are not used, speeding up the sin and cos: + bool doSimple = true; +}; +#endif