wirtz
/
OpenAudio_ArduinoLibrary
mirror of https://github.com/chipaudette/OpenAudio_ArduinoLibrary.git
1
0
Fork 0
Code Issues Releases Wiki Activity

Update real input 1024 fft

Browse Source
pull/6/merge
boblark 3 years ago
parent 7b8cda0bde
commit bc7ad85336
3 changed files with 344 additions and 110 deletions
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Show Stats Download Patch File Download Diff File
  1. 132
      analyze_fft1024_F32.cpp
  2. 100
      analyze_fft1024_F32.h
  3. 222
      examples/TestFFT1024/TestFFT1024.ino

132
analyze_fft1024_F32.cpp
Unescape View File

@ -1,5 +1,6 @@
/* analyze_fft1024_F32.cpp Converted from Teensy I16 Audio Library /* analyze_fft1024_F32.cpp Converted from Teensy I16 Audio Library
* This version uses float F32 inputs. See comments at analyze_fft1024_F32.h * 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 * Conversion parts copyright (c) Bob Larkin 2021
* *
@ -33,36 +34,32 @@
#include "analyze_fft1024_F32.h" #include "analyze_fft1024_F32.h"
// Move audio data from an audio_block_f32_t to the FFT instance buffer. // Move audio data from an audio_block_f32_t to the FFT instance buffer.
// This is for 128 numbers per block // This is for 128 numbers per block, only.
static void copy_to_fft_buffer(void *destination, const void *source) static void copy_to_fft_buffer(void *destination, const void *source) {
{
const float *src = (const float *)source; const float *src = (const float *)source;
float *dst = (float *)destination; float *dst = (float *)destination;
for (int i=0; i < 128; i++) { for (int i=0; i < 128; i++) {
*dst++ = *src++; // real sample *dst++ = *src++; // real sample for half-length FFT
*dst++ = 0.0f; // 0 for Imag }
} }
}
static void apply_window_to_fft_buffer(void *buffer, const void *window) static void apply_window_to_fft_buffer(void *buffer, const void *window) {
{ float *buf = (float *)buffer;
float *buf = (float *)buffer; // 0th entry is real (do window) 1th is imag
const float *win = (float *)window; const float *win = (float *)window;
for (int i=0; i < 1024; i++) for(int i=0; i<NFFT; i++)
buf[2*i] *= *win++; buf[i] *= *win++;
} }
void AudioAnalyzeFFT1024_F32::update(void) { void AudioAnalyzeFFT1024_F32::update(void) {
audio_block_f32_t *block; audio_block_f32_t *block;
uint32_t tt; float outDC=0.0f;
float magsq=0.0f;
block = AudioStream_F32::receiveReadOnly_f32(); block = AudioStream_F32::receiveReadOnly_f32();
if (!block) return; if (!block) return;
// tt=micros();
switch (state) { switch (state) {
case 0: case 0:
blocklist[0] = block; blocklist[0] = block;
@ -82,6 +79,55 @@ void AudioAnalyzeFFT1024_F32::update(void) {
break; break;
case 4: case 4:
blocklist[4] = block; 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; state = 5;
break; break;
case 5: case 5:
@ -96,55 +142,32 @@ void AudioAnalyzeFFT1024_F32::update(void) {
blocklist[7] = block; blocklist[7] = block;
// We have 4 previous blocks pointed to by blocklist[]: // 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+0x000, blocklist[0]->data);
copy_to_fft_buffer(fft_buffer+0x100, blocklist[1]->data); copy_to_fft_buffer(fft_buffer+0x080, blocklist[1]->data);
copy_to_fft_buffer(fft_buffer+0x200, blocklist[2]->data); copy_to_fft_buffer(fft_buffer+0x100, blocklist[2]->data);
copy_to_fft_buffer(fft_buffer+0x300, blocklist[3]->data); copy_to_fft_buffer(fft_buffer+0x180, blocklist[3]->data);
// and 4 new blocks, just gathered: // and 4 new blocks, just gathered:
copy_to_fft_buffer(fft_buffer+0x400, blocklist[4]->data); copy_to_fft_buffer(fft_buffer+0x200, blocklist[4]->data);
copy_to_fft_buffer(fft_buffer+0x500, blocklist[5]->data); copy_to_fft_buffer(fft_buffer+0x280, blocklist[5]->data);
copy_to_fft_buffer(fft_buffer+0x600, blocklist[6]->data); copy_to_fft_buffer(fft_buffer+0x300, blocklist[6]->data);
copy_to_fft_buffer(fft_buffer+0x700, blocklist[7]->data); copy_to_fft_buffer(fft_buffer+0x380, blocklist[7]->data);
if (pWin) if (pWin)
apply_window_to_fft_buffer(fft_buffer, window); 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__) #if defined(__IMXRT1062__)
// Teensyduino core for T4.x supports arm_cfft_f32 // 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 (const arm_cfft_instance_f32 *S, float32_t *p1, uint8_t ifftFlag, uint8_t bitReverseFlag)
arm_cfft_f32 (&Sfft, fft_buffer, 0, 1); arm_cfft_f32 (&Sfft, fft_buffer, 0, 1);
#else #else
// For T3.x go back to old (deprecated) style // For T3.x go back to old (deprecated) style (check radix2/radix4)<<<
arm_cfft_radix4_f32(&fft_inst, fft_buffer); arm_cfft_radix2_f32(&fft_inst, fft_buffer);
#endif #endif
// FFT output is now in fft_buffer. Pick up processing at state==4.
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
AudioStream_F32::release(blocklist[0]); AudioStream_F32::release(blocklist[0]);
AudioStream_F32::release(blocklist[1]); AudioStream_F32::release(blocklist[1]);
@ -158,5 +181,4 @@ void AudioAnalyzeFFT1024_F32::update(void) {
state = 4; state = 4;
break; break;
} // End switch(state) } // End switch(state)
// Serial.print("uSec: "); Serial.println(micros()-tt);
} // End update() } // End update()

100
analyze_fft1024_F32.h
Unescape View File

@ -26,7 +26,7 @@
* THE SOFTWARE. * 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 * 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), * 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 * Power or dBFS (full scale). Output can be bin by bin or a pointer to
@ -45,13 +45,13 @@
* void setOutputType(int _type) * void setOutputType(int _type)
* void setNAverage(int nAverage) * void setNAverage(int nAverage)
* *
* Timing, max is longest update() time: * Timing, max is longest update() time. Comparison is using full complex FFT
* T3.6 Windowed, Power Out, 975 uSec * and no load sharing on "states".
* T3.6 Windowed, RMS out, 1016 uSec * T3.6 Windowed, Power Out, 682 uSec (was 975 w/ 1024 FFT)
* T3.6 Windowed, dBFS out, 1591 uSec * T3.6 Windowed, dBFS out, 834 uSec (was 1591 w/1024 FFT)
* No Window saves 60 uSec on T3.6 for any output. * No Window saves 60 uSec on T3.6 for any output.
* T4.0 Windowed, Ave=1, Power Out, 156 uSec * T4.0 Windowed, Power Out, 54 uSec (was 156 w/1024 FFT)
* T4.0 Windowed, Ave=1, dBFS Out, 302 uSec * T4.0 Windowed, dBFS Out, 203 uSec (was 302 w/1024 FFT)
* Scaling: * Scaling:
* Full scale for floating point DSP is a nebulous concept. Normally the * 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 * 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. * when building the INO.
*/ */
// Fixed float/int problem in read(first, last). RSL 3 Mar 21 // 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_ #ifndef analyze_fft1024_F32_h_
#define analyze_fft1024_F32_h_ #define analyze_fft1024_F32_h_
@ -76,6 +79,15 @@
#include "arm_const_structs.h" #include "arm_const_structs.h"
#endif #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_RMS 0
#define FFT_POWER 1 #define FFT_POWER 1
#define FFT_DBFS 2 #define FFT_DBFS 2
@ -96,13 +108,17 @@ public:
#if defined(__IMXRT1062__) #if defined(__IMXRT1062__)
// Teensy4 core library has the right files for new FFT // Teensy4 core library has the right files for new FFT
// arm CMSIS library has predefined structures of type arm_cfft_instance_f32 // 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 #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 #endif
// This is always 128 block size. Any sample rate. No use of "setings" // This class is always 128 block size. Any sample rate. No use of "settings"
// Bob: Revisit this to use CMSIS fast fft in place of cfft. faster?
useHanningWindow(); 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() // Inform that the output is available for read()
@ -116,7 +132,7 @@ public:
// Output data from a single bin is transferred // Output data from a single bin is transferred
float read(unsigned int binNumber) { 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]; return output[binNumber];
} }
@ -128,8 +144,8 @@ public:
binLast = binFirst; binLast = binFirst;
binFirst = tmp; binFirst = tmp;
} }
if (binFirst > 511) return 0.0; if (binFirst > NFFT_D2M1) return 0.0;
if (binLast > 511) binLast = 511; if (binLast > NFFT_D2M1) binLast = NFFT_D2M1;
float sum = 0.0f; float sum = 0.0f;
do { do {
sum += output[binFirst++]; sum += output[binFirst++];
@ -138,7 +154,7 @@ public:
} }
int windowFunction(int wNum) { int windowFunction(int wNum) {
if(wNum == AudioWindowKaiser1024) if(wNum == AudioWindowKaiser1024) // Changes with size <<<
return -1; // Kaiser needs the kdb return -1; // Kaiser needs the kdb
windowFunction(wNum, 0.0f); windowFunction(wNum, 0.0f);
return 0; return 0;
@ -149,14 +165,14 @@ public:
pWin = window; pWin = window;
if(wNum == NO_WINDOW) if(wNum == NO_WINDOW)
pWin = NULL; pWin = NULL;
else if (wNum == AudioWindowKaiser1024) { else if (wNum == AudioWindowKaiser1024) { // Changes with size <<<
if(_kdb<20.0f) if(_kdb<20.0f)
kd = 20.0f; kd = 20.0f;
else else
kd = _kdb; kd = _kdb;
useKaiserWindow(kd); useKaiserWindow(kd);
} }
else if (wNum == AudioWindowBlackmanHarris1024) else if (wNum == AudioWindowBlackmanHarris1024) // Changes with size <<<
useBHWindow(); useBHWindow();
else else
useHanningWindow(); // Default useHanningWindow(); // Default
@ -177,7 +193,7 @@ public:
// Bring custom window from the INO // Bring custom window from the INO
void putWindow(float *pwin) { void putWindow(float *pwin) {
float *p = window; float *p = window;
for(int i=0; i<1024; i++) for(int i=0; i<NFFT; i++)
*p++ = *pwin++; *p++ = *pwin++;
} }
@ -194,21 +210,27 @@ public:
virtual void update(void); virtual void update(void);
private: private:
float output[512]; float output[NFFT_D2];
float sumsq[512]; float sumsq[NFFT_D2]; // Accumulates averages of outputs
float window[1024]; float window[NFFT];
float *pWin = window; 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]; audio_block_f32_t *blocklist[8];
float fft_buffer[2048];
uint8_t state = 0; uint8_t state = 0;
bool outputflag = false; bool outputflag = false;
audio_block_f32_t *inputQueueArray[1]; audio_block_f32_t *inputQueueArray[1];
#if defined(__IMXRT1062__) #if defined(__IMXRT1062__)
// For T4.x // For T4.x, 512 length for real 1024 input, etc.
// const static arm_cfft_instance_f32 arm_cfft_sR_f32_len1024; // const static arm_cfft_instance_f32 arm_cfft_sR_f32_len512;
arm_cfft_instance_f32 Sfft; arm_cfft_instance_f32 Sfft;
#else #else
arm_cfft_radix4_instance_f32 fft_inst; arm_cfft_radix2_instance_f32 fft_inst; // Check radix2/radix4 <<<
#endif #endif
int outputType = FFT_RMS; //Same type as I16 version init int outputType = FFT_RMS; //Same type as I16 version init
int nAverage = 1; int nAverage = 1;
@ -216,22 +238,22 @@ private:
// The Hann window is a good all-around window // The Hann window is a good all-around window
void useHanningWindow(void) { void useHanningWindow(void) {
for (int i=0; i < 1024; i++) { for (int i=0; i < NFFT; i++) {
// 2*PI/1023 = 0.006141921 float kx = FFT_2PI/((float)NFFT_M1); // 0.006141921 for 1024
window[i] = 0.5*(1.0 - cosf(0.006141921f*(float)i)); window[i] = 0.5f*(1.0f - cosf(kx*(float)i));
} }
} }
// Blackman-Harris produces a first sidelobe more than 90 dB down. // Blackman-Harris produces a first sidelobe more than 90 dB down.
// The price is a bandwidth of about 2 bins. Very useful at times. // The price is a bandwidth of about 2 bins. Very useful at times.
void useBHWindow(void) { void useBHWindow(void) {
for (int i=0; i < 1024; i++) { for (int i=0; i < NFFT; i++) {
float kx = 0.006141921; // 2*PI/1023 float kx = FFT_2PI/((float)NFFT_M1); // 0.006141921 for 1024
int ix = (float) i; int ix = (float) i;
window[i] = 0.35875 - window[i] = 0.35875f -
0.48829*cosf( kx*ix) + 0.48829f*cosf( kx*ix) +
0.14128*cosf(2.0f*kx*ix) - 0.14128f*cosf(2.0f*kx*ix) -
0.01168*cosf(3.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) // 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 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; xn2 = 0.5f+(float32_t)n;
// 4/(1023^2)=0.00000382215877f float kx = 4.0f/(((float)NFFT_M1) * ((float)NFFT_M1)); // 0.00000382215877f for 1024
xn2 = 0.00000382215877f*xn2*xn2; xn2 = kx*xn2*xn2;
window[511 - n]=kbes*(mathEqualizer.i0f(beta*sqrtf(1.0-xn2))); window[NFFT_D2M1 - n]=kbes*(mathEqualizer.i0f(beta*sqrtf(1.0-xn2)));
window[512 + n] = window[511 - n]; window[NFFT_D2 + n] = window[NFFT_D2M1 - n];
} }
} }
}; };

222
examples/TestFFT1024/TestFFT1024.ino
Unescape View File

@ -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 // TestFFT1024.ino Bob Larkin W7PUA
// Started from PJRC Teensy Examples/Audio/Analysis/FFT // Started from PJRC Teensy Examples/Audio/Analysis/FFT
// //
@ -16,9 +206,6 @@
#include <Audio.h> #include <Audio.h>
#include "OpenAudio_ArduinoLibrary.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 // Create the Audio components. These should be created in the
// order data flows, inputs/sources -> processing -> outputs // 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(audioInput, 0, myFFT, 0);
AudioConnection_F32 patchCord1(sinewave, 0, myFFT, 0); AudioConnection_F32 patchCord1(sinewave, 0, myFFT, 0);
AudioControlSGTL5000 audioShield; AudioControlSGTL5000 audioShield;
float xxx[1024];
uint32_t ct = 0; uint32_t ct = 0;
uint32_t count = 0; uint32_t count = 0;
@ -42,15 +229,15 @@ void setup() {
// Enable the audio shield and set the output volume. // Enable the audio shield and set the output volume.
audioShield.enable(); audioShield.enable();
audioShield.inputSelect(myInput); audioShield.inputSelect(AUDIO_INPUT_LINEIN);
audioShield.volume(0.5);
// Create a synthetic sine wave, for testing // Create a synthetic sine wave, for testing
// To use this, edit the connections above // To use this, edit the connections above
// sinewave.frequency(1033.99f); // Bin 24 T3.x // sinewave.frequency(1033.99f); // Bin 24 T3.x
// sinewave.frequency(1033.59375f); // Bin 24 T4.x at 44100 // sinewave.frequency(1033.59375f); // Bin 24 T4.x at 44100
// sinewave.frequency(1055.0f); // Bin 24.5, demonstrates windowing // 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); sinewave.amplitude(1.0f);
@ -73,16 +260,19 @@ void setup() {
void loop() { void loop() {
if (myFFT.available() /*&& ++ct == 4*/ ) { if (myFFT.available() /*&& ++ct == 4*/ ) {
// each time new FFT data is available // each time new FFT data is available
// print it all to the Arduino Serial Monitor // print it all to the Arduino Serial Monitor
Serial.println("FFT Output: ");
for (int i=0; i<512; i++) { 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(i);
Serial.print(","); Serial.print(", ");
Serial.println(myFFT.read(i), 3); Serial.println(xxx[i], 8);
} }
Serial.println(); Serial.println();
} }
/* /*
if(count++<200) { if(count++<200) {
@ -92,5 +282,5 @@ void loop() {
Serial.println(AudioMemoryUsageMax_F32()); Serial.println(AudioMemoryUsageMax_F32());
} }
*/ */
delay(500);
} }

Write Preview
Loading…
Cancel
Save