From bc7ad85336d033232235727f06568ee6fcef4ff9 Mon Sep 17 00:00:00 2001
From: boblark <bob@janbob.com>
Date: Thu, 18 Mar 2021 13:55:02 -0700
Subject: [PATCH] Update real input 1024 fft

---
 analyze_fft1024_F32.cpp              | 132 +++++++++-------
 analyze_fft1024_F32.h                | 100 +++++++-----
 examples/TestFFT1024/TestFFT1024.ino | 222 +++++++++++++++++++++++++--
 3 files changed, 344 insertions(+), 110 deletions(-)

diff --git a/analyze_fft1024_F32.cpp b/analyze_fft1024_F32.cpp
index d7b9076..0cdb094 100644
--- a/analyze_fft1024_F32.cpp
+++ b/analyze_fft1024_F32.cpp
@@ -1,5 +1,6 @@
 /* analyze_fft1024_F32.cpp      Converted from Teensy I16 Audio Library
  * This version uses float F32 inputs.  See comments at analyze_fft1024_F32.h
+ * Converted to use half-length FFT 17 March 2021  RSL
  *
  * Conversion parts copyright (c) Bob Larkin 2021
  *
@@ -33,36 +34,32 @@
 #include "analyze_fft1024_F32.h"
 
 // Move audio data from an audio_block_f32_t to the FFT instance buffer.
-// This is for 128 numbers per block
-static void copy_to_fft_buffer(void *destination, const void *source)
-{
+// This is for 128 numbers per block, only.
+static void copy_to_fft_buffer(void *destination, const void *source)  {
     const float *src = (const float *)source;
     float *dst = (float *)destination;
 
     for (int i=0; i < 128; i++) {
-        *dst++ = *src++;     // real sample
-        *dst++ = 0.0f;       // 0 for Imag
+        *dst++ = *src++;     // real sample for half-length FFT
+        }
     }
-}
 
-static void apply_window_to_fft_buffer(void *buffer, const void *window)
-{
-    float *buf = (float *)buffer;      // 0th entry is real (do window) 1th is imag
+static void apply_window_to_fft_buffer(void *buffer, const void *window) {
+    float *buf = (float *)buffer; 
     const float *win = (float *)window;
 
-    for (int i=0; i < 1024; i++)
-        buf[2*i] *= *win++; 
-}
+    for(int i=0; i<NFFT; i++)
+       buf[i] *= *win++;
+    }
 
 void AudioAnalyzeFFT1024_F32::update(void)  {
     audio_block_f32_t *block;
-    uint32_t tt;
+    float outDC=0.0f;
+    float magsq=0.0f;
 
     block = AudioStream_F32::receiveReadOnly_f32();
     if (!block) return;
 
-    // tt=micros();
-
     switch (state) {
     case 0:
         blocklist[0] = block;
@@ -82,6 +79,55 @@ void AudioAnalyzeFFT1024_F32::update(void)  {
         break;
     case 4:
         blocklist[4] = block;
+        // Now the post FT processing for using half-length transform
+        // FFT was in state==7, but it loops around to here.  Does some
+        // zero data at startup that is harmless.
+        count++;      // Next do non-coherent averaging
+        for(int i=0; i<NFFT_D2; i++)  {
+            if(i>0) {
+               float rns = 0.5f*(fft_buffer[2*i] +   fft_buffer[NFFT-2*i]);
+               float ins = 0.5f*(fft_buffer[2*i+1] + fft_buffer[NFFT-2*i+1]);
+               float rnd = 0.5f*(fft_buffer[2*i] -   fft_buffer[NFFT-2*i]); 
+               float ind = 0.5f*(fft_buffer[2*i+1] - fft_buffer[NFFT-2*i+1]);
+
+               float xr = rns + cosN[i]*ins - sinN[i]*rnd;
+               float xi = ind - sinN[i]*ins - cosN[i]*rnd;
+               magsq = xr*xr + xi*xi;
+               }
+            else  {
+               magsq = outDC*outDC;    // Do the DC term
+               }
+
+            if(count==1)               // Starting new average
+               sumsq[i] = magsq;
+            else if (count<=nAverage)  // Adding on to average
+               sumsq[i] += magsq;
+            }
+
+        if (count >= nAverage) {       // Average is finished
+        // Set outputflag false here to minimize reads of output[]  data
+        // when it is being updated. 
+        outputflag = false;
+            count = 0;
+            float inAf = 1.0f/(float)nAverage;
+            float kMaxDB = 20.0*log10f((float)NFFT_D2);  // 54.1854 for 1024
+
+            for(int i=0; i<NFFT_D2; i++)  {
+               if(outputType==FFT_RMS)
+                  output[i] = sqrtf(inAf*sumsq[i]);
+               else if(outputType==FFT_POWER)
+                  output[i] = inAf*sumsq[i];
+               else if(outputType==FFT_DBFS)  {
+                  if(sumsq[i]>0.0f)
+                     output[i] = 10.0f*log10f(inAf*sumsq[i]) - kMaxDB;  // Scaled to FS sine wave
+                  else
+                     output[i] = -193.0f;   // lsb for 23 bit mantissa
+                  }
+               else
+                  output[i] = 0.0f;
+               }    // End, set output[i] over all NFFT_D2
+            outputflag = true;
+            }    // End of average is finished
         state = 5;
         break;
     case 5:
@@ -96,55 +142,32 @@ void AudioAnalyzeFFT1024_F32::update(void)  {
         blocklist[7] = block;
         // We have 4 previous blocks pointed to by blocklist[]:
         copy_to_fft_buffer(fft_buffer+0x000, blocklist[0]->data);
-        copy_to_fft_buffer(fft_buffer+0x100, blocklist[1]->data);
-        copy_to_fft_buffer(fft_buffer+0x200, blocklist[2]->data);
-        copy_to_fft_buffer(fft_buffer+0x300, blocklist[3]->data);
+        copy_to_fft_buffer(fft_buffer+0x080, blocklist[1]->data);
+        copy_to_fft_buffer(fft_buffer+0x100, blocklist[2]->data);
+        copy_to_fft_buffer(fft_buffer+0x180, blocklist[3]->data);
         // and 4 new blocks, just gathered:
-        copy_to_fft_buffer(fft_buffer+0x400, blocklist[4]->data);
-        copy_to_fft_buffer(fft_buffer+0x500, blocklist[5]->data);
-        copy_to_fft_buffer(fft_buffer+0x600, blocklist[6]->data);
-        copy_to_fft_buffer(fft_buffer+0x700, blocklist[7]->data);
- 
+        copy_to_fft_buffer(fft_buffer+0x200, blocklist[4]->data);
+        copy_to_fft_buffer(fft_buffer+0x280, blocklist[5]->data);
+        copy_to_fft_buffer(fft_buffer+0x300, blocklist[6]->data);
+        copy_to_fft_buffer(fft_buffer+0x380, blocklist[7]->data);
+
         if (pWin)
            apply_window_to_fft_buffer(fft_buffer, window);
 
+        outDC = 0.0f;
+        for(int i=0; i<NFFT; i++)
+            outDC += fft_buffer[i];
+        outDC /= ((float)NFFT);
+
 #if defined(__IMXRT1062__)
         // Teensyduino core for T4.x supports arm_cfft_f32
         // arm_cfft_f32 (const arm_cfft_instance_f32 *S, float32_t *p1, uint8_t ifftFlag, uint8_t bitReverseFlag)
         arm_cfft_f32 (&Sfft, fft_buffer, 0, 1);
 #else
-        // For T3.x go back to old (deprecated) style
-        arm_cfft_radix4_f32(&fft_inst, fft_buffer);
+        // For T3.x go back to old (deprecated) style (check radix2/radix4)<<<
+        arm_cfft_radix2_f32(&fft_inst, fft_buffer);
 #endif
-
-        count++;      // Next do non-coherent averaging
-        for (int i=0; i < 512; i++) {
-            float magsq = fft_buffer[2*i]*fft_buffer[2*i] + fft_buffer[2*i+1]*fft_buffer[2*i+1];
-            if(count==1)               // Starting new average
-               sumsq[i] = magsq;
-            else if (count<=nAverage)  // Adding on to average
-               sumsq[i] += magsq;
-            }
-        if (count >= nAverage) {       // Average is finished
-            count = 0;
-            float inAf = 1.0f/(float)nAverage;
-            for(int i=0; i<512; i++)  {
-               if(outputType==FFT_RMS)
-                  output[i] = sqrtf(inAf*sumsq[i]);
-               else if(outputType==FFT_POWER)
-                  output[i] = inAf*sumsq[i];
-               else if(outputType==FFT_DBFS)  {
-                  if(sumsq[i]>0.0f)
-                     output[i] = 10.0f*log10f(inAf*sumsq[i])-54.1854f;  // Scaled to FS sine wave
-                  else
-                     output[i] = -193.0f;   // lsb for 23 bit mantissa
-                  }
-               else
-                  output[i] = 0.0f;
-               }    // End, set output[i] over all 512
-
-            outputflag = true;
-            }    // End of average is finished
+        // FFT output is now in fft_buffer.  Pick up processing at state==4.
 
         AudioStream_F32::release(blocklist[0]);
         AudioStream_F32::release(blocklist[1]);
@@ -158,5 +181,4 @@ void AudioAnalyzeFFT1024_F32::update(void)  {
         state = 4;
         break;
      }     // End switch(state)
-     // Serial.print("uSec: "); Serial.println(micros()-tt);
   }        // End update()
diff --git a/analyze_fft1024_F32.h b/analyze_fft1024_F32.h
index 728e80b..b0d5cc3 100644
--- a/analyze_fft1024_F32.h
+++ b/analyze_fft1024_F32.h
@@ -26,7 +26,7 @@
  * THE SOFTWARE.
  */
 
-/* Moved directly I16 to F32. Bob Larkin 16 Feb 2021
+/* Translated from  I16 to F32. Bob Larkin 16 Feb 2021
  * Does real input FFT of 1024 points. Output is not audio, and is magnitude
  * only.  Multiple output formats of RMS (same as I16 version, and default),
  * Power or dBFS (full scale).  Output can be bin by bin or a pointer to
@@ -45,13 +45,13 @@
  *   void setOutputType(int _type)
  *   void setNAverage(int nAverage)
  *
- * Timing, max is longest update() time:
- *   T3.6 Windowed, Power Out, 975 uSec
- *   T3.6 Windowed, RMS out,  1016 uSec
- *   T3.6 Windowed, dBFS out, 1591 uSec
+ * Timing, max is longest update() time. Comparison is using full complex FFT
+ * and no load sharing on "states".
+ *   T3.6 Windowed, Power Out, 682 uSec (was 975 w/ 1024 FFT)
+ *   T3.6 Windowed, dBFS out, 834 uSec (was 1591 w/1024 FFT)
  *   No Window saves 60 uSec on T3.6 for any output.
- *   T4.0 Windowed, Ave=1, Power Out, 156 uSec
- *   T4.0 Windowed, Ave=1, dBFS Out,  302 uSec
+ *   T4.0 Windowed, Power Out, 54 uSec (was 156 w/1024  FFT)
+ *   T4.0 Windowed, dBFS Out, 203 uSec (was 302 w/1024 FFT)
  * Scaling:
  *   Full scale for floating point DSP is a nebulous concept.  Normally the
  *   full scale is -1.0 to +1.0.  This is an unscaled FFT and for a sine
@@ -64,6 +64,9 @@
  *   when building the INO.
  */
 // Fixed float/int problem in read(first, last).  RSL 3 Mar 21
+// Converted to using half-size FFT for real input, with no zero inputs.
+// See E. Oran Brigham and many other FFT references.  16 March 2021 RSL
+// Moved post-FFT calculations to state 4 to load share.  RSL 18 Mar 2021
 
 #ifndef analyze_fft1024_F32_h_
 #define analyze_fft1024_F32_h_
@@ -76,6 +79,15 @@
 #include "arm_const_structs.h"
 #endif
 
+// Doing an FFT with NFFT real inputs
+#define NFFT 1024
+#define NFFT_M1 NFFT-1
+#define NFFT_D2  NFFT/2
+#define NFFT_D2M1 (NFFT/2)-1
+#define NFFT_X2  NFFT*2
+#define FFT_PI  3.14159265359f
+#define FFT_2PI 6.28318530718f
+
 #define FFT_RMS 0
 #define FFT_POWER 1
 #define FFT_DBFS 2
@@ -96,13 +108,17 @@ public:
 #if defined(__IMXRT1062__)
         // Teensy4 core library has the right files for new FFT
         // arm CMSIS library has predefined structures of type arm_cfft_instance_f32
-        Sfft = arm_cfft_sR_f32_len1024;   // This is one of the structures
+        Sfft = arm_cfft_sR_f32_len512;   // Like this.  Changes with size <<<
 #else
-        arm_cfft_radix4_init_f32(&fft_inst, 1024, 0, 1); // for T3.x
+        arm_cfft_radix2_init_f32(&fft_inst, NFFT_D2, 0, 1);   // for T3.x (check radix2/radix4)<<<
 #endif
-        // This is always 128 block size.  Any sample rate.  No use of "setings"
-        // Bob:  Revisit this to use CMSIS fast fft in place of cfft.  faster?
+        // This class is always 128 block size.  Any sample rate.  No use of "settings"
         useHanningWindow();
+        // Factors for using half size complex FFT
+        for(int n=0; n<NFFT_D2; n++)  {
+           sinN[n] = sinf(FFT_PI*((float)n)/((float)NFFT_D2));
+           cosN[n] = cosf(FFT_PI*((float)n)/((float)NFFT_D2));
+           }
     }
 
     // Inform that the output is available for read()
@@ -116,7 +132,7 @@ public:
 
     // Output data from a single bin is transferred
     float read(unsigned int binNumber) {
-        if (binNumber>511 || binNumber<0) return 0.0;
+        if (binNumber>NFFT_D2M1 || binNumber<0) return 0.0;
         return output[binNumber];
         }
 
@@ -128,8 +144,8 @@ public:
             binLast = binFirst;
             binFirst = tmp;
         }
-        if (binFirst > 511) return 0.0;
-        if (binLast > 511) binLast = 511;
+        if (binFirst > NFFT_D2M1) return 0.0;
+        if (binLast > NFFT_D2M1) binLast = NFFT_D2M1;
         float sum = 0.0f;
         do {
             sum += output[binFirst++];
@@ -138,7 +154,7 @@ public:
     }
 
     int windowFunction(int wNum) {
-      if(wNum == AudioWindowKaiser1024)
+      if(wNum == AudioWindowKaiser1024)  // Changes with size <<<
          return -1;                 // Kaiser needs the kdb
       windowFunction(wNum, 0.0f);
       return 0;
@@ -149,14 +165,14 @@ public:
       pWin = window;
       if(wNum == NO_WINDOW)
          pWin = NULL;
-      else if (wNum == AudioWindowKaiser1024)  {
+      else if (wNum == AudioWindowKaiser1024)  {  // Changes with size <<<
          if(_kdb<20.0f)
             kd = 20.0f;
          else
             kd = _kdb;
          useKaiserWindow(kd);
          }
-      else if (wNum == AudioWindowBlackmanHarris1024)
+      else if (wNum == AudioWindowBlackmanHarris1024)  // Changes with size <<<
          useBHWindow();
      else
          useHanningWindow();   // Default
@@ -177,7 +193,7 @@ public:
     // Bring custom window from the INO
     void putWindow(float *pwin)  {
        float *p = window;
-       for(int i=0; i<1024; i++)
+       for(int i=0; i<NFFT; i++)
           *p++ = *pwin++;
        }
 
@@ -194,21 +210,27 @@ public:
     virtual void update(void);
 
 private:
-    float output[512];
-    float sumsq[512];
-    float window[1024];
+    float output[NFFT_D2];
+    float sumsq[NFFT_D2];  // Accumulates averages of outputs
+    float window[NFFT];
     float *pWin = window;
+    float fft_buffer[NFFT];
+
+    // The cosN and sinN would seem to be twidddle factors.  Someday
+    // look at this and see if they can be stolen from arm math/DSP.
+    float cosN[NFFT_D2];
+    float sinN[NFFT_D2];
+
     audio_block_f32_t *blocklist[8];
-    float fft_buffer[2048];
     uint8_t state = 0;
     bool outputflag = false;
     audio_block_f32_t *inputQueueArray[1];
 #if defined(__IMXRT1062__)
-  // For T4.x
-  // const static arm_cfft_instance_f32   arm_cfft_sR_f32_len1024;
+  // For T4.x, 512 length for real 1024 input, etc.
+  // const static arm_cfft_instance_f32   arm_cfft_sR_f32_len512;
   arm_cfft_instance_f32 Sfft;
 #else
-  arm_cfft_radix4_instance_f32 fft_inst;
+  arm_cfft_radix2_instance_f32 fft_inst;  // Check radix2/radix4 <<<
 #endif
     int outputType = FFT_RMS;  //Same type as I16 version init
     int nAverage = 1;
@@ -216,22 +238,22 @@ private:
 
     // The Hann window is a good all-around window
     void useHanningWindow(void) {
-        for (int i=0; i < 1024; i++) {
-           // 2*PI/1023 = 0.006141921
-           window[i] = 0.5*(1.0 - cosf(0.006141921f*(float)i));
+        for (int i=0; i < NFFT; i++) {
+           float kx = FFT_2PI/((float)NFFT_M1); // 0.006141921 for 1024
+           window[i] = 0.5f*(1.0f - cosf(kx*(float)i));
         }
     }
 
     // Blackman-Harris produces a first sidelobe more than 90 dB down.
     // The price is a bandwidth of about 2 bins.  Very useful at times.
     void useBHWindow(void) {
-        for (int i=0; i < 1024; i++) {
-           float kx = 0.006141921;  // 2*PI/1023
+        for (int i=0; i < NFFT; i++) {
+           float kx = FFT_2PI/((float)NFFT_M1);  // 0.006141921 for 1024
            int ix = (float) i;
-           window[i] = 0.35875 -
-                       0.48829*cosf(     kx*ix) +
-                       0.14128*cosf(2.0f*kx*ix) -
-                       0.01168*cosf(3.0f*kx*ix);
+           window[i] = 0.35875f -
+                       0.48829f*cosf(     kx*ix) +
+                       0.14128f*cosf(2.0f*kx*ix) -
+                       0.01168f*cosf(3.0f*kx*ix);
         }
     }
 
@@ -252,12 +274,12 @@ private:
 
        // 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
-       for (int n=0; n<512; n++) {
+       for (int n=0; n<NFFT_D2; n++) {
           xn2 = 0.5f+(float32_t)n;
-          // 4/(1023^2)=0.00000382215877f
-          xn2 = 0.00000382215877f*xn2*xn2;
-          window[511 - n]=kbes*(mathEqualizer.i0f(beta*sqrtf(1.0-xn2)));
-          window[512 + n] = window[511 - n];
+          float kx = 4.0f/(((float)NFFT_M1) * ((float)NFFT_M1));  // 0.00000382215877f for 1024
+          xn2 = kx*xn2*xn2;
+          window[NFFT_D2M1 - n]=kbes*(mathEqualizer.i0f(beta*sqrtf(1.0-xn2)));
+          window[NFFT_D2 + n] = window[NFFT_D2M1 - n];
           }
        }
     };
diff --git a/examples/TestFFT1024/TestFFT1024.ino b/examples/TestFFT1024/TestFFT1024.ino
index 52e52be..7789232 100644
--- a/examples/TestFFT1024/TestFFT1024.ino
+++ b/examples/TestFFT1024/TestFFT1024.ino
@@ -1,3 +1,193 @@
+
+#if 0
+
+#include "Arduino.h"
+#include "AudioStream_F32.h"
+#include "arm_math.h"
+#include "mathDSP_F32.h"
+#if defined(__IMXRT1062__)
+#include "arm_const_structs.h"
+#endif
+
+#define NFFT 1024
+#define NFFT_D2 NFFT/2
+
+#define FFT_PI 3.14159265359f
+#define PI2 2.0f*FFT_PI
+
+void setup(void) {
+    float x[NFFT];       // Real DFT input
+    float Xre[NFFT], Xim[NFFT]; // DFT of x
+    float P[NFFT];           // power spectrum of x
+    float kf, nf;
+    float fft_buffer[2*NFFT];  // 2 is fo CMSIS FFT
+    float sinN[NFFT_D2];
+    float cosN[NFFT_D2];
+
+    uint32_t tt;
+    // Instantiate FFT, T4.x ONLY
+    arm_cfft_instance_f32   Sfft;
+    Sfft = arm_cfft_sR_f32_len1024;
+
+   // Instantiate FFT, T4.x ONLY
+    arm_cfft_instance_f32   Sfft_128;
+    Sfft_128 = arm_cfft_sR_f32_len512;
+
+    Serial.begin(300); // Any speed works
+    delay(1000);
+
+    // Factors for using half size complex FFT
+    for(int n=0; n<NFFT_D2; n++)  {
+       sinN[n] = sinf(FFT_PI*((float)n)/((float)NFFT_D2));
+       cosN[n] = cosf(FFT_PI*((float)n)/((float)NFFT_D2));
+       }
+
+    // The signal, a sine wave that fits integer number (2 here) of times
+    Serial.println("The input waveform, NFFT points");
+    for (int n=0; n<NFFT; n++) {
+        x[n] = sinf(5.0f*FFT_PI*((float)n)/((float)NFFT_D2));
+        Serial.println(x[n], 8);
+        }
+    Serial.println();
+
+    Serial.println("The DFT by full NxN DFT, k, Real, Imag, Power");
+    // Calculate DFT of x[n] with NFFT point DFT
+    for (int k=0 ; k<NFFT ; ++k) {
+        kf = (float)k;
+        // Real part of X[k]
+        Xre[k] = 0.0f;
+        Xim[k] = 0.0f;
+        for (int n=0 ; n<NFFT ; ++n)  {
+           nf = (float)n;
+           Xre[k] += x[n] * cos(nf * kf * PI2 / ((float)NFFT));
+           Xim[k] -= x[n] * sin(nf * kf * PI2 / ((float)NFFT));
+           }
+        // Power at kth frequency bin
+        P[k] = 10.0f*log10f(Xre[k]*Xre[k] + Xim[k]*Xim[k]);
+        Serial.print(k); Serial.print(",");
+        Serial.print(Xre[k],8); Serial.print(",");
+        Serial.print(Xim[k],8); Serial.print(",");
+        Serial.println(P[k],3);
+        }
+
+
+    Serial.println("And now the same thing with FFT size NFFT/2");
+
+    for (int k = 0; k < NFFT_D2; k++) {  // For each output element
+        kf = (float)k;
+        float sumreal = 0;
+        float sumimag = 0;
+        for (int n = 0; n < NFFT_D2; n++) {  // For each input element
+            nf = (float)n;
+            float angle = PI2 * nf * kf / ((float)NFFT_D2);
+            sumreal +=  x[2*n] * cos(angle) + x[2*n+1] * sin(angle);
+            sumimag += -x[2*n] * sin(angle) + x[2*n+1] * cos(angle);
+        }
+        Xre[k] = sumreal;
+        Xim[k] = sumimag;
+    }
+
+    // Rearrange the outputs to look like the FFT
+    for(int k=0; k<NFFT_D2; k++) {
+       fft_buffer[2*k] = Xre[k];
+       fft_buffer[2*k+1] = Xim[k];
+       }
+
+
+    // Now the post FT processing for using half-length transform
+    Xre[0] = 0.0f;
+    for(int n=0; n<NFFT; n++)
+       Xre[0] += x[n]/((float)NFFT);    // DC real
+    Xim[0] = 0.0f;                      // DC Imag
+    P[0] = 10.0f*log10f(Xre[0]*Xre[0]);
+    // And the non-DC values
+    for(int i=1; i<NFFT_D2; i++)  {
+        float rns = 0.5f*(fft_buffer[2*i] + fft_buffer[NFFT-2*i]);
+        float ins = 0.5f*(fft_buffer[2*i+1] + fft_buffer[NFFT-2*i+1]);
+        float rnd = 0.5f*(fft_buffer[2*i] - fft_buffer[NFFT-2*i]);
+        float ind = 0.5f*(fft_buffer[2*i+1] - fft_buffer[NFFT-2*i+1]);
+        Xre[i] = rns + cosN[i]*ins - sinN[i]*rnd;
+        Xim[i] = ind - sinN[i]*ins - cosN[i]*rnd;
+        P[i] = 10.0f*log10f(Xre[i]*Xre[i] + Xim[i]*Xim[i]);
+        }
+
+    for(int k=0; k<NFFT_D2; k++)  {
+        Serial.print(k); Serial.print(",");
+        Serial.print(Xre[k],8); Serial.print(",");
+        Serial.print(Xim[k],8); Serial.print(",");
+        Serial.println(P[k],3);
+        }
+    Serial.println();
+
+
+    // And the same data with the CMSIS FFT
+    Serial.println("And now the same thing with CMSIS FFT size NFFT, 0.0 input for imag");
+        // Teensyduino core for T4.x supports arm_cfft_f32
+        // arm_cfft_f32 (const arm_cfft_instance_f32 *S, float32_t *p1, uint8_t ifftFlag, uint8_t bitReverseFlag)
+
+        for(int k=0; k<NFFT; k++)  {
+            fft_buffer[2*k] = x[k];
+            fft_buffer[2*k+1] = 0.0f;
+            }
+        tt = micros();
+        arm_cfft_f32 (&Sfft, fft_buffer, 0, 1);
+        Serial.print("NFFT length FFT uSec = ");
+        Serial.println(micros()-tt);
+        for (int i=0; i < NFFT; i++) {
+            float magsq = fft_buffer[2*i]*fft_buffer[2*i] + fft_buffer[2*i+1]*fft_buffer[2*i+1];
+            Serial.print(i); Serial.print(",");
+            // Serial.print(Xre[k],8); Serial.print(",");
+            // Serial.print(Xim[k],8); Serial.print(",");
+            Serial.println(10.0f*log10f(magsq), 3);
+            }
+
+
+    // And the same data with the CMSIS FFT of size NFFT/2 and interleaved Trig sorting
+    // Time: for NFFT=256, this method is 18 uSec, down from 23 uSec using 256 length
+    // and 0's in the imag inputs.  Memory use is about half and accuracy is the same.
+    Serial.println("And now with CMSIS FFT size NFFT/2, 0.0 input for imag");
+    // Teensyduino core for T4.x supports arm_cfft_f32
+    // arm_cfft_f32 (const arm_cfft_instance_f32 *S, float32_t *p1, uint8_t ifftFlag, uint8_t bitReverseFlag)
+
+    for(int k=0; k<NFFT; k++)  {
+        fft_buffer[k] = x[k];
+        }
+    arm_cfft_f32 (&Sfft_128, fft_buffer, 0, 1);
+
+    // Now the post FT processing for using half-length transform
+    Xre[0] = 0.0f;
+    for(int n=0; n<NFFT; n++)
+       Xre[0] += x[n]/((float)NFFT);    // DC real
+    Xim[0] = 0.0f;                      // DC Imag
+    P[0] = 10.0f*log10f(Xre[0]*Xre[0]);
+    // And the non-DC values
+    for(int i=1; i<NFFT_D2; i++)  {
+        float rns = 0.5f*(fft_buffer[2*i] + fft_buffer[NFFT-2*i]);
+        float ins = 0.5f*(fft_buffer[2*i+1] + fft_buffer[NFFT-2*i+1]);
+        float rnd = 0.5f*(fft_buffer[2*i] - fft_buffer[NFFT-2*i]);
+        float ind = 0.5f*(fft_buffer[2*i+1] - fft_buffer[NFFT-2*i+1]);
+        Xre[i] = rns + cosN[i]*ins - sinN[i]*rnd;
+        Xim[i] = ind - sinN[i]*ins - cosN[i]*rnd;
+        P[i] = 10.0f*log10f(Xre[i]*Xre[i] + Xim[i]*Xim[i]);
+        }
+    for(int k=0; k<NFFT_D2; k++)  {
+        Serial.print(k); Serial.print(",");
+        Serial.print(Xre[k],8); Serial.print(",");
+        Serial.print(Xim[k],8); Serial.print(",");
+        Serial.println(P[k],3);
+        }
+    Serial.println();
+    }
+
+void loop(void)  {
+
+}
+#endif
+
+
+// *****************************************************************
+// *****************************************************************
+
 // TestFFT1024.ino       Bob Larkin  W7PUA
 // Started from PJRC Teensy  Examples/Audio/Analysis/FFT
 //
@@ -16,9 +206,6 @@
 #include <Audio.h>
 #include "OpenAudio_ArduinoLibrary.h"
 
-const int myInput = AUDIO_INPUT_LINEIN;
-//const int myInput = AUDIO_INPUT_MIC;
-
 // Create the Audio components.  These should be created in the
 // order data flows, inputs/sources -> processing -> outputs
 //
@@ -30,7 +217,7 @@ AudioOutputI2S_F32         audioOutput; // audio shield: headphones & line-out N
 // AudioConnection_F32        patchCord1(audioInput, 0, myFFT, 0);
 AudioConnection_F32        patchCord1(sinewave, 0, myFFT, 0);
 AudioControlSGTL5000       audioShield;
-
+float xxx[1024];
 uint32_t ct = 0;
 uint32_t count = 0;
 
@@ -42,15 +229,15 @@ void setup() {
 
   // Enable the audio shield and set the output volume.
   audioShield.enable();
-  audioShield.inputSelect(myInput);
-  audioShield.volume(0.5);
+  audioShield.inputSelect(AUDIO_INPUT_LINEIN);
 
   // Create a synthetic sine wave, for testing
   // To use this, edit the connections above
   // sinewave.frequency(1033.99f);   // Bin 24  T3.x
   // sinewave.frequency(1033.59375f);   // Bin 24  T4.x at 44100
   // sinewave.frequency(1055.0f);  // Bin 24.5, demonstrates windowing
-  sinewave.frequency(1076.0f);
+  // Or some random frequency:
+  sinewave.frequency(1234.5f);
 
   sinewave.amplitude(1.0f);
 
@@ -73,16 +260,19 @@ void setup() {
 
 void loop() {
   if (myFFT.available() /*&& ++ct == 4*/ ) {
-    // each time new FFT data is available
-    // print it all to the Arduino Serial Monitor
-    Serial.println("FFT Output: ");
-    for (int i=0; i<512; i++) {
+     // each time new FFT data is available
+     // print it all to the Arduino Serial Monitor
+
+     float* pin = myFFT.getData();
+     for (int  gg=0; gg<512; gg++)
+        xxx[gg]= *(pin + gg);
+     for (int i=0; i<512; i++) {
         Serial.print(i);
-        Serial.print(",");
-        Serial.println(myFFT.read(i), 3);
+        Serial.print(",  ");
+        Serial.println(xxx[i], 8);
         }
-    Serial.println();
-    }
+     Serial.println();
+     }
 
   /*
   if(count++<200) {
@@ -92,5 +282,5 @@ void loop() {
      Serial.println(AudioMemoryUsageMax_F32());
      }
    */
-  delay(500);
   }
+