Update real input 1024 fft

pull/6/merge
boblark 3 years ago
parent bc7ad85336
commit bc44235e00
  1. 239
      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());
}
}
*/
}

Loading…
Cancel
Save