From f917ff4f53d51ca26bd7cf2eefc705bc57e21002 Mon Sep 17 00:00:00 2001 From: boblark Date: Fri, 11 Feb 2022 10:48:48 -0800 Subject: [PATCH] Draft for adding ConvFilt's --- AudioFilterConvolution_F32.cpp | 293 +++++++++--------- AudioFilterConvolution_F32.h | 179 ++++++----- .../TestConvolutinalFilter.ino | 74 +++-- 3 files changed, 298 insertions(+), 248 deletions(-) diff --git a/AudioFilterConvolution_F32.cpp b/AudioFilterConvolution_F32.cpp index caf5049..f66cd0f 100644 --- a/AudioFilterConvolution_F32.cpp +++ b/AudioFilterConvolution_F32.cpp @@ -2,8 +2,8 @@ ****************************************************************************** * @file AudioFilterConvolution_F32.cpp * @author Giuseppe Callipo - IK8YFW - ik8yfw@libero.it - * @version V1.0.0 - * @date 02-05-2021 + * @version V2.0.0 + * @date 06-02-2021 * @brief F32 Filter Convolution * ****************************************************************************** @@ -17,144 +17,130 @@ 1) Class refactoring, change some methods visibility; 2) Filter coefficients calculation included into class; 3) Change the class for running in both with F32 - OpenAudio_ArduinoLibrary for Teensy; + OpenAudio_ArduinoLibrary for Teensy; 4) Added initFilter method for single anf fast initialization and on the fly reinititializzation; - 5) Optimize it to use as output audio filter on SDR receiver. + 5) Optimize it to use as output audio filter on SDR receiver. + 6) Optimize the time execution *******************************************************************/ + // Revised for OpenAudio_Arduino Teensy F32 library, 8 Feb 2022 #include "AudioFilterConvolution_F32.h" -boolean AudioFilterConvolution_F32::begin(int status) -{ - enabled = status; - return(true); -} - void AudioFilterConvolution_F32::passThrough(int stat) { passThru=stat; } -void AudioFilterConvolution_F32::impulse(float32_t *FIR_coef) { -// arm_q15_to_float(coefs, FIR_coef, 513); // convert int_buffer to float 32bit - int k = 0; - int i = 0; - enabled = 0; // shut off audio stream while impulse is loading - for (i = 0; i < (FFT_length / 2) + 1; i++) - { - FIR_filter_mask[k++] = FIR_coef[i]; - FIR_filter_mask[k++] = 0; - } +// Function to pre-calculate the multiplying frequency function, the "mask." +void AudioFilterConvolution_F32::impulse(float32_t *FIR_coef) +{ + uint32_t k = 0; + uint32_t i = 0; + enabled = 0; // shut off audio stream while impulse is loading + for (i = 0; i < (FFT_length / 2) + 1; i++) + { + FIR_filter_mask[k++] = FIR_coef[i]; + FIR_filter_mask[k++] = 0; + } - for (i = FFT_length + 1; i < FFT_length * 2; i++) - { - FIR_filter_mask[i] = 0.0; - } - arm_cfft_f32( &arm_cfft_sR_f32_len1024, FIR_filter_mask, 0, 1); - for (int i = 0; i < 1024; i++) { - // Serial.println(FIR_filter_mask[i] * 32768); - } - // for 1st time thru, zero out the last sample buffer to 0 - memset(last_sample_buffer_L, 0, sizeof(last_sample_buffer_L)); - state = 0; - enabled = 1; //enable audio stream again + for (i = FFT_length + 1; i < FFT_length * 2; i++) + { + FIR_filter_mask[i] = 0.0; + } + arm_cfft_f32( &arm_cfft_sR_f32_len1024, FIR_filter_mask, 0, 1); + + // for 1st time thru, zero out the last sample buffer to 0 + arm_fill_f32(0, last_sample_buffer_L, BUFFER_SIZE *4); + + state = 0; + enabled = 1; //enable audio stream again } void AudioFilterConvolution_F32::update(void) { - audio_block_f32_t *block; - float32_t *bp; + audio_block_f32_t *block; + float32_t *bp; - if (enabled != 1 ) return; - block = receiveWritable_f32(0); // MUST be Writable, as convolution results are written into block - if (block) { - switch (state) { - case 0: - bp = block->data; - for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { - buffer[i] = *bp++; - } - bp = block->data; - for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { - *bp++ = tbuffer[i]; // tbuffer contains results of last FFT/multiply/iFFT processing (convolution filtering) - } - AudioStream_F32::transmit(block); - AudioStream_F32::release(block); - state = 1; - break; - case 1: - bp = block->data; - for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { - buffer[128+i] = *bp++; - } - bp = block->data; - for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { - *bp++ = tbuffer[i+128]; // tbuffer contains results of last FFT/multiply/iFFT processing (convolution filtering) - } - AudioStream_F32::transmit(block); - AudioStream_F32::release(block); - state = 2; - break; - case 2: - bp = block->data; - for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { - buffer[256 + i] = *bp++; - } - bp = block->data; - for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { - *bp++ = tbuffer[i+256]; // tbuffer contains results of last FFT/multiply/iFFT processing (convolution filtering) - } - AudioStream_F32::transmit(block); - AudioStream_F32::release(block); - state = 3; - break; - case 3: - bp = block->data; - for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { - buffer[384 + i] = *bp++; - } - bp = block->data; - for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { - *bp++ = tbuffer[i + 384]; // tbuffer contains results of last FFT/multiply/iFFT processing (convolution filtering) - } - AudioStream_F32::transmit(block); - AudioStream_F32::release(block); - state = 0; - // 4 blocks are in- now do the FFT1024,complex multiply and iFFT1024 on 512samples of data - // using the overlap/add method - // 1st convert Q15 samples to float - //arm_q15_to_float(buffer, float_buffer_L, 512); - arm_copy_f32 (buffer, float_buffer_L, 512); - // float_buffer samples are now standardized from > -1.0 to < 1.0 - if (passThru ==0) { - memset(FFT_buffer + 1024, 0, sizeof(FFT_buffer) / 2); // zero pad last half of array- necessary to prevent aliasing in FFT - //fill FFT_buffer with current audio samples - k = 0; - for (i = 0; i < 512; i++) - { - FFT_buffer[k++] = float_buffer_L[i]; // real - FFT_buffer[k++] = float_buffer_L[i]; // imag - } - // calculations are performed in-place in FFT routines - arm_cfft_f32(&arm_cfft_sR_f32_len1024, FFT_buffer, 0, 1);// perform complex FFT - arm_cmplx_mult_cmplx_f32(FFT_buffer, FIR_filter_mask, iFFT_buffer, FFT_length); // complex multiplication in Freq domain = convolution in time domain - arm_cfft_f32(&arm_cfft_sR_f32_len1024, iFFT_buffer, 1, 1); // perform complex inverse FFT - k = 0; - l = 1024; - for (int i = 0; i < 512; i++) { - float_buffer_L[i] = last_sample_buffer_L[i] + iFFT_buffer[k++]; // this performs the "ADD" in overlap/Add - last_sample_buffer_L[i] = iFFT_buffer[l++]; // this saves 512 samples (overlap) for next time around - k++; - l++; + if (enabled != 1 ) return; + block = receiveWritable_f32(0); // MUST be Writable, as convolution results are written into block + if (block) { + switch (state) { + case 0: + if (passThru ==0) { + arm_cmplx_mult_cmplx_f32(FFT_buffer, FIR_filter_mask, iFFT_buffer, FFT_length); // complex multiplication in Freq domain = convolution in time domain + arm_cfft_f32(&arm_cfft_sR_f32_len1024, iFFT_buffer, 1, 1); // perform complex inverse FFT + k = 0; + l = 1024; + for (int i = 0; i < 512; i++) { + buffer[i] = last_sample_buffer_L[i] + iFFT_buffer[k++]; // this performs the "ADD" in overlap/Add + last_sample_buffer_L[i] = iFFT_buffer[l++]; // this saves 512 samples (overlap) for next time around + k++; + l++; + } + // convert floats to Q15 and save in temporary array tbuffer + arm_copy_f32 (&buffer[0], &tbuffer[0], BUFFER_SIZE*4); + } + bp = block->data; + for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { + buffer[i] = *bp; + *bp++ = tbuffer[i]; + } + AudioStream_F32::transmit(block); + AudioStream_F32::release(block); + state = 1; + break; + + case 1: + bp = block->data; + for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { + buffer[128+i] = *bp; + *bp++ = tbuffer[i+128]; + } + AudioStream_F32::transmit(block); + AudioStream_F32::release(block); + state = 2; + break; + + case 2: + bp = block->data; + for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { + buffer[256 + i] = *bp; + *bp++ = tbuffer[i+256]; // tbuffer contains results of last FFT/multiply/iFFT processing (convolution filtering) + } + + AudioStream_F32::transmit(block); + AudioStream_F32::release(block); + // zero pad last half of array- necessary to prevent aliasing in FFT + arm_fill_f32(0, FFT_buffer + 1024, FFT_length); + state = 3; + break; + + case 3: + bp = block->data; + for (int i = 0; i < AUDIO_BLOCK_SAMPLES; i++) { + buffer[384 + i] = *bp; + *bp++ = tbuffer[i + 384]; // tbuffer contains results of last FFT/multiply/iFFT processing (convolution filtering) + } + AudioStream_F32::transmit(block); + AudioStream_F32::release(block); + state = 0; + // 4 blocks are in- now do the FFT1024,complex multiply and iFFT1024 on 512samples of data + // using the overlap/add method + if (passThru ==0) { + //fill FFT_buffer with current audio samples + k = 0; + for (i = 0; i < 512; i++) + { + FFT_buffer[k++] = buffer[i]; // real + FFT_buffer[k++] = buffer[i]; // imag + } + // calculations are performed in-place in FFT routines + arm_cfft_f32(&arm_cfft_sR_f32_len1024, FFT_buffer, 0, 1);// perform complex FFT + } //end if passTHru + break; } - } //end if passTHru - // convert floats to Q15 and save in temporary array tbuffer - //arm_float_to_q15(&float_buffer_L[0], &tbuffer[0], BUFFER_SIZE*4); - arm_copy_f32 (&float_buffer_L[0], &tbuffer[0], BUFFER_SIZE*4); - break; } - } } float32_t AudioFilterConvolution_F32::Izero (float32_t x) @@ -177,7 +163,7 @@ float32_t AudioFilterConvolution_F32::Izero (float32_t x) } // END Izero float AudioFilterConvolution_F32::m_sinc(int m, float fc) -{ // fc is f_cut/(Fsamp/2) +{ // fc is f_cut/(Fsamp/2) // m is between -M and M step 2 // float x = m*PIH; @@ -187,20 +173,25 @@ float AudioFilterConvolution_F32::m_sinc(int m, float fc) return sinf(x*fc)/(fc*x); } -void AudioFilterConvolution_F32::calc_FIR_coeffs (float * coeffs, int numCoeffs, float32_t fc, float32_t Astop, int type, float dfc, float Fsamprate) - // pointer to coefficients variable, no. of coefficients to calculate, frequency where it happens, stopband attenuation in dB, - // filter type, half-filter bandwidth (only for bandpass and notch) - { // modified by WMXZ and DD4WH after - // Wheatley, M. (2011): CuteSDR Technical Manual. www.metronix.com, pages 118 - 120, FIR with Kaiser-Bessel Window - // assess required number of coefficients by +void AudioFilterConvolution_F32::calc_FIR_coeffs + (float *coeffs, int numCoeffs, float32_t fc, float32_t Astop, + int type, float dfc, float Fsamprate){ + // pointer to coefficients variable, no. of coefficients to calculate, + // frequency where it happens, stopband attenuation in dB, + // filter type, half-filter bandwidth (only for bandpass and notch). + // Modified by WMXZ and DD4WH after + // Wheatley, M. (2011): CuteSDR Technical Manual. www.metronix.com, + // pages 118 - 120, FIR with Kaiser-Bessel Window. + // Assess required number of coefficients by // numCoeffs = (Astop - 8.0) / (2.285 * TPI * normFtrans); - // selecting high-pass, numCoeffs is forced to an even number for better frequency response + // selecting high-pass, numCoeffs is forced to an even number for + // better frequency response int ii,jj; float32_t Beta; float32_t izb; float fcf = fc; int nc = numCoeffs; - fc = fc / Fsamprate; + fc = 2.0f * fc / Fsamprate; // Corrected dfc = dfc / Fsamprate; // calculate Kaiser-Bessel window shape factor beta from stop-band attenuation if (Astop < 20.96) @@ -211,38 +202,40 @@ void AudioFilterConvolution_F32::calc_FIR_coeffs (float * coeffs, int numCoeffs, Beta = 0.5842 * powf((Astop - 20.96), 0.4) + 0.07886 * (Astop - 20.96); izb = Izero (Beta); - if(type == 0) // low pass filter -// { fcf = fc; - { fcf = fc * 2.0; + if(type == LOWPASS) + { fcf = fc; nc = numCoeffs; } - else if(type == 1) // high-pass filter + else if(type == HIGHPASS) { fcf = -fc; nc = 2*(numCoeffs/2); } - else if ((type == 2) || (type==3)) // band-pass filter + else if ((type == BANDPASS) || (type==BANDREJECT)) { fcf = dfc; nc = 2*(numCoeffs/2); // maybe not needed } - else if (type==4) // Hilbert transform - { - nc = 2*(numCoeffs/2); + +/* + else if (type==HILBERT) + { + nc = 2*(numCoeffs/2); // clear coefficients for(ii=0; ii< 2*(nc-1); ii++) coeffs[ii]=0; // set real delay - coeffs[nc]=1; - + coeffs[nc]=1; // <<<<2000.0f ? 100 : 5)); +// toneDet1.frequency(f, (f>2000.0f ? 100 : 5)); toneDet2.frequency(f, (f>2000.0f ? 100 : 5)); delay(50); - while(!toneDet1.available())delay(2) ; - toneDet1.read(); +// while(!toneDet1.available())delay(2) ; +// toneDet1.read(); while(!toneDet2.available())delay(2) ; toneDet2.read(); while(!toneDet2.available())delay(2) ; Serial.print(f); Serial.print(", "); while(!toneDet2.available())delay(2) ; - Serial.print(20.0f*log10f(toneDet2.read() ) ); - Serial.print(", "); - while(!toneDet1.available())delay(2) ; - Serial.println(20.0f*log10f(toneDet1.read() ) ); + Serial.println(20.0f*log10f(toneDet2.read() ) ); } }