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

---
 examples/TestFFT1024/TestFFT1024.ino | 239 +++------------------------
 1 file changed, 24 insertions(+), 215 deletions(-)

diff --git a/examples/TestFFT1024/TestFFT1024.ino b/examples/TestFFT1024/TestFFT1024.ino
index 7789232..cb39eb5 100644
--- a/examples/TestFFT1024/TestFFT1024.ino
+++ b/examples/TestFFT1024/TestFFT1024.ino
@@ -1,193 +1,4 @@
-
-#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
 //
@@ -199,7 +10,7 @@ void loop(void)  {
 // the Arduino Serial Monitor.  The format is selectable.
 // Output power averaging is an option
 //
-// T4.0: Uses 11.5% processor and 9 F32 memory blocks, both max.
+// T4.0: Uses 7.8% processor and 9 F32 memory blocks, both max.
 //
 // This example code is in the public domain.
 
@@ -217,15 +28,14 @@ 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;
+
+float saveDat[512];
 
 void setup() {
   Serial.begin(300); // Any speed works
   delay(1000);
 
-  AudioMemory_F32(20);
+  AudioMemory_F32(50);
 
   // Enable the audio shield and set the output volume.
   audioShield.enable();
@@ -256,31 +66,30 @@ void setup() {
   myFFT.setNAverage(1);
 
   myFFT.setOutputType(FFT_DBFS);   // FFT_RMS or FFT_POWER or FFT_DBFS
-}
+
+  Serial.println("1024 point real FFT output in dB relative to full scale sine wave");
+  }
 
 void loop() {
-  if (myFFT.available() /*&& ++ct == 4*/ ) {
+  static uint32_t nTimes = 0;
+  if ( myFFT.available() ) {
      // 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(xxx[i], 8);
-        }
-     Serial.println();
-     }
-
-  /*
-  if(count++<200) {
-     Serial.print("CPU: Max Percent Usage: ");
-     Serial.println(AudioProcessorUsageMax());
-     Serial.print("   Max Float 32 Memory: ");
-     Serial.println(AudioMemoryUsageMax_F32());
+     for (int  kk=0; kk<512; kk++)
+        saveDat[kk]= *(pin + kk);
+     if(++nTimes>4 && nTimes<6)  {
+         for (int i=0; i<512; i++) {
+            Serial.print(i);
+            Serial.print(",  ");
+            Serial.println(saveDat[i], 8);
+            }
+         Serial.println();
+         Serial.print("CPU: Max Percent Usage: ");
+         Serial.println(AudioProcessorUsageMax());
+         Serial.print("   Max Float 32 Memory: ");
+         Serial.println(AudioMemoryUsageMax_F32());
+         }
      }
-   */
   }