pull/11/head
parent
ae99656d29
commit
0af24d2bef
@ -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 <math.h> |
||||
#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; j<block0->length; 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; i<block_size; i++) { |
||||
block0->data[i] = mult1 * acosf(mult0 * block0->data[i]);
|
||||
} |
||||
} |
||||
else if ((pdConfig & ACOS_MASK) == 0b00100) { // Fast acos from polynomial
|
||||
for (uint i = 0; i<block_size; i++) { |
||||
block0->data[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; i<block_size; i++) { |
||||
block0->data[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()
|
@ -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 <math.h> |
||||
|
||||
class AudioAnalyzePhase_F32 : public AudioStream_F32 { |
||||
//GUI: inputs:2, outputs:1 //this line used for automatic generation of GUI node
|
||||
//GUI: shortName: AnalyzePhase
|
||||
public: |
||||
// Option of AudioSettings_F32 change to block size or sample rate:
|
||||
AudioAnalyzePhase_F32(void) : AudioStream_F32(2, inputQueueArray_f32) { // default block_size and sampleRate_Hz
|
||||
// Initialize BiQuad IIR instance (ARM DSP Math Library)
|
||||
arm_biquad_cascade_df1_init_f32(&iir_inst, N_STAGES, &iir_coeffs[0], &IIRStateF32[0]); |
||||
} |
||||
// Constructor including new block_size and/or sampleRate_Hz
|
||||
AudioAnalyzePhase_F32(const AudioSettings_F32 &settings) : AudioStream_F32(2, inputQueueArray_f32) { |
||||
block_size = settings.audio_block_samples; |
||||
sampleRate_Hz = settings.sample_rate_Hz;
|
||||
// Initialize BiQuad IIR instance (ARM DSP Math Library)
|
||||
arm_biquad_cascade_df1_init_f32(&iir_inst, N_STAGES, &iir_coeffs[0], &IIRStateF32[0]); |
||||
} |
||||
|
||||
// Set AnalyzePhaseConfig while leaving pdConfig as is
|
||||
void setAnalyzePhaseConfig(const uint16_t _LPType, float32_t *_pCoeffs, uint16_t _nCoeffs) { |
||||
setAnalyzePhaseConfig( _LPType, _pCoeffs, _nCoeffs, pdConfig); |
||||
} |
||||
// Set AnalyzePhaseConfig in full generality
|
||||
void setAnalyzePhaseConfig(const uint16_t _LPType, float32_t *_pCoeffs, uint16_t _nCoeffs, uint16_t _pdConfig) { |
||||
AudioNoInterrupts(); // No interrupts while changing parameters
|
||||
LPType = _LPType; |
||||
if (LPType == NO_LP_FILTER) { |
||||
//Serial.println("Advice: in AnalyzePhase, for NO_LP_FILTER the output contains 2nd harmonics");
|
||||
//Serial.println(" that need external filtering.");
|
||||
} |
||||
else if (LPType == IIR_LP_FILTER) { |
||||
if(_pCoeffs != NULL){ |
||||
pIirCoeffs = _pCoeffs; |
||||
nIirCoeffs = _nCoeffs; |
||||
} |
||||
if (nIirCoeffs != 20){ |
||||
//Serial.println("Error, in AnalyzePhase, for IIR_LP_FILTER there must be 20 coefficients.");
|
||||
nIirCoeffs = 20; |
||||
} |
||||
arm_biquad_cascade_df1_init_f32(&iir_inst, N_STAGES, pIirCoeffs, &IIRStateF32[0]); |
||||
}
|
||||
else if (LPType==FIR_LP_FILTER) { |
||||
if(_pCoeffs != NULL){ |
||||
pFirCoeffs = _pCoeffs; |
||||
nFirCoeffs = _nCoeffs; |
||||
} |
||||
if (nFirCoeffs<4 || nFirCoeffs>NFIR_MAX) { // Too many or too few
|
||||
//Serial.print("Error, in AnalyzePhase, for FIR_LP_FILTER there must be >4 and <=");
|
||||
//Serial.print(NFIR_MAX);
|
||||
//Serial.println(" coefficients.");
|
||||
//Serial.println(" Restoring default IIR Filter.");
|
||||
LPType = IIR_LP_FILTER; |
||||
pIirCoeffs = &iir_coeffs[0]; |
||||
nIirCoeffs = 20; // Number of coefficients 20
|
||||
pdConfig = 0b11100; |
||||
LPType = IIR_LP_FILTER; // Variables were set in setup() above
|
||||
} |
||||
else { //Acceptable number, so initialize it
|
||||
arm_fir_init_f32(&fir_inst, nFirCoeffs, pFirCoeffs, &FIRStateF32[0], block_size); |
||||
} |
||||
} |
||||
pdConfig = _pdConfig; |
||||
AudioInterrupts();
|
||||
} |
||||
|
||||
void showError(uint16_t e) { |
||||
errorPrint = e; |
||||
} |
||||
|
||||
void update(void); |
||||
|
||||
private: |
||||
float32_t sampleRate_Hz = AUDIO_SAMPLE_RATE_EXACT; |
||||
uint16_t block_size = AUDIO_BLOCK_SAMPLES; |
||||
// Two input data pointers
|
||||
audio_block_f32_t *inputQueueArray_f32[2]; |
||||
// Variables controlling the configuration
|
||||
uint16_t LPType = IIR_LP_FILTER; // NO_LP_FILTER, IIR_LP_FILTER or FIR_LP_FILTER
|
||||
float32_t *pIirCoeffs = &iir_coeffs[0]; // Coefficients for IIR
|
||||
float32_t *pFirCoeffs = &fir_coeffs[0]; // Coefficients for FIR
|
||||
uint16_t nIirCoeffs = 20; // Number of coefficients 20
|
||||
uint16_t nFirCoeffs = 53; // Number of coefficients <=200
|
||||
uint16_t pdConfig = 0b11100; // No limiter, fast acos, scale multiplier, radians out;
|
||||
// Control error printing in update(). Should never be enabled
|
||||
// until all audio objects have been initialized.
|
||||
// Only used as 0 or 1 now, but 16 bits are available.
|
||||
uint16_t errorPrint = 0; |
||||
|
||||
// *Temporary* - TEST_TIME allows measuring time in microseconds for each part of the update()
|
||||
#if TEST_TIME |
||||
elapsedMicros tElapse; |
||||
int32_t iitt = 998000; // count up to a million during startup
|
||||
#endif |
||||
|
||||
/* FIR filter designed with http://t-filter.appspot.com
|
||||
* Sampling frequency: 44100 Hz |
||||
* 0 Hz - 4000 Hz gain = 1.0, ripple = 0.101 dB |
||||
* 7000 - 22000 Hz attenuation >= 81.8 dB |
||||
* Suitable for measuring phase in communications systems with linear phase. |
||||
*/ |
||||
float32_t fir_coeffs[53] = { |
||||
-0.000206064,-0.000525129,-0.00083518, -0.000774011, 2.5925E-05, |
||||
0.001614912, 0.003431897, 0.004335125, 0.003127158, -0.000566047, |
||||
-0.005566484,-0.009192163,-0.008417443,-0.001801824, 0.008839149, |
||||
0.018273049, 0.019879265, 0.009349346,-0.011696836, -0.034389317, |
||||
-0.045008839,-0.030706279, 0.013824834, 0.082060266, 0.156328996, |
||||
0.213799940, 0.235420817, 0.213799940, 0.156328996, 0.082060266, |
||||
0.013824834,-0.030706279,-0.045008839,-0.034389317, -0.011696836, |
||||
0.009349346, 0.019879265, 0.018273049, 0.008839149, -0.001801824, |
||||
-0.008417443,-0.009192163,-0.005566484,-0.000566047, 0.003127158, |
||||
0.004335125, 0.003431897, 0.001614912, 2.5925E-05, -0.000774011, |
||||
-0.000835180,-0.000525129,-0.000206064 }; |
||||
|
||||
// 8-pole Biquad fc=0.0025fs, -80 dB Iowa Hills
|
||||
// This is roughly the narrowest that doesn't have
|
||||
// artifacts from numerical errors more than about
|
||||
// 0.001 radians (0.06 deg), per experiments using F32.
|
||||
// b0,b1,b2,a1,a2 for each BiQuad. Start with stage 0
|
||||
float32_t iir_coeffs[5 * N_STAGES]={ |
||||
0.08686551007982608, |
||||
-0.1737214710369926, |
||||
0.08686551007982608, |
||||
1.9951804375779567, |
||||
-0.9951899867006161, |
||||
// and stage 1
|
||||
0.20909791845765324, |
||||
-0.4181667739705088, |
||||
0.20909791845765324, |
||||
1.9965910753714984, |
||||
-0.9966201383162961, |
||||
// stage 2
|
||||
0.18360046797931723, |
||||
-0.3671514768697197, |
||||
0.18360046797931723, |
||||
1.9981966389027592, |
||||
-0.998246097991674, |
||||
// stage 3
|
||||
0.03079484444321144, |
||||
-0.061529427044071175, |
||||
0.03079484444321144, |
||||
1.999421284937329, |
||||
-0.9994815467796806}; |
||||
|
||||
// ARM DSP Math library IIR filter instance
|
||||
arm_biquad_casd_df1_inst_f32 iir_inst; |
||||
|
||||
// And a FIR type, as either can be used via begin()
|
||||
arm_fir_instance_f32 fir_inst; |
||||
|
||||
// Delay line space for the FIR
|
||||
float32_t FIRStateF32[AUDIO_BLOCK_SAMPLES + NFIR_MAX]; |
||||
|
||||
// Delay line space for the Biquad, each arranged as {x[n-1], x[n-2], y[n-1], y[n-2]}
|
||||
float32_t IIRStateF32[4 * N_STAGES]; |
||||
}; |
||||
#endif |
@ -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; kk<nFIR; kk++) // To be sure, zero the coefficients
|
||||
cf32f[kk] = 0.0f; |
||||
|
||||
// Convert dB to Voltage ratios, frequencies to fractions of sampling freq
|
||||
if(nBands <2 || nBands>50) return ERR_EQ_BANDS; |
||||
for (i=0; i<nBands; i++) { |
||||
aVolts[i]=powf(10.0, (0.05*adb[i])); |
||||
fNorm[i]=feq[i]/sample_rate_Hz; |
||||
} |
||||
|
||||
/* Find FIR coefficients, the Fourier transform of the frequency
|
||||
* response. This is done by dividing the response into a sequence |
||||
* of nBands rectangular frequency blocks, each of a different level. |
||||
* We can precalculate the Fourier transform for each rectangular band. |
||||
* The linearity of the Fourier transform allows us to sum the transforms |
||||
* of the individual blocks to get pre-windowed coefficients. As follows |
||||
* |
||||
* Numbering example for nFIR==199: |
||||
* Subscript 0 to 98 is 99 taps; 100 to 198 is 99 taps; 99+1+99=199 taps |
||||
* The center coef ( for nFIR=199 taps, nHalfFIR=99 ) is a
|
||||
* special case that comes from sin(0)/0 and treated first: |
||||
*/ |
||||
cf32f[nHalfFIR] = 2.0f*(aVolts[0]*fNorm[0]); // Coefficient "99"
|
||||
for(i=1; i<nBands; i++) { |
||||
cf32f[nHalfFIR] += 2.0f*aVolts[i]*(fNorm[i]-fNorm[i-1]); |
||||
} |
||||
for (j=1; j<=nHalfFIR; j++) { // Coefficients "100 to 198"
|
||||
q = MF_PI*(float32_t)j; |
||||
// First, deal with the zero frequency end band that is "low-pass."
|
||||
cf32f[j+nHalfFIR] = aVolts[0]*sinf(fNorm[0]*2.0*q)/q; |
||||
// and then the rest of the bands that have low and high frequencies
|
||||
for(i=1; i<nBands; i++) |
||||
cf32f[j+nHalfFIR] += aVolts[i]*( (sinf(fNorm[i]*2.0*q)/q) - (sinf(fNorm[i-1]*2.0*q)/q) ); |
||||
} |
||||
|
||||
/* At this point, the cf32f[] coefficients are simply truncated sin(x)/x shapes, creating
|
||||
* very high sidelobe responses. To reduce the sidelobes, a windowing function is applied. |
||||
* This has the side affect of increasing the rate of cutoff for sharp frequency changes. |
||||
* The only windowing function available here is that of James Kaiser. This has a number |
||||
* of desirable features. The tradeoff of sidelobe level versus cutoff rate is variable.
|
||||
* We specify it in terms of kdb, the highest sidelobe, in dB, next to a sharp cutoff. For |
||||
* calculating the windowing vector, we need a parameter beta, found as follows: |
||||
*/ |
||||
if (kdb<0) return ERR_EQ_SIDELOBES; |
||||
if (kdb>50) |
||||
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; i<nFreq; i++) { |
||||
bt = cf32f[nHalfFIR];//bt = 0.5f*cf32f[nHalfFIR]; // Center coefficient
|
||||
for (j=0; j<nHalfFIR; j++) // Add in the others twice, as they are symmetric
|
||||
bt += 2.0f*cf32f[j]*cosf(piOnNfreq*(float32_t)((nHalfFIR-j)*i)); |
||||
rdb[i] = 20.0f*log10f(fabsf(bt)); // Convert to dB
|
||||
} |
||||
} |
@ -0,0 +1,178 @@ |
||||
/*
|
||||
* AudioFilterEqualizer_F32 |
||||
* |
||||
* Created: Bob Larkin W7PUA 8 May 2020 |
||||
*
|
||||
* This is a direct translation of the receiver audio equalizer written |
||||
* by this author for the open-source DSP-10 radio in 1999. See |
||||
* http://www.janbob.com/electron/dsp10/dsp10.htm and
|
||||
* http://www.janbob.com/electron/dsp10/uhf3_35a.zip
|
||||
*
|
||||
* Credit and thanks to PJRC, Paul Stoffregen and Teensy products for the audio |
||||
* system and library that this is built upon as well as the float32 |
||||
* work of Chip Audette embodied in the OpenAudio_ArduinoLibrary. Many thanks |
||||
* for the library structures and wonderful Teensy products. |
||||
* |
||||
* This equalizer is specified by an array of 'nBands' frequency bands |
||||
* each of of arbitrary frequency span. The first band always starts at |
||||
* 0.0 Hz, and that value is not entered. Each band is specified by the upper |
||||
* frequency limit to the band. |
||||
* The last band always ends at half of the sample frequency, which for 44117 Hz |
||||
* sample frequency would be 22058.5. Each band is specified by its upper |
||||
* frequency in an .INO supplied array feq[]. The dB level of that band is |
||||
* specified by a value, in dB, arranged in an .INO supplied array |
||||
* aeq[]. Thus a trivial bass/treble control might look like: |
||||
* nBands = 3; |
||||
* feq[] = {300.0, 1500.0, 22058.5}; |
||||
* float32_t bass = -2.5; // in dB, relative to anything
|
||||
* float32_t treble = 6.0; |
||||
* aeq[] = {bass, 0.0, treble}; |
||||
* |
||||
* It may be obvious that this equalizer is a more general case of the common |
||||
* functions such as low-pass, band-pass, notch, etc. For instance, a pair |
||||
* of band pass filters would look like: |
||||
* nBands = 5; |
||||
* feq[] = {500.0, 700.0, 2000.0, 2200.0, 22058.5}; |
||||
* aeq[] = {-100.0, 0.0, -100.0, 2.0, -100.0}; |
||||
* where we added 2 dB of gain to the 2200 to 2400 Hz filter, relative to the 500 |
||||
* to 700 Hz band. |
||||
* |
||||
* An octave band equalizer is made by starting at some low frequency, say 40 Hz for the |
||||
* first band. The lowest frequency band will be from 0.0 Hz up to that first frequency. |
||||
* Next multiply the first frequency by 2, creating in our example, a band from 40.0 |
||||
* to 80 Hz. This is continued until the last frequency is about 22058 Hz. |
||||
* This works out to require 10 bands, as follows: |
||||
* nBands = 10; |
||||
* feq[] = { 40.0, 80.0, 160.0, 320.0, 640.0, 1280.0, 2560.0, 5120.0, 10240.0, 22058.5}; |
||||
* aeq[] = { 5.0, 4.0, 2.0, -3.0, -4.0, -1.0, 3.0, 6.0, 3.0, 0.5 }; |
||||
* |
||||
* For a "half octave" equalizer, multiply each upper band limit by the square root of 2 = 1.414 |
||||
* to get the next band limit. For that case, feq[] would start with a sequence |
||||
* like 40, 56.56, 80.00, 113.1, 160.0, ... for a total of about 20 bands. |
||||
* |
||||
* How well all of this is achieved depends on the number of FIR coefficients |
||||
* being used. In the Teensy 3.6 / 4.0 the resourses allow a hefty number, |
||||
* say 201, of coefficients to be used without stealing all the processor time |
||||
* (see Timing, below). The coefficient and FIR memory is sized for a maximum of |
||||
* 250 coefficients, but can be recompiled for bigger with the define FIR_MAX_COEFFS. |
||||
* To simplify calculations, the number of FIR coefficients should be odd. If not |
||||
* odd, the number will be reduced by one, quietly. |
||||
* |
||||
* If you try to make the bands too narrow for the number of FIR coeffficients, |
||||
* the approximation to the desired curve becomes poor. This can all be evaluated |
||||
* by the function getResponse(nPoints, pResponse) which fills an .INO-supplied array |
||||
* pResponse[nPoints] with the frequency response of the equalizer in dB. The nPoints |
||||
* are spread evenly between 0.0 and half of the sample frequency. |
||||
* |
||||
* Initialization is a 2-step process. This makes it practical to change equalizer |
||||
* levels on-the-fly. The constructor starts up with a 4-tap FIR setup for direct |
||||
* pass through. Then the setup() in the .INO can specify the equalizer. |
||||
* The newEqualizer() function has several parameters, the number of equalizer bands, |
||||
* the frequencies of the bands, and the sidelobe level. All of these can be changed |
||||
* dynamically. This function can be changed dynamically, but it may be desireable to |
||||
* mute the audio during the change to prevent clicks. |
||||
*
|
||||
* This 16-bit integer version adjusts the maximum coefficient size to scale16 in the calls |
||||
* to both equalizerNew() and getResponse(). Broadband equalizers can work with full-scale |
||||
* 32767.0f sorts of levels, where narrow band filtering may need smaller values to |
||||
* prevent overload. Experiment and check carefully. Use lower values if there are doubts. |
||||
*
|
||||
* For a pass-through function, something like this (which can be intermixed with fancy equalizers): |
||||
* float32_t fBand[] = {10000.0f, 22058.5f};
|
||||
* float32_t dbBand[] = {0.0f, 0.0f}; |
||||
* equalize1.equalizerNew(2, &fBand[0], &dbBand[0], 4, &equalizeCoeffs[0], 30.0f, 32767.0f); |
||||
*
|
||||
* Measured timing of update() for a 128 sample block, Teensy 3.6: |
||||
* Fixed time 13 microseconds |
||||
* Per FIR Coefficient time 2.5 microseconds |
||||
* Total for 199 FIR Coefficients = 505 microseconds (17.4% of 44117 Hz available time) |
||||
*
|
||||
* Per FIR Coefficient, Teensy 4.0, 0.44 microseconds |
||||
*
|
||||
* Copyright (c) 2020 Bob Larkin |
||||
* Any snippets of code from PJRC or Chip Audette used here brings with it |
||||
* the associated license. |
||||
*
|
||||
* In addition, work here is covered by MIT LIcense: |
||||
* |
||||
* 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 _filter_equalizer_f32_h |
||||
#define _filter_equalizer_f32_h |
||||
|
||||
#include "Arduino.h" |
||||
#include "AudioStream_F32.h" |
||||
#include "arm_math.h" |
||||
#include "mathDSP_F32.h" |
||||
|
||||
#ifndef MF_PI |
||||
#define MF_PI 3.1415926f |
||||
#endif |
||||
|
||||
// Temporary timing test
|
||||
#define TEST_TIME_EQ 0 |
||||
|
||||
#define EQUALIZER_MAX_COEFFS 251 |
||||
|
||||
#define ERR_EQ_BANDS 1 |
||||
#define ERR_EQ_SIDELOBES 2 |
||||
#define ERR_EQ_NFIR 3 |
||||
|
||||
class AudioFilterEqualizer_F32 : public AudioStream_F32 |
||||
{ |
||||
//GUI: inputs:1, outputs:1 //this line used for automatic generation of GUI node
|
||||
//GUI: shortName:filter_Equalizer
|
||||
public: |
||||
AudioFilterEqualizer_F32(void): AudioStream_F32(1,inputQueueArray) { |
||||
// Initialize FIR instance (ARM DSP Math Library) with default simple passthrough FIR
|
||||
arm_fir_init_f32(&fir_inst, nFIR, (float32_t *)cf32f, &StateF32[0], (uint32_t)block_size); |
||||
} |
||||
AudioFilterEqualizer_F32(const AudioSettings_F32 &settings): AudioStream_F32(1,inputQueueArray) { |
||||
block_size = settings.audio_block_samples; |
||||
sample_rate_Hz = settings.sample_rate_Hz; |
||||
arm_fir_init_f32(&fir_inst, nFIR, (float32_t *)cf32f, &StateF32[0], (uint32_t)block_size); |
||||
} |
||||
|
||||
uint16_t equalizerNew(uint16_t _nBands, float32_t *feq, float32_t *adb, |
||||
uint16_t _nFIR, float32_t *_cf32f, float32_t kdb); |
||||
void getResponse(uint16_t nFreq, float32_t *rdb); |
||||
void update(void); |
||||
|
||||
private: |
||||
audio_block_f32_t *inputQueueArray[1]; |
||||
uint16_t block_size = AUDIO_BLOCK_SAMPLES; |
||||
float32_t firStart[4] = {0.0, 1.0, 0.0, 0.0}; // Initialize to passthrough
|
||||
float32_t* cf32f = firStart; // pointer to current coefficients
|
||||
uint16_t nFIR = 4; // Number of coefficients
|
||||
uint16_t nBands = 2; |
||||
float32_t sample_rate_Hz = AUDIO_SAMPLE_RATE; |
||||
|
||||
// *Temporary* - TEST_TIME allows measuring time in microseconds for each part of the update()
|
||||
#if TEST_TIME_EQ |
||||
elapsedMicros tElapse; |
||||
int32_t iitt = 999000; // count up to a million during startup
|
||||
#endif |
||||
// ARM DSP Math library filter instance
|
||||
arm_fir_instance_f32 fir_inst; |
||||
float32_t StateF32[AUDIO_BLOCK_SAMPLES + EQUALIZER_MAX_COEFFS]; // max, max
|
||||
}; |
||||
#endif |
||||
|
||||
|
@ -0,0 +1,259 @@ |
||||
/* AudioFilterFIRGeneralr_F32.cpp
|
||||
* |
||||
* Bob Larkin, W7PUA 20 May 2020 (c) |
||||
* |
||||
* 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. |
||||
*/ |
||||
/* A math note - The design of the FIR filter from the frequency response requires
|
||||
* taking a discrete Fourier transform. The calculation of the achieved |
||||
* frequency response requires another Fourier transform. Good news is |
||||
* that this is arithmetically simple, because of |
||||
* working with real (not complex) inputs and requiring the FIR be linear |
||||
* phase. One of a number of references is C.S. Burris and T. W. Parks, |
||||
* "Digital Filter Design." That book also has excellent sections on |
||||
* Chebychev error (Parks-McClellan-Rabiner) FIR filters and IIR Filters. |
||||
*/ |
||||
|
||||
#include "AudioFilterFIRGeneral_F32.h" |
||||
|
||||
void AudioFilterFIRGeneral_F32::update(void) { |
||||
audio_block_f32_t *blockIn, *blockOut; |
||||
|
||||
#if TEST_TIME_FIRG |
||||
if (iitt++ >1000000) 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<nFIR; i++) // To be sure, zero the coefficients
|
||||
cf32f[i] = 0.0f; |
||||
for(i=0; i<(nFIR+AUDIO_BLOCK_SAMPLES); i++) // And the storage
|
||||
pStateArray[i] = 0.0f; |
||||
|
||||
// Convert dB to "Voltage" ratios, frequencies to fractions of sampling freq
|
||||
// Borrow pStateArray, as it has not yet gone to ARM math
|
||||
for(i=0; i<=nHalfFIR; i++) |
||||
pStateArray[i] = powf(10.0, 0.05f*adb[i]); |
||||
|
||||
/* Find FIR coefficients, the Fourier transform of the frequency
|
||||
* response. This general frequency description has half as many |
||||
* frequency dB points as FIR coefficients. If nFIR is odd, |
||||
* the number of frequency points is nHalfFIR = (nFIR - 1)/2. |
||||
* |
||||
* Numbering example for nFIR==21: |
||||
* Odd nFIR: Subscript 0 to 9 is 10 taps; 11 to 20 is 10 taps; 10+1+10=21 taps |
||||
* |
||||
* Numbering example for nFIR==20:
|
||||
* Even nFIR: Subscript 0 to 9 is 10 taps; 10 to 19 is 10 taps; 10+10=20 taps |
||||
*/ |
||||
float32_t oneOnNFIR = 1.0f / (float32_t)nFIR; |
||||
if(even) { |
||||
for(n=0; n<nHalfFIR; n++) { // Over half of even nFIR FIR coeffs
|
||||
cf32f[n] = pStateArray[0]; // DC term
|
||||
for(k=1; k<nHalfFIR; k++) { // Over all frequencies, half of nFIR
|
||||
num = (M - (float32_t)n) * (float32_t)k;
|
||||
cf32f[n] += 2.0f*pStateArray[k]*cosf(MF_TWOPI*num*oneOnNFIR); |
||||
} |
||||
cf32f[n] *= oneOnNFIR; // Divide by nFIR
|
||||
}
|
||||
} |
||||
else { // nFIR is odd
|
||||
for(n=0; n<=nHalfFIR; n++) { // Over half of FIR coeffs, nFIR odd
|
||||
cf32f[n] = pStateArray[0]; |
||||
for(k=1; k<=nHalfFIR; k++) { |
||||
num = (float32_t)((nHalfFIR - n)*k); |
||||
cf32f[n] += 2.0f*pStateArray[k]*cosf(MF_TWOPI*num*oneOnNFIR); |
||||
} |
||||
cf32f[n] *= oneOnNFIR; |
||||
}
|
||||
} |
||||
/* At this point, the cf32f[] coefficients are simply truncated, creating
|
||||
* high sidelobe responses. To reduce the sidelobes, a windowing function is applied. |
||||
* This has the side affect of increasing the rate of cutoff for sharp frequency transition. |
||||
* The only windowing function available here is that of James Kaiser. This has a number |
||||
* of desirable features. The sidelobes drop off as the frequency away from a transition. |
||||
* Also, the tradeoff of sidelobe level versus cutoff rate is variable. |
||||
* Here we specify it in terms of kdb, the highest sidelobe, in dB, next to a sharp cutoff. For |
||||
* calculating the windowing vector, we need a parameter beta, found as follows: |
||||
*/ |
||||
if (kdb<0.0f) return ERR_FIRGEN_SIDELOBES; |
||||
if (kdb < 20.0f) |
||||
beta = 0.0; |
||||
else |
||||
beta = -2.17+0.17153*kdb-0.0002841*kdb*kdb; // Within a dB or so
|
||||
// Note: i0f is the fp zero'th order modified Bessel function (see mathDSP_F32.h)
|
||||
kbes = 1.0f / mathEqualizer.i0f(beta); // An additional derived parameter used in loop
|
||||
scaleXn2 = 4.0f / ( ((float32_t)nFIR - 1.0f)*((float32_t)nFIR - 1.0f) ); // Needed for even & odd
|
||||
if (even) { |
||||
// nHalfFIR = nFIR/2; for nFIR even
|
||||
for (n=0; n<nHalfFIR; n++) { // For 20 Taps, this is 0 to 9
|
||||
xn2 = 0.5f+(float32_t)n; |
||||
xn2 = scaleXn2*xn2*xn2; |
||||
WindowWt=kbes*(mathEqualizer.i0f(beta*sqrt(1.0-xn2))); |
||||
if(kdb > 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; i<nFreq; i++) { |
||||
bt = 0.0; |
||||
for (n=0; n<nHalfFIR; n++) // Add in terms twice, as they are symmetric
|
||||
bt += 2.0f*cf32f[n]*cosf(piOnNfreq*((M-(float32_t)n)*(float32_t)i)); |
||||
rdb[i] = 20.0f*log10f(fabsf(bt)); // Convert to dB
|
||||
} |
||||
} |
||||
else { // it is odd
|
||||
nHalfFIR = (nFIR - 1)/2; |
||||
for (i=0; i<nFreq; i++) { |
||||
bt = cf32f[nHalfFIR]; // Center coefficient
|
||||
for (n=0; n<nHalfFIR; n++) // Add in the others twice, as they are symmetric
|
||||
bt += 2.0f*cf32f[n]*cosf(piOnNfreq*(float32_t)((nHalfFIR-n)*i)); |
||||
rdb[i] = 20.0f*log10f(fabsf(bt)); |
||||
} |
||||
}
|
||||
} |
@ -0,0 +1,188 @@ |
||||
/*
|
||||
* AudioFilterFIRGeneral_F32 |
||||
* |
||||
* Created: Bob Larkin W7PUA 20 May 2020 |
||||
* |
||||
* Credit and thanks to PJRC, Paul Stoffregen and Teensy products for the audio |
||||
* system and library that this is built upon as well as the float32 |
||||
* work of Chip Audette embodied in the OpenAudio_ArduinoLibrary. Many thanks |
||||
* for the library structures and wonderful Teensy products. |
||||
* |
||||
* There are enough different FIR filter Audio blocks to need a summary. Here goes: |
||||
*
|
||||
* AudioFilterFIR (Teensy Audio Library by PJRC) handles 16-bit integer data, |
||||
* and a maximum of 200 FIR coefficients, even only. (taps). For Teensy Audio. |
||||
* AudioFilterFIR_F32 (OpenAudio_ArduinoLibrary by Chip Audette) handles 32-bit floating point |
||||
* data and a maximum of 200 taps. Can be used with a mix of Teensy Audio |
||||
* and 32-bit OpenAudio_Arduino blocks. |
||||
* AudioFilterFIRGeneral_F32 (This block) handles 32-bit floating point |
||||
* data and very large number of taps, even or odd. Can be used with a mix of Teensy Audio |
||||
* and 32-bit OpenAudio_Arduino blocks. This includes the design of an |
||||
* arbitrary frequency response using Kaiser windows. |
||||
* AudioFilterFIRGeneral_I16 Same as this block, but data is 16-bit integer and |
||||
* the number of taps must be even. |
||||
* AudioFilterEqualizer_F32 handles 32-bit floating point data and 250 maximum taps |
||||
* even or odd. Can be used with a mix of Teensy Audio |
||||
* and 32-bit OpenAudio_Arduino blocks. This includes the design of an |
||||
* arbitrary frequency response using multiple flat response steps. A Kaiser windows |
||||
* is used. |
||||
* AudioFilterEqualizer_I16 handles 16-bit integer data and 250 maximum taps, |
||||
* even only. |
||||
* |
||||
* FIR filters suffer from needing considerable computation of the multiply-and-add |
||||
* sort. This limits the number of taps that can be used, but less so as time goes by. |
||||
* In particular, the Teensy 4.x, if it *did nothing but* FIR calculations, could |
||||
* use about 6000 taps inmonaural, which is a huge number. But, this also
|
||||
* suggests that if the filtering task is an important function of a project, |
||||
* using, say 2000 taps is practical. |
||||
*
|
||||
* FIR filters can be (and are here) implemented to have symmetrical coefficients. This |
||||
* results in constant delay at all frequencies (linear phase). For some applications this can |
||||
* be an important feature. Sometimes it is suggested that the FIR should not be |
||||
* used because of the latency it creates. Note that if constant delay is needed, the FIR |
||||
* implementation does this with minimum latency. |
||||
*
|
||||
* For this block, AudioFilterFIRGeneral_F32, memory storage for the FIR |
||||
* coefficiients as well as working storage for the ARM FIR routine is provided |
||||
* by the calling .INO. This allows large FIR sizes without always tying up a |
||||
* big memory block. |
||||
*
|
||||
* This block normally calculates the FIR coefficients using a Fourier transform |
||||
* of the desired amplitude response and a Kaiser window. This flexability requires |
||||
* the calling .INO to provide an array of response levels, in relative dB, |
||||
* that is half the length of the number of FIR taps. An example of this entry is a |
||||
* 300 point low-pass filter with a cutoff at 10% of the 44.1 kHz sample frequency: |
||||
* for(int i=0; i<150; i++) { |
||||
* if (i<30) dbA[i] = 0.0f; |
||||
* else dbA[i] = -140.0f; |
||||
* } |
||||
* firg1.FIRGeneralNew(&dbA[0], 300, &equalizeCoeffs[0], 50.0f, &workSpace[0]); |
||||
*
|
||||
* As an alternate to inputting the response function, the FIR coefficients can be |
||||
* entered directly using LoadCoeffs(nFIR, cf32f, *pStateArray). This is a very quick |
||||
* operation as only pointers to coefficients are involved. Several filters can be |
||||
* stored in arrays and switched quickly this way. If this is done, pStateArray[] |
||||
* as initially setup should be large enough for all filters. There will be "clicks" |
||||
* associated with filter changes and these may need to be muted.
|
||||
* |
||||
* How well the desired response is achieved depends on the number of FIR coefficients |
||||
* being used. As noted above, for some applications it may be desired to use |
||||
* large numbers of taps. The achieved response can be evaluated
|
||||
* by the function getResponse(nPoints, pResponse) which fills an .INO-supplied array |
||||
* pResponse[nPoints] with the frequency response of the equalizer in dB. The nPoints |
||||
* are spread evenly between 0.0 and half of the sample frequency. |
||||
* |
||||
* Initialization is a 2-step process. This makes it practical to change equalizer |
||||
* levels on-the-fly. The constructor starts up with a 4-tap FIR setup for direct |
||||
* pass through. Then the setup() in the .INO can specify the equalizer. |
||||
* The newEqualizer() function has several parameters, the number of equalizer bands, |
||||
* the frequencies of the bands, and the sidelobe level. All of these can be changed |
||||
* dynamically. This function can be changed dynamically, but it may be desireable to |
||||
* mute the audio during the change to prevent clicks. |
||||
* |
||||
* Measured timing of update() for L or R 128 sample block, Teensy 3.6: |
||||
* Fixed time 13 microseconds |
||||
* Per FIR Coefficient time 2.5 microseconds (about 51E6 multiply and accumulates/sec) |
||||
* Total for 199 FIR Coefficients = 505 microseconds (17.4% of 44117 Hz available time) |
||||
* Total for stereo is twice those numbers. |
||||
* Measured timing of update() for L or R 128 sample block, Teensy 4.0: |
||||
* Fixed time 1.4 microseconds |
||||
* Per FIR Coefficient time 0.44 microseconds (about 290E6 multiply and accumulate/sec) |
||||
* Total for 199 FIR Coefficients = 90 microseconds (3.1% of 44117 Hz available time) |
||||
* Total for stereo is twice those numbers |
||||
* Measured for FIRGeneralNew(), T4.0, to design a 4000 tap FIR is 14 sec. This goes |
||||
* with the square of the number of taps. |
||||
* Measured for getResponse() for nFIR=4000 and nFreq=5000, T4.0, is about a minute. |
||||
* |
||||
* Functions for the AudioFilterFIRGeneral_F32 object are |
||||
* FIRGeneralNew(*adb, nFIR, cf32f, kdb, *pStateArray); // to design and use an adb[]
|
||||
* frequency response. |
||||
* LoadCoeffs(nFIR, cf32f, *pStateArray); // To directly load FIR coefficients cf32f[]
|
||||
* getResponse(nFreq, *rdb); // To obtain the amplitude response in dB, rdb[]
|
||||
*
|
||||
* Status: Tested T3.6, T4.0 No known bugs |
||||
*
|
||||
* Examples: TestFIRGeneralLarge4.ino TestFIRGeneralLarge5.ino |
||||
*
|
||||
* Copyright (c) 2020 Bob Larkin |
||||
* Any snippets of code from PJRC or Chip Audette used here brings with it |
||||
* the associated license. |
||||
* |
||||
* In addition, work here is covered by MIT LIcense: |
||||
* |
||||
* 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 _filter_FIR_general_f32_h |
||||
#define _filter_FIR_general_f32_h |
||||
|
||||
#include "Arduino.h" |
||||
#include "AudioStream_F32.h" |
||||
#include "arm_math.h" |
||||
#include "mathDSP_F32.h" |
||||
|
||||
#ifndef MF_PI |
||||
#define MF_PI 3.1415926f |
||||
#endif |
||||
|
||||
// Temporary timing test
|
||||
#define TEST_TIME_FIRG 0 |
||||
|
||||
#define ERR_FIRGEN_BANDS 1 |
||||
#define ERR_FIRGEN_SIDELOBES 2 |
||||
#define ERR_FIRGEN_NFIR 3 |
||||
|
||||
class AudioFilterFIRGeneral_F32 : public AudioStream_F32 |
||||
{ |
||||
//GUI: inputs:1, outputs:1 //this line used for automatic generation of GUI node
|
||||
//GUI: shortName:filter_Equalizer
|
||||
public: |
||||
AudioFilterFIRGeneral_F32(void): AudioStream_F32(1,inputQueueArray) { |
||||
// Initialize FIR instance (ARM DSP Math Library) with default simple passthrough FIR
|
||||
arm_fir_init_f32(&fir_inst, nFIR, (float32_t *)cf32f, &StateF32[0], (uint32_t)block_size); |
||||
} |
||||
AudioFilterFIRGeneral_F32(const AudioSettings_F32 &settings): AudioStream_F32(1,inputQueueArray) { |
||||
block_size = settings.audio_block_samples; |
||||
sample_rate_Hz = settings.sample_rate_Hz; |
||||
arm_fir_init_f32(&fir_inst, nFIR, (float32_t *)cf32f, &StateF32[0], (uint32_t)block_size); |
||||
} |
||||
|
||||
uint16_t FIRGeneralNew(float32_t *adb, uint16_t _nFIR, float32_t *_cf32f, float32_t kdb, float32_t *pStateArray); |
||||
uint16_t LoadCoeffs(uint16_t _nFIR, float32_t *_cf32f, float32_t *pStateArray); |
||||
void getResponse(uint16_t nFreq, float32_t *rdb); |
||||
void update(void |