diff --git a/examples/TestFFT4096iq/TestFFT4096iq.ino b/examples/TestFFT4096iq/TestFFT4096iq.ino index c11d615..768fe98 100644 --- a/examples/TestFFT4096iq/TestFFT4096iq.ino +++ b/examples/TestFFT4096iq/TestFFT4096iq.ino @@ -6,11 +6,15 @@ // Serial Print out powers of all 4096 bins in // dB relative to Sine Wave Full Scale +// Rev- added i to print; added filters on oscillators. 23 Feb 2022 + // Public Domain #include "OpenAudio_ArduinoLibrary.h" #include "AudioStream_F32.h" +int jj = 0; + // GUItool: begin automatically generated code AudioSynthSineCosine_F32 sine_cos1; //xy=76,532 AudioAnalyzeFFT4096_IQ_F32 FFT4096iq1; //xy=243,532 @@ -37,10 +41,13 @@ void setup(void) { // or pick any old frequency sine_cos1.frequency(1000.0f); - // elect the output format + // Engage the identical BP Filters on sine/cosine outputs (true). + sine_cos1.pureSpectrum(true); + + // Select the output format FFT_POWER, FFT_RMS or FFT_DBFS FFT4096iq1.setOutputType(FFT_DBFS); - // Select the wndow function + // Select the wndow function from these four: //FFT4096iq1.windowFunction(AudioWindowNone); //FFT4096iq1.windowFunction(AudioWindowHanning4096); //FFT4096iq1.windowFunction(AudioWindowKaiser4096, 55.0f); @@ -50,24 +57,37 @@ void setup(void) { // float* pw = FFT4096iq1.getWindow(); // Print window // for (int i=0; i<4096; i++) Serial.println(pw[i], 7); - // xAxis, bit 0 left/right; bit 1 low to high; default 0X03 + // x-Axis direction and offset for sine to I and cosine to Q. + // If xAxis=0 f=fs/2 in middle, f=0 on right edge + // If xAxis=1 f=fs/2 in middle, f=0 on left edge + // If xAxis=2 f=fs/2 on left edge, f=0 in middle + // If xAxis=3 f=fs/2 on right edge, f=0 in middle + // xAxis=3 is mathematically proper -fs/2 start to +fs/2 stop. FFT4096iq1.setXAxis(0X03); - FFT4096iq1.setNAverage(1); - delay(100); + FFT4096iq1.setNAverage(5); } void loop(void) { static bool doPrint=true; float *pPwr; // Print output, once - if( FFT4096iq1.available() && doPrint ) { + if( FFT4096iq1.available() && jj++>2 && doPrint ) + { pPwr = FFT4096iq1.getData(); for(int i=0; i<4096; i++) + { + // Serial.print((int)((float32_t)i * 44100.0/4096.0)); // Print freq + Serial.print(i); // FFT Output index (0, 4095) + Serial.print(" "); Serial.println(*(pPwr + i), 8 ); + } doPrint = false; } - Serial.print(" Audio MEM Float32 Peak: "); - Serial.println(AudioMemoryUsageMax_F32()); - delay(500); + if(doPrint) + { + Serial.print(" Audio MEM Float32 Peak: "); + Serial.println(AudioMemoryUsageMax_F32()); + } + delay(100); } diff --git a/synth_sin_cos_f32.cpp b/synth_sin_cos_f32.cpp index c55bd8a..7d01aca 100644 --- a/synth_sin_cos_f32.cpp +++ b/synth_sin_cos_f32.cpp @@ -1,5 +1,5 @@ /* synth_sin_cos_f32.cpp - * + * * SynthSinCos_F32 Bob Larkin April 17, 2020 * * Based on Chip Audette's OpenAudio sine(), that was @@ -44,13 +44,13 @@ void AudioSynthSineCosine_F32::update(void) { 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++) { @@ -62,15 +62,16 @@ void AudioSynthSineCosine_F32::update(void) { /* 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 */ - + // Corrected + // blockS->data[i] = a + 0.001953125*(b-a)*deltaPhase; /* Linear interpolation process */ + blockS->data[i] = a+(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 */ + blockC->data[i] = a +(b-a)*deltaPhase; /* Linear interpolation process */ } } else { // Do a more flexible update, i.e., not doSimple @@ -82,8 +83,8 @@ void AudioSynthSineCosine_F32::update(void) { /* 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 */ - + blockS->data[i] = amplitude_pk*(a +(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; @@ -92,13 +93,21 @@ void AudioSynthSineCosine_F32::update(void) { /* 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 */ + blockC->data[i] = amplitude_pk*(a +(b-a)*deltaPhase); /* Linear interpolation process */ } } + + // For higher frequencies, an optional bandpass filter the output + // This does a pass through for lower frequencies + if(doPureSpectrum) + { + arm_biquad_cascade_df1_f32(&bq_instS, blockS->data, blockS->data, 128); + arm_biquad_cascade_df1_f32(&bq_instC, blockC->data, blockC->data, 128); + } + AudioStream_F32::transmit(blockS, 0); - AudioStream_F32::release (blockS); + 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 index 19ba124..1d492b3 100644 --- a/synth_sin_cos_f32.h +++ b/synth_sin_cos_f32.h @@ -1,8 +1,10 @@ /* synth_sin_cos_f32.h * AudioSynthSineCosine_F32 - * + * * Status: Checked for function and accuracy. 19 April 2020 * 10 March 2021 Corrected Interpolation equations. Bob L + * 23 Feb 2022 Added "pure spectrum" filtering. RSL + * This adds the function pureSpectrum(bool _setPure); // Default: false * * Created: Bob Larkin 15 April 2020 * @@ -19,19 +21,19 @@ * (-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. - * + * 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 @@ -79,16 +81,47 @@ class AudioSynthSineCosine_F32 : public AudioStream_F32 { 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) { + 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; + + // Find coeff for 2 stages of BPF to remove harmoncs + // Always compute these in case pureSpectrum is enabled later. + if(freq > 0.003f*sample_rate_Hz) + { + float32_t q = 20.0f; + float32_t w0 = freq * (2.0f * 3.141592654f / sample_rate_Hz); + float32_t alpha = sin(w0) / (q * 2.0); + float32_t scale = 1.0f / (1.0f + alpha); + /* b0 */ coeff32[0] = alpha * scale; + /* b1 */ coeff32[1] = 0; + /* b2 */ coeff32[2] = (-alpha) * scale; + /* a1 */ coeff32[3] = -(-2.0 * cos(w0)) * scale; + /* a2 */ coeff32[4] = -(1.0 - alpha) * scale; + /* b0 */ coeff32[5] = coeff32[0]; + /* b1 */ coeff32[6] = coeff32[1]; + /* b2 */ coeff32[7] = coeff32[2]; + /* a1 */ coeff32[8] = coeff32[3]; + /* a2 */ coeff32[9] = coeff32[4]; + arm_biquad_cascade_df1_init_f32( &bq_instS, 2, coeff32, state32S ); + arm_biquad_cascade_df1_init_f32( &bq_instC, 2, coeff32, state32C ); + } + else + { + for(int ii=0; ii<10; ii++) // Coeff for BiQuad BPF + coeff32[ii] = 0.0; + coeff32[0] = 1.0; // b0 = 1 for pass through + coeff32[5] = 1.0; + arm_biquad_cascade_df1_init_f32( &bq_instS, 2, coeff32, state32S ); + arm_biquad_cascade_df1_init_f32( &bq_instC, 2, coeff32, state32C ); + } } /* Externally, phase comes in the range (0,2*M_PI) keeping with C math functions @@ -102,7 +135,7 @@ public: 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. @@ -123,7 +156,7 @@ public: return; } - // Speed up calculations by setting phaseS_C=90deg, amplitude=1 + // 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; @@ -132,7 +165,7 @@ public: amplitude_pk = 1.0f; } return; - } + } void setSampleRate_Hz(float32_t fs_Hz) { // Check freq range @@ -145,7 +178,9 @@ public: if(bl > 128) bl = 128; block_length = bl; } - + + void pureSpectrum(bool _setPure) { doPureSpectrum = _setPure; } + virtual void update(void); private: @@ -159,5 +194,11 @@ private: // 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; + bool doPureSpectrum = false; // Adds bandpass filter (not normally needed) + float32_t coeff32[10]; // 2 biquad stages for filtering output shared S&C + float32_t state32S[8]; + arm_biquad_casd_df1_inst_f32 bq_instS; // ARM DSP Math library filter instance. + float32_t state32C[8]; + arm_biquad_casd_df1_inst_f32 bq_instC; // ARM DSP Math library filter instance. }; #endif