FFT with arm math / CMSIS DSP

Hello together,

I am pretty new to this whole topic.

What I would like to achieve:
Basically I would like to achieve this https://learn.adafruit.com/fft-fun-with-fourier-transforms/hardware without the led’s.

So I would like to analyse analog audio signals with a fft.

Current Progress:

I was able to import arm_math and create a arm_cfft_radix4_instance_f32 instance.

#define ARM_MATH_CM3

#undef A0
#undef A1
#undef A2

#include <arm_math.h>

… in loop

arm_cfft_radix4_instance_f32 fft_inst;

My Issue:

If I now try to perform the fft I receive the following error:

    arm_cfft_radix4_init_f32(&fft_inst, FFT_SIZE, ifftFlag, bitReverseFlag);
    arm_cfft_radix4_f32(&fft_inst, samples);
    arm_cmplx_mag_f32(samples, magnitudes, FFT_SIZE);

image

Do you have any idea on how to progress with this issue? I am open for any suggestions!

THANK YOU!

I found some related threads from a few years ago, however they didn’t helped me alot.

https://community.particle.io/t/trying-to-use-cmsis-dsp/9268
https://community.particle.io/t/integrate-cmsis-library-into-firmware/40698
https://community.particle.io/t/fast-fourier-transform-library/3784

Complete coding if required:

Code
/*
 * Project fun-with-fft
 * Description:
 * Author:
 * Date:
 */
#define ARM_MATH_CM3

#undef A0
#undef A1
#undef A2

#include <arm_math.h>


#include "core_cm3.h"

#include "SparkIntervalTimer.h"

int SAMPLE_RATE_HZ = 9000;             // Sample rate of the audio in hertz.
float SPECTRUM_MIN_DB = 30.0;          // Audio intensity (in decibels) that maps to low LED brightness.
float SPECTRUM_MAX_DB = 60.0;          // Audio intensity (in decibels) that maps to high LED brightness.

const uint16_t FFT_SIZE = 256;              // Size of the FFT.  Realistically can only be at most 256
const uint8_t ifftFlag = 0;
const uint8_t bitReverseFlag = 1;
                                       // without running out of memory for buffers and other state.
const int AUDIO_INPUT_PIN = A3;        // Input ADC pin for audio data.
const int ANALOG_READ_RESOLUTION = 10; // Bits of resolution for the ADC.
const int ANALOG_READ_AVERAGING = 16;  // Number of samples to average with each ADC reading.

const int MAX_CHARS = 65;              // Max size of the input command buffer


////////////////////////////////////////////////////////////////////////////////
// INTERNAL STATE
// These shouldn't be modified unless you know what you're doing.
////////////////////////////////////////////////////////////////////////////////

IntervalTimer samplingTimer;
float samples[FFT_SIZE*2];
float magnitudes[FFT_SIZE];
int sampleCounter = 0;
char commandBuffer[MAX_CHARS];



// setup() runs once, when the device is first turned on.
void setup() {

  Serial.begin(38400);
  memset(commandBuffer, 0, sizeof(commandBuffer));


  samplingBegin();

}


void loop() {
  // Calculate FFT if a full sample is available.
  if (samplingIsDone()) {
    // Run FFT on sample data.
    arm_cfft_radix4_instance_f32 fft_inst;
    arm_cfft_radix4_init_f32(&fft_inst, FFT_SIZE, ifftFlag, bitReverseFlag);
    arm_cfft_radix4_f32(&fft_inst, samples);
    arm_cmplx_mag_f32(samples, magnitudes, FFT_SIZE);

    //arm_cfft_f32(&fft_inst, FFT_SIZE, 0, 1)

    //arm_cfft_radix4_f32(&fft_inst, samples);

    // Calculate magnitude of complex numbers output by the FFT.
    //arm_cmplx_mag_f32(samples, magnitudes, FFT_SIZE);

    // Restart audio sampling.
    samplingBegin();
  }

  // Parse any pending commands.
  parserLoop();

}



////////////////////////////////////////////////////////////////////////////////
// UTILITY FUNCTIONS
////////////////////////////////////////////////////////////////////////////////

// Compute the average magnitude of a target frequency window vs. all other frequencies.
void windowMean(float* magnitudes, int lowBin, int highBin, float* windowMean, float* otherMean) {
    *windowMean = 0;
    *otherMean = 0;
    // Notice the first magnitude bin is skipped because it represents the
    // average power of the signal.
    for (int i = 1; i < FFT_SIZE/2; ++i) {
      if (i >= lowBin && i <= highBin) {
        *windowMean += magnitudes[i];
      }
      else {
        *otherMean += magnitudes[i];
      }
    }
    *windowMean /= (highBin - lowBin) + 1;
    *otherMean /= (FFT_SIZE / 2 - (highBin - lowBin));
}

// Convert a frequency to the appropriate FFT bin it will fall within.
int frequencyToBin(float frequency) {
  float binFrequency = float(SAMPLE_RATE_HZ) / float(FFT_SIZE);
  return int(frequency / binFrequency);
}


////////////////////////////////////////////////////////////////////////////////
// SAMPLING FUNCTIONS
////////////////////////////////////////////////////////////////////////////////

void samplingCallback() {
  // Read from the ADC and store the sample data
  samples[sampleCounter] = (float32_t)analogRead(AUDIO_INPUT_PIN);
  // Complex FFT functions require a coefficient for the imaginary part of the input.
  // Since we only have real data, set this coefficient to zero.
  samples[sampleCounter+1] = 0.0;
  // Update sample buffer position and stop after the buffer is filled
  sampleCounter += 2;
  if (sampleCounter >= FFT_SIZE*2) {
    samplingTimer.end();
  }
}

void samplingBegin() {
  // Reset sample buffer position and start callback at necessary rate.
  sampleCounter = 0;
  samplingTimer.begin(samplingCallback, 1000000/SAMPLE_RATE_HZ, uSec);
}

boolean samplingIsDone() {
  return sampleCounter >= FFT_SIZE*2;
}


////////////////////////////////////////////////////////////////////////////////
// COMMAND PARSING FUNCTIONS
// These functions allow parsing simple commands input on the serial port.
// Commands allow reading and writing variables that control the device.
//
// All commands must end with a semicolon character.
//
// Example commands are:
// GET SAMPLE_RATE_HZ;
// - Get the sample rate of the device.
// SET SAMPLE_RATE_HZ 400;
// - Set the sample rate of the device to 400 hertz.
//
////////////////////////////////////////////////////////////////////////////////

void parserLoop() {
  // Process any incoming characters from the serial port
  while (Serial.available() > 0) {
    char c = Serial.read();
    // Add any characters that aren't the end of a command (semicolon) to the input buffer.
    if (c != ';') {
      c = toupper(c);
      strncat(commandBuffer, &c, 1);
    }
    else
    {
      // Parse the command because an end of command token was encountered.
      parseCommand(commandBuffer);
      // Clear the input buffer
      memset(commandBuffer, 0, sizeof(commandBuffer));
    }
  }
}


// Macro used in parseCommand function to simplify parsing get and set commands for a variable
#define GET_AND_SET(variableName) \
  else if (strcmp(command, "GET " #variableName) == 0) { \
    Serial.println(variableName); \
  } \
  else if (strstr(command, "SET " #variableName " ") != NULL) { \
    variableName = (typeof(variableName)) atof(command+(sizeof("SET " #variableName " ")-1)); \
  }

void parseCommand(char* command) {
  if (strcmp(command, "GET MAGNITUDES") == 0) {
    for (int i = 0; i < FFT_SIZE; ++i) {
      Serial.println(magnitudes[i]);
    }
  }
  else if (strcmp(command, "GET SAMPLES") == 0) {
    for (int i = 0; i < FFT_SIZE*2; i+=2) {
      Serial.println(samples[i]);
    }
  }
  else if (strcmp(command, "GET FFT_SIZE") == 0) {
    Serial.println(FFT_SIZE);
  }
  GET_AND_SET(SAMPLE_RATE_HZ)
  GET_AND_SET(SPECTRUM_MIN_DB)
  GET_AND_SET(SPECTRUM_MAX_DB)


}

1 Like

Hi @jan1

You need to go back re-read this thread:

This is a fairly advanced maneuver and involves using a local gcc compile (probably a given already) and changing some makefiles plus using what is called a monolithic or one-part build. The poster at the end of that thread got it working.

With monolithic build you give up the possibility of over-the-air flashing.

This seems like a lot of work to me to get the highly optimized ARM FFT library. Have you considered just using a plain C FFT function instead? There are lots of examples on the web.

1 Like

The PlainFFT library is available for Particle too - not worked with it myself tho’.

Thank you a lot for your fast response @bko and @ScruffR .

It definitely sounds like a lot of work in order to get something working which I probably do not need at all :wink: I think I was too focused to get it working… I wasn’t even considering using a plain C library / function.

So as a next step I will look into the mentioned library and probably do a little bit of research for alternative C libraries.

I will post an update as soon as I get it working (or if I have another issue :wink: ).

If someone already worked with such a lib please let me know.

Again, thank you! Your help is much appreciated.

2 Likes

@jan1 where did your research into FFT libraries land?

Hi everyone,
I don’t wanna open the new thread in this matter I, just want to share with what I achieved with simple arduinoFFT mixed with adafruit example by using the method “copy-paste-try”. My goal was simple make the 3 LED’s reacting just low(red), mid(green), and treb(blue) spectrum, just like old style “disco pod lights”. The final effects are impress me as everybody are talking no is not possible !!! no you can’t analyse audio spectrum without dedicated processor with special super mega DSP and et cetera et cetera … So I tested entire toy with some sin wave gen app and small Sony boombox. and is working just as expected !

What is weird for me, is that all examples on entire internet shows FFT spectrum analiser with sample rates 9000Hz even the adafruit example is just analyser of spectrum from 0-4500Hz and then output results is showing blinking leds but without any sense for me something is blinking there but what exactly I don’t know. I’m able to analyse the spectrum from 50Hz to 18kHz and even if I “pretend” (with my mouth) low (hummmmm), mid(whistling), treb(the rustling)
the leds reacting appropriately.

My Led’s are not neo pixels is my DIY project which I make couple years ago just for learning. On the board I have 3 x AD5241 I2C Compatible, 256-Position Digital Potentiometers , 3 x NE555 as adjustable PWM and 3 x CAT4101 as a power LED’s driver. I do not use any special lib. for AD5241 just simple “registers write” regarding to the data from manual. Photon is communicating via I2C with AD5241 and AD5241 is adjusting NE555 PWM output then NE555 is driving CAT4101 that’s it. Microphone is connected to A0 i used adafruit Electret Microphone Amplifier - MAX9814 with Auto Gain Control

Regarding the code i modify couple things also ArduinoFFT has to be modified a litle and attached manually on web IDE. And I have some questions.

  1. Why we need calculate weird x
    double x = chroma*(1.0-fabs(fmod(h1, 2.0)-1.0));
    when for hue 0, 120.0, 240.0 will always be 0 ?
  2. why we need to spreed windows size evenly ? in this approach we gonna have low mixed with mid and mid mixed with treb and (something is blinking but what is this ? i don’t know)
  3. why my approach is working perfect ? regarding to many examples should’t even work.
    4 why we need special fast “super-duper” FFT lib. ? when is working with the simple one ?

Maybe my code doesn’t make any sense and is working just “by accident” :slight_smile: but my goal is reached !

Here is the code:

// Audio Spectrum Display
// Copyright 2013 Tony DiCola (tony@tonydicola.com)
// borrowed and modified by Arek Otreba
// This code is part of the guide at http://learn.adafruit.com/fft-fun-with-fourier-transforms/


#include "Particle.h"
#include <SparkIntervalTimer.h>
#include "arduinoFFT.h"


void samplingCbISR();                   // forward declaration


////////////////////////////////////////////////////////////////////////////////
// CONIFIGURATION 
// These values can be changed to alter the behavior of the spectrum display.
////////////////////////////////////////////////////////////////////////////////

IntervalTimer samplingTimer;

#define CHANNEL A0                     // Input ADC pin for audio data.
#define INDICATOR_LED_PIN D7           // feedback indicator 
#define blue  0x2C                     // blue LED address 
#define red   0x2D                     // red LED address 
#define green 0x2F                     // green LED address 
#define reg1  0x00                     //  from AD541X manual
#define reg2  0X00
#define AD524X_O1_HIGH  0x10           // AD5241 output O1
#define AD524X_O2_HIGH  0x08           // AD5241 output O2



int SAMPLE_RATE_HZ = 54000;            // Sample rate of the audio in hertz.
float SPECTRUM_MIN_DB = 60.0;          // Audio intensity (in decibels) that maps to low LED brightness.
float SPECTRUM_MAX_DB = 90.0;          // Audio intensity (in decibels) that maps to high LED brightness.                               
const uint16_t FFT_SIZE = 256;         // Size of the FFT.  Realistically can only be at most 256 
int LED_COUNT = 3;                     // Number of Aro pixels.
uint8_t win_type = FFT_WIN_TYP_RECTANGLE;
uint8_t dir = FFT_FORWARD;
double vReal[FFT_SIZE];
double vImag[FFT_SIZE];
int sampleCounter = 0;
double frequencyWindow[4] = {50.0,900.0,5000.0,18000.0};
double hues[3];
bool windowing = true;
double saturation = 1.0;
double low = -2.0;               
double mid = 0.2;
double treb = 0.6;

arduinoFFT FFT = arduinoFFT();         // Create FFT object 

////////////////////////////////////////////////////////////////////////////////
// MAIN SKETCH FUNCTIONS
////////////////////////////////////////////////////////////////////////////////

void setup() {

  SYSTEM_MODE(AUTOMATIC);
  setADCSampleTime(ADC_SampleTime_3Cycles);
  Wire.begin();
  delay(20);
  writeReg(red, reg2, 0);
  writeReg(green, reg2, 0);
  writeReg(blue, reg2, 0);
  Particle.function("CONTROLL", controll);
  Serial.begin(115200);
  Serial.println("Ready");
   
  
  pinMode(INDICATOR_LED_PIN, OUTPUT);
    
  // Initialize spectrum display
  spectrumSetup();
  
  // Begin sampling audio
  samplingBegin();
}

void loop() {
  // Calculate FFT if a full sample is available.
  if (samplingIsDone()) {
    // Run FFT on sample data.
    if(windowing){
        FFT.Windowing(vReal, FFT_SIZE, win_type, dir);
      }
    FFT.Compute(vReal, vImag, FFT_SIZE, FFT_FORWARD); /* Compute FFT */
    // Calculate magnitude of complex numbers output by the FFT.
    FFT.ComplexToMagnitude(vReal, vImag, FFT_SIZE); /* Compute magnitudes */
    
    spectrumLoop();
    
    // Restart audio sampling.
    samplingBegin();
  }
}


////////////////////////////////////////////////////////////////////////////////
// UTILITY FUNCTIONS
////////////////////////////////////////////////////////////////////////////////

// Compute the average magnitude of a target frequency window vs. all other frequencies.
void windowMean(double* vReal, int lowBin, int highBin, double* windowMean, double* otherMean) {
    *windowMean = 0;
    *otherMean = 0;
    // Notice the first magnitude bin is skipped because it represents the
    // average power of the signal.
    for (int i = 1; i < FFT_SIZE/2; ++i) {
      if (i >= lowBin && i <= highBin) {
        *windowMean += vReal[i];
      }
      else {
        *otherMean += vReal[i];
      }
    }
    *windowMean /= (highBin - lowBin) + 1;
    *otherMean /= (FFT_SIZE / 2 - (highBin - lowBin));
}

// Convert a frequency to the appropriate FFT bin it will fall within.
int frequencyToBin(double frequency) {
  double binFrequency = double(SAMPLE_RATE_HZ) / double(FFT_SIZE);
  return int(frequency / binFrequency);
}

// Convert from HSV values (in floating point 0 to 1.0) to RGB colors usable
// by neo pixel functions.
uint32_t pixelHSVtoRGBColor(double hue, double saturation, double value) {
  // Implemented from algorithm at http://en.wikipedia.org/wiki/HSL_and_HSV#From_HSV
  double chroma = value * saturation;
  double h1 = double(hue)/60.0;
  double x = 0.0; //chroma*(1.0-fabs(fmod(h1, 2.0)-1.0)); no need to be calculated for hue 0,120,240 will always be 0
  double r = 0.0;
  double g = 0.0;
  double b = 0.0;
  /*
  Serial.print("HUE");          
  Serial.println(hue, 4);
  Serial.print("CHROMA");          
  Serial.println(chroma, 4);
  Serial.print("WEIRD_X");          
  Serial.println(x, 4);
  
  int R = 0;
  int G = 0;
  int B = 0;
  */

  if (h1 < 1.0) {
    r = chroma;
    g = x;
  }
  else if (h1 < 2.0) {
    r = x;
    g = chroma;
  }
  else if (h1 < 3.0) {
    g = chroma;
    b = x;
  }
  else if (h1 < 4.0) {
    g = x;
    b = chroma;
  }
  else if (h1 < 5.0) {
    r = x;
    b = chroma;
  }
  else // h1 <= 6.0
  {
    r = chroma;
    b = x;
  }
  float m = value - chroma;
  r += m;
  g += m;
  b += m;

 
 if(0.0 >= h1 && h1 <= 2.0){
  if(r == 0){
    writeReg(red, reg2, int(255*r)); 
    }
    
  else{
    writeReg(red, AD524X_O1_HIGH, int(255*r));
    }
 }
 if(2.0 >= h1 && h1 <= 4.0 ){ 
  if(g == 0){
    writeReg(green, reg2, int(255*g)); 
    }
    
  else{
    writeReg(green, AD524X_O1_HIGH, int(255*g));
    }
 }
 
 if(4.0 >= h1 && h1 <= 5.0){ 
   if(b == 0){
    writeReg(blue, reg2, int(255*b));
    }
    
  else{
    writeReg(blue, AD524X_O1_HIGH, int(255*b));
    }
 }
  /*
  Serial.print("R1 ");
  Serial.println(r,3);
  Serial.print("G1 ");
  Serial.println(g, 3);
  Serial.print("B1 ");
  Serial.println(b, 3); 
 */
  
}


////////////////////////////////////////////////////////////////////////////////
// SPECTRUM DISPLAY FUNCTIONS
///////////////////////////////////////////////////////////////////////////////

void spectrumSetup() {
  // Set the frequency window values by evenly dividing the possible frequency
  // spectrum across the number of neo pixels.
  /*
  double windowSize = (SAMPLE_RATE_HZ / 2.0) / double(LED_COUNT);
  for (int i = 0; i < LED_COUNT+1; ++i) {
    frequencyWindow[i] = i*windowSize;
   // making windows size evenly doesn't make any sense for me.
  }
  */
  // Evenly spread hues across all pixels.
  for (int i = 0; i < LED_COUNT; ++i) {
    hues[i] = 360.0*(double(i)/double(LED_COUNT));
  }
}

void spectrumLoop() {
  // Update each LED based on the intensity of the audio 
  // in the associated frequency window.
  double intensity, otherMean;
  for (int i = 0; i < LED_COUNT; ++i) {
    windowMean(vReal, frequencyToBin(frequencyWindow[i]), frequencyToBin(frequencyWindow[i+1]), &intensity, &otherMean);
 // Convert intensity to decibels.
    if(i == 0){
       intensity = (20.0 + low)*log10(intensity);    // adiust intensivity red
    }
    else if(i == 1){
        intensity = (20.0 + mid)*log10(intensity);   // adjust intensivity green
    }
    else if(i == 2){
        intensity = (20.0 + treb)*log10(intensity);  // adjust intensivity blue
    }
    // Scale the intensity and clamp between 0 and 1.0.
    intensity -= SPECTRUM_MIN_DB;
    intensity = intensity < 0.0 ? 0.0 : intensity;
    intensity /= (SPECTRUM_MAX_DB-SPECTRUM_MIN_DB);
    intensity = intensity > 1.0 ? 1.0 : intensity;
    pixelHSVtoRGBColor(hues[i], saturation, intensity);
  }
  
}

////////////////////////////////////////////////////////////////////////////////
// communication with AD524X
///////////////////////////////////////////////////////////////////////////////

void writeReg(byte device, byte reg, int Data)
{
    Wire.beginTransmission(device);
    Wire.write(reg);
    Wire.write(Data);
  
  Wire.endTransmission();
}

////////////////////////////////////////////////////////////////////////////////
// Paticle controll function 
///////////////////////////////////////////////////////////////////////////////


int controll(String message){
    // ex: setrelais:on
    // ex: setrelais:off
    // ex: setalarm:10
    int colonPos = message.indexOf(":") ;
    if (colonPos < 0) return -1 ;                 // syntax error not found
    String command = message.substring(0,colonPos) ;
    String argument = message.substring(colonPos+1) ;
    
    if(command.equals("windowing")){
        int pointer = argument.toInt();
        if(pointer == 0){
          windowing = false;
        }
        if(pointer == 1){
          windowing = true;
        }
       
       return windowing;
    }
    else if(command.equals("samplehz")){
        int pointer = argument.toInt();
        SAMPLE_RATE_HZ = pointer;
        return SAMPLE_RATE_HZ;
        
    }
    else if(command.equals("spectrummin")){
        double tempolary = argument.toFloat();
        SPECTRUM_MIN_DB = tempolary;
        return SPECTRUM_MIN_DB;
    }
    
    else if (command.equals("spectrummax")){
        double tempolary = argument.toFloat();
        SPECTRUM_MAX_DB = tempolary;
        
        return SPECTRUM_MAX_DB;
    }  
    
    else if (command.equals("saturation")){
        double tempolary = argument.toFloat();
        saturation = tempolary;
        return saturation;
    } 
    else if (command.equals("low")){
        double tempolary = argument.toFloat();
        low = tempolary;
        return low;
    }
    else if (command.equals("mid")){
        double tempolary = argument.toFloat();
        mid = tempolary;
        return mid;
    }
    else if (command.equals("treb")){
        double tempolary = argument.toFloat();
        treb = tempolary;
        return treb;
    }
    
    else if (command.equals("ledcount")){
        int pointer = argument.toInt();
        LED_COUNT = pointer;
        return LED_COUNT;
    } 
    else if (command.equals("wintype")){
      int pointer = argument.toInt();
      if(pointer == 0){
         win_type = FFT_WIN_TYP_WELCH;  /* welch */
         return FFT_WIN_TYP_WELCH;
         }
      if(pointer == 1 ){
         win_type = FFT_WIN_TYP_RECTANGLE;  /* rectangle (Box car) */
         return FFT_WIN_TYP_RECTANGLE;
      }
      if(pointer == 2){
        win_type = FFT_WIN_TYP_HAMMING;  /* hamming */
        return FFT_WIN_TYP_HAMMING;
      }
      if(pointer == 3){
        win_type = FFT_WIN_TYP_HANN;  /* hann */
        return FFT_WIN_TYP_HANN;
      }
      if(pointer == 4){
        win_type = FFT_WIN_TYP_TRIANGLE; /* triangle (Bartlett) */
        return FFT_WIN_TYP_TRIANGLE;
      }
      if(pointer == 5){
        win_type = FFT_WIN_TYP_NUTTALL;  /* nuttall */
        return FFT_WIN_TYP_NUTTALL;
      }
      if(pointer == 6){
         win_type =FFT_WIN_TYP_BLACKMAN;  /* blackman */
         return FFT_WIN_TYP_BLACKMAN;
      }
      if(pointer == 7){
        win_type = FFT_WIN_TYP_BLACKMAN_NUTTALL;  /* blackman nuttall */
        return FFT_WIN_TYP_BLACKMAN_NUTTALL;
      }
      if(pointer == 8){
         win_type = FFT_WIN_TYP_BLACKMAN_HARRIS;  /* blackman harris*/
         return FFT_WIN_TYP_BLACKMAN_HARRIS;
      }
      if(pointer == 9){
         win_type = FFT_WIN_TYP_FLT_TOP;  /* flat top */
         return FFT_WIN_TYP_FLT_TOP;
      }
      
      }
      else if (command.equals("dir")){
      int pointer = argument.toInt();
       if(pointer == 0){
         dir = FFT_REVERSE;  
         return FFT_REVERSE;
         }
      if(pointer == 1 ){
         dir = FFT_FORWARD;
         return FFT_FORWARD;
         }
      }
    
     else if (command.equalsIgnoreCase("reset")){
        System.reset();
        return 0;
    } else {
        return -2 ;                         // command not found
    }

}

////////////////////////////////////////////////////////////////////////////////
// SAMPLING FUNCTIONS
////////////////////////////////////////////////////////////////////////////////

void samplingCbISR() {
 
  // Read from the ADC and store the sample data
  vReal[sampleCounter] =  analogRead(CHANNEL);
  // Complex FFT functions require a coefficient for the imaginary part of the input.
  // Since we only have real data, set this coefficient to zero.
  vImag[sampleCounter] = 0.0;
  // Update sample buffer position and stop after the buffer is filled
  sampleCounter += 1;
  if (sampleCounter >= FFT_SIZE) {
    samplingTimer.end();
    digitalWrite(INDICATOR_LED_PIN, HIGH);
  }
}

void samplingBegin() {
  // Reset sample buffer position and start callback at necessary rate.
  digitalWrite(INDICATOR_LED_PIN, LOW);
  sampleCounter = 0;
  samplingTimer.begin(samplingCbISR, 1000000/SAMPLE_RATE_HZ, uSec);
}

boolean samplingIsDone() {
  return sampleCounter >= FFT_SIZE;
}

and modified ArduinoFFT.cpp:



#include "arduinoFFT.h"

arduinoFFT::arduinoFFT(void)
{ // Constructor
	#warning("This method is deprecated and will be removed on future revisions.")
}

arduinoFFT::arduinoFFT(double *vReal, double *vImag, uint16_t samples, double samplingFrequency)
{// Constructor
	this->_vReal = vReal;
	this->_vImag = vImag;
	this->_samples = samples;
	this->_samplingFrequency = samplingFrequency;
	this->_power = Exponent(samples);
}

arduinoFFT::~arduinoFFT(void)
{
// Destructor
}

uint8_t arduinoFFT::Revision(void)
{
	return(FFT_LIB_REV);
}

void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t dir)
{
	#warning("This method is deprecated and will be removed on future revisions.")
	Compute(vReal, vImag, samples, Exponent(samples), dir);
}

void arduinoFFT::Compute(uint8_t dir)
{// Computes in-place complex-to-complex FFT /
	// Reverse bits /
	uint16_t j = 0;
	for (uint16_t i = 0; i < (this->_samples - 1); i++) {
		if (i < j) {
			Swap(&this->_vReal[i], &this->_vReal[j]);
			if(dir==FFT_REVERSE)
				Swap(&this->_vImag[i], &this->_vImag[j]);
		}
		uint16_t k = (this->_samples >> 1);
		while (k <= j) {
			j -= k;
			k >>= 1;
		}
		j += k;
	}
	// Compute the FFT  /
	double c1 = -1.0;
	double c2 = 0.0;
	uint16_t l2 = 1;
	for (uint8_t l = 0; (l < this->_power); l++) {
		uint16_t l1 = l2;
		l2 <<= 1;
		double u1 = 1.0;
		double u2 = 0.0;
		for (j = 0; j < l1; j++) {
			 for (uint16_t i = j; i < this->_samples; i += l2) {
					uint16_t i1 = i + l1;
					double t1 = u1 * this->_vReal[i1] - u2 * this->_vImag[i1];
					double t2 = u1 * this->_vImag[i1] + u2 * this->_vReal[i1];
					this->_vReal[i1] = this->_vReal[i] - t1;
					this->_vImag[i1] = this->_vImag[i] - t2;
					this->_vReal[i] += t1;
					this->_vImag[i] += t2;
			 }
			 double z = ((u1 * c1) - (u2 * c2));
			 u2 = ((u1 * c2) + (u2 * c1));
			 u1 = z;
		}
		c2 = sqrt((1.0 - c1) / 2.0);
		if (dir == FFT_FORWARD) {
			c2 = -c2;
		}
		c1 = sqrt((1.0 + c1) / 2.0);
	}
	// Scaling for reverse transform /
	if (dir != FFT_FORWARD) {
		for (uint16_t i = 0; i < this->_samples; i++) {
			 this->_vReal[i] /= this->_samples;
			 this->_vImag[i] /= this->_samples;
		}
	}
}

void arduinoFFT::Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir)
{	// Computes in-place complex-to-complex FFT
	// Reverse bits
	#warning("This method is deprecated and will be removed on future revisions.")
	uint16_t j = 0;
	for (uint16_t i = 0; i < (samples - 1); i++) {
		if (i < j) {
			Swap(&vReal[i], &vReal[j]);
			if(dir==FFT_REVERSE)
				Swap(&vImag[i], &vImag[j]);
		}
		uint16_t k = (samples >> 1);
		while (k <= j) {
			j -= k;
			k >>= 1;
		}
		j += k;
	}
	// Compute the FFT
	double c1 = -1.0;
	double c2 = 0.0;
	uint16_t l2 = 1;
	for (uint8_t l = 0; (l < power); l++) {
		uint16_t l1 = l2;
		l2 <<= 1;
		double u1 = 1.0;
		double u2 = 0.0;
		for (j = 0; j < l1; j++) {
			 for (uint16_t i = j; i < samples; i += l2) {
					uint16_t i1 = i + l1;
					double t1 = u1 * vReal[i1] - u2 * vImag[i1];
					double t2 = u1 * vImag[i1] + u2 * vReal[i1];
					vReal[i1] = vReal[i] - t1;
					vImag[i1] = vImag[i] - t2;
					vReal[i] += t1;
					vImag[i] += t2;
			 }
			 double z = ((u1 * c1) - (u2 * c2));
			 u2 = ((u1 * c2) + (u2 * c1));
			 u1 = z;
		}
		c2 = sqrt((1.0 - c1) / 2.0);
		if (dir == FFT_FORWARD) {
			c2 = -c2;
		}
		c1 = sqrt((1.0 + c1) / 2.0);
	}
	// Scaling for reverse transform
	if (dir != FFT_FORWARD) {
		for (uint16_t i = 0; i < samples; i++) {
			 vReal[i] /= samples;
			 vImag[i] /= samples;
		}
	}
}

void arduinoFFT::ComplexToMagnitude()
{ // vM is half the size of vReal and vImag
	for (uint16_t i = 0; i < this->_samples; i++) {
		this->_vReal[i] = sqrt(sq(this->_vReal[i]) + sq(this->_vImag[i]));
	}
}

void arduinoFFT::ComplexToMagnitude(double *vReal, double *vImag, uint16_t samples)
{	// vM is half the size of vReal and vImag
	#warning("This method is deprecated and will be removed on future revisions.")
	for (uint16_t i = 0; i < samples; i++) {
		vReal[i] = sqrt(sq(vReal[i]) + sq(vImag[i]));
	}
}

void arduinoFFT::DCRemoval()
{
	// calculate the mean of vData
	double mean = 0;
	for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++)
	{
		mean += this->_vReal[i];
	}
	mean /= this->_samples;
	// Subtract the mean from vData
	for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++)
	{
		this->_vReal[i] -= mean;
	}
}

void arduinoFFT::DCRemoval(double *vData, uint16_t samples)
{
	// calculate the mean of vData
	//#warning("This method is deprecated and will be removed on future revisions.")
	double mean = 0;
	for (uint16_t i = 1; i < ((samples >> 1) + 1); i++)
	{
		mean += vData[i];
	}
	mean /= samples;
	// Subtract the mean from vData
	for (uint16_t i = 1; i < ((samples >> 1) + 1); i++)
	{
		vData[i] -= mean;
	}
}

void arduinoFFT::Windowing(uint8_t windowType, uint8_t dir)
{// Weighing factors are computed once before multiple use of FFT
// The weighing function is symetric; half the weighs are recorded
	double samplesMinusOne = (double(this->_samples) - 1.0);
	for (uint16_t i = 0; i < (this->_samples >> 1); i++) {
		double indexMinusOne = double(i);
		double ratio = (indexMinusOne / samplesMinusOne);
		double weighingFactor = 1.0;
		// Compute and record weighting factor
		switch (windowType) {
		case FFT_WIN_TYP_RECTANGLE: // rectangle (box car)
			weighingFactor = 1.0;
			break;
		case FFT_WIN_TYP_HAMMING: // hamming
			weighingFactor = 0.54 - (0.46 * cos(twoPi * ratio));
			break;
		case FFT_WIN_TYP_HANN: // hann
			weighingFactor = 0.54 * (1.0 - cos(twoPi * ratio));
			break;
		case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett)
			weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne);
			break;
		case FFT_WIN_TYP_NUTTALL: // nuttall
			weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) + (0.144232 * (cos(fourPi * ratio))) - (0.012604 * (cos(sixPi * ratio)));
			break;
		case FFT_WIN_TYP_BLACKMAN: // blackman
			weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) + (0.07922 * (cos(fourPi * ratio)));
			break;
		case FFT_WIN_TYP_BLACKMAN_NUTTALL: // blackman nuttall
			weighingFactor = 0.3635819 - (0.4891775 * (cos(twoPi * ratio))) + (0.1365995 * (cos(fourPi * ratio))) - (0.0106411 * (cos(sixPi * ratio)));
			break;
		case FFT_WIN_TYP_BLACKMAN_HARRIS: // blackman harris
			weighingFactor = 0.35875 - (0.48829 * (cos(twoPi * ratio))) + (0.14128 * (cos(fourPi * ratio))) - (0.01168 * (cos(sixPi * ratio)));
			break;
		case FFT_WIN_TYP_FLT_TOP: // flat top
			weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) + (0.1980399 * cos(fourPi * ratio));
			break;
		case FFT_WIN_TYP_WELCH: // welch
			weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0));
			break;
		}
		if (dir == FFT_FORWARD) {
			this->_vReal[i] *= weighingFactor;
			this->_vReal[this->_samples - (i + 1)] *= weighingFactor;
		}
		else {
			this->_vReal[i] /= weighingFactor;
			this->_vReal[this->_samples - (i + 1)] /= weighingFactor;
		}
	}
}


void arduinoFFT::Windowing(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir)
{// Weighing factors are computed once before multiple use of FFT
// The weighing function is symetric; half the weighs are recorded
	#warning("This method is deprecated and will be removed on future revisions.")
	double samplesMinusOne = (double(samples) - 1.0);
	for (uint16_t i = 0; i < (samples >> 1); i++) {
		double indexMinusOne = double(i);
		double ratio = (indexMinusOne / samplesMinusOne);
		double weighingFactor = 1.0;
		// Compute and record weighting factor
		switch (windowType) {
		case FFT_WIN_TYP_RECTANGLE: // rectangle (box car)
			weighingFactor = 1.0;
			break;
		case FFT_WIN_TYP_HAMMING: // hamming
			weighingFactor = 0.54 - (0.46 * cos(twoPi * ratio));
			break;
		case FFT_WIN_TYP_HANN: // hann
			weighingFactor = 0.54 * (1.0 - cos(twoPi * ratio));
			break;
		case FFT_WIN_TYP_TRIANGLE: // triangle (Bartlett)
			weighingFactor = 1.0 - ((2.0 * abs(indexMinusOne - (samplesMinusOne / 2.0))) / samplesMinusOne);
			break;
		case FFT_WIN_TYP_NUTTALL: // nuttall
			weighingFactor = 0.355768 - (0.487396 * (cos(twoPi * ratio))) + (0.144232 * (cos(fourPi * ratio))) - (0.012604 * (cos(sixPi * ratio)));
			break;
		case FFT_WIN_TYP_BLACKMAN: // blackman
			weighingFactor = 0.42323 - (0.49755 * (cos(twoPi * ratio))) + (0.07922 * (cos(fourPi * ratio)));
			break;
		case FFT_WIN_TYP_BLACKMAN_NUTTALL: // blackman nuttall
			weighingFactor = 0.3635819 - (0.4891775 * (cos(twoPi * ratio))) + (0.1365995 * (cos(fourPi * ratio))) - (0.0106411 * (cos(sixPi * ratio)));
			break;
		case FFT_WIN_TYP_BLACKMAN_HARRIS: // blackman harris
			weighingFactor = 0.35875 - (0.48829 * (cos(twoPi * ratio))) + (0.14128 * (cos(fourPi * ratio))) - (0.01168 * (cos(sixPi * ratio)));
			break;
		case FFT_WIN_TYP_FLT_TOP: // flat top
			weighingFactor = 0.2810639 - (0.5208972 * cos(twoPi * ratio)) + (0.1980399 * cos(fourPi * ratio));
			break;
		case FFT_WIN_TYP_WELCH: // welch
			weighingFactor = 1.0 - sq((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0));
			break;
		}
		if (dir == FFT_FORWARD) {
			vData[i] *= weighingFactor;
			vData[samples - (i + 1)] *= weighingFactor;
		}
		else {
			vData[i] /= weighingFactor;
			vData[samples - (i + 1)] /= weighingFactor;
		}
	}
}

double arduinoFFT::MajorPeak()
{
	double maxY = 0;
	uint16_t IndexOfMaxY = 0;
	//If sampling_frequency = 2 * max_frequency in signal,
	//value would be stored at position samples/2
	for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) {
		if ((this->_vReal[i-1] < this->_vReal[i]) && (this->_vReal[i] > this->_vReal[i+1])) {
			if (this->_vReal[i] > maxY) {
				maxY = this->_vReal[i];
				IndexOfMaxY = i;
			}
		}
	}
	double delta = 0.5 * ((this->_vReal[IndexOfMaxY-1] - this->_vReal[IndexOfMaxY+1]) / (this->_vReal[IndexOfMaxY-1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY+1]));
	double interpolatedX = ((IndexOfMaxY + delta)  * this->_samplingFrequency) / (this->_samples-1);
	if(IndexOfMaxY==(this->_samples >> 1)) //To improve calculation on edge values
		interpolatedX = ((IndexOfMaxY + delta)  * this->_samplingFrequency) / (this->_samples);
	// returned value: interpolated frequency peak apex
	return(interpolatedX);
}

void arduinoFFT::MajorPeak(double *f, double *v)
{
	double maxY = 0;
	uint16_t IndexOfMaxY = 0;
	//If sampling_frequency = 2 * max_frequency in signal,
	//value would be stored at position samples/2
	for (uint16_t i = 1; i < ((this->_samples >> 1) + 1); i++) {
		if ((this->_vReal[i - 1] < this->_vReal[i]) && (this->_vReal[i] > this->_vReal[i + 1])) {
			if (this->_vReal[i] > maxY) {
				maxY = this->_vReal[i];
				IndexOfMaxY = i;
			}
		}
	}
	double delta = 0.5 * ((this->_vReal[IndexOfMaxY - 1] - this->_vReal[IndexOfMaxY + 1]) / (this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]));
	double interpolatedX = ((IndexOfMaxY + delta)  * this->_samplingFrequency) / (this->_samples - 1);
	if (IndexOfMaxY == (this->_samples >> 1)) //To improve calculation on edge values
		interpolatedX = ((IndexOfMaxY + delta)  * this->_samplingFrequency) / (this->_samples);
	// returned value: interpolated frequency peak apex
	*f = interpolatedX;
	*v = abs(this->_vReal[IndexOfMaxY - 1] - (2.0 * this->_vReal[IndexOfMaxY]) + this->_vReal[IndexOfMaxY + 1]);
}

double arduinoFFT::MajorPeak(double *vD, uint16_t samples, double samplingFrequency)
{
	#warning("This method is deprecated and will be removed on future revisions.")
	double maxY = 0;
	uint16_t IndexOfMaxY = 0;
	//If sampling_frequency = 2 * max_frequency in signal,
	//value would be stored at position samples/2
	for (uint16_t i = 1; i < ((samples >> 1) + 1); i++) {
		if ((vD[i-1] < vD[i]) && (vD[i] > vD[i+1])) {
			if (vD[i] > maxY) {
				maxY = vD[i];
				IndexOfMaxY = i;
			}
		}
	}
	double delta = 0.5 * ((vD[IndexOfMaxY-1] - vD[IndexOfMaxY+1]) / (vD[IndexOfMaxY-1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY+1]));
	double interpolatedX = ((IndexOfMaxY + delta)  * samplingFrequency) / (samples-1);
	if(IndexOfMaxY==(samples >> 1)) //To improve calculation on edge values
		interpolatedX = ((IndexOfMaxY + delta)  * samplingFrequency) / (samples);
	// returned value: interpolated frequency peak apex
	return(interpolatedX);
}

void arduinoFFT::MajorPeak(double *vD, uint16_t samples, double samplingFrequency, double *f, double *v)
{
	#warning("This method is deprecated and will be removed on future revisions.")
	double maxY = 0;
	uint16_t IndexOfMaxY = 0;
	//If sampling_frequency = 2 * max_frequency in signal,
	//value would be stored at position samples/2
	for (uint16_t i = 1; i < ((samples >> 1) + 1); i++) {
		if ((vD[i - 1] < vD[i]) && (vD[i] > vD[i + 1])) {
			if (vD[i] > maxY) {
				maxY = vD[i];
				IndexOfMaxY = i;
			}
		}
	}
	double delta = 0.5 * ((vD[IndexOfMaxY - 1] - vD[IndexOfMaxY + 1]) / (vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]));
	double interpolatedX = ((IndexOfMaxY + delta)  * samplingFrequency) / (samples - 1);
	//double popo =
	if (IndexOfMaxY == (samples >> 1)) //To improve calculation on edge values
		interpolatedX = ((IndexOfMaxY + delta)  * samplingFrequency) / (samples);
	// returned value: interpolated frequency peak apex
	*f = interpolatedX;
	*v = abs(vD[IndexOfMaxY - 1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY + 1]);
}

uint8_t arduinoFFT::Exponent(uint16_t value)
{
	#warning("This method will not be accessible on future revisions.")
	// Calculates the base 2 logarithm of a value
	uint8_t result = 0;
	while (((value >> result) & 1) != 1) result++;
	return(result);
}

// Private functions

void arduinoFFT::Swap(double *x, double *y)
{
	double temp = *x;
	*x = *y;
	*y = temp;
}

and modified ArduinoFFT.h:



#ifndef arduinoFFT_h /* Prevent loading library twice */
#define arduinoFFT_h
	#include <stdlib.h>
	#include <stdio.h>
	#ifdef __AVR__
		#include <avr/io.h>
	#endif
	#include <math.h>
	//#include "defs.h"
	//#include "types.h"
//#endif

#define FFT_LIB_REV 0x14
/* Custom constants */
#define FFT_FORWARD 0x01
#define FFT_REVERSE 0x00

/* Windowing type */
#define FFT_WIN_TYP_RECTANGLE 0x00 /* rectangle (Box car) */
#define FFT_WIN_TYP_HAMMING 0x01 /* hamming */
#define FFT_WIN_TYP_HANN 0x02 /* hann */
#define FFT_WIN_TYP_TRIANGLE 0x03 /* triangle (Bartlett) */
#define FFT_WIN_TYP_NUTTALL 0x04 /* nuttall */
#define FFT_WIN_TYP_BLACKMAN 0x05 /* blackman */
#define FFT_WIN_TYP_BLACKMAN_NUTTALL 0x06 /* blackman nuttall */
#define FFT_WIN_TYP_BLACKMAN_HARRIS 0x07 /* blackman harris*/
#define FFT_WIN_TYP_FLT_TOP 0x08 /* flat top */
#define FFT_WIN_TYP_WELCH 0x09 /* welch */
/*Mathematial constants*/
#define twoPi 6.28318531
#define fourPi 12.56637061
#define sixPi 18.84955593
#define sq(x) ((x)*(x))

class arduinoFFT {
public:
	/* Constructor */
	arduinoFFT(void);
	arduinoFFT(double *vReal, double *vImag, uint16_t samples, double samplingFrequency);
	/* Destructor */
	~arduinoFFT(void);
	/* Functions */
	uint8_t Revision(void);
	uint8_t Exponent(uint16_t value);
	void ComplexToMagnitude(double *vReal, double *vImag, uint16_t samples);
	void Compute(double *vReal, double *vImag, uint16_t samples, uint8_t dir);
	void Compute(double *vReal, double *vImag, uint16_t samples, uint8_t power, uint8_t dir);
	void DCRemoval(double *vData, uint16_t samples);
	double MajorPeak(double *vD, uint16_t samples, double samplingFrequency);
	void Windowing(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir);
	void ComplexToMagnitude();
	void Compute(uint8_t dir);
	void DCRemoval();
	double MajorPeak();
	void Windowing(uint8_t windowType, uint8_t dir);

	void MajorPeak(double *f, double *v);
	void MajorPeak(double *vD, uint16_t samples, double samplingFrequency, double *f, double *v);


private:
	/* Variables */
	uint16_t _samples;
	double _samplingFrequency;
	double *_vReal;
	double *_vImag;
	uint8_t _power;
	/* Functions */
	void Swap(double *x, double *y);
};

#endif

http://www.ti.com/lit/ds/symlink/na555.pdf

5 Likes

I made some video is annoying i know :slight_smile: but shows that all is working pretty good. this was just to check if will work at all but the main purpose is that I would like to try use this on the floor at work to control bearings and belts splices on conveyors and all weird noises to prevent down time. All data can can be posted to Particle cloud. It will be great if SparkIntervalTimer.h will work on xenon’s as I would like to try make a mesh network and attach as many xenon’s as possible because we have a lot of conveyors.

3 Likes