Added pureSpectrum filtering option.

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);
   }

synth_sin_cos_f32.cpp

@@ -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,7 +83,7 @@ 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;
@@ -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::transmit(blockC, 1);
     AudioStream_F32::release (blockC);
     return;
 }
-

synth_sin_cos_f32.h

@@ -3,6 +3,8 @@
  *
  * 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
  *
@@ -89,6 +91,37 @@ public:
         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
@@ -146,6 +179,8 @@ public:
       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