Fft

frohr

Well-known member
Hi all,
I want read ADXL from analog pin, for example 16384 samples and save values to buffer. Then I want apply FFT (hanning window) and get average value. But I think I have problem because Arduino library FFT has only Magnitude and no amplitude. So I have wrong data to get mV or mm/s.
Do you have any suggestion how make FFT with amplitude output for example in mV?
Thanks!
Michal

Code:
void rms()
{ 
   suma = 0;
   microseconds = micros();
  for(int i=0; i<BUFFER_SIZE; i++)
  {
      vReal[i] = analogRead(A1);
      vImag[i] = 0;
      while(micros() - microseconds < sampling_period_us){
      }
      microseconds += sampling_period_us;
  }
  
  FFT.Windowing(vReal, BUFFER_SIZE, FFT_WIN_TYP_HANN, FFT_FORWARD);
  FFT.Compute(vReal, vImag, BUFFER_SIZE, FFT_FORWARD);
  FFT.ComplexToMagnitude(vReal, vImag, BUFFER_SIZE);  

   for(int i= 0;i < BUFFER_SIZE - 1; i++)    
   {
    frequency = (i * 1.0 * 20000) / BUFFER_SIZE;    
    magnitude = vReal[i];  
    suma = suma + magnitude;
   }  
   suma = (suma / BUFFER_SIZE);
   Serial.println(suma, 0);
 
To me, the terms magnitude and amplitude, as applied to Fourier coefficients, mean the same thing. What do they mean to you?

Code:
mag[i] = sqrt(vReal[i]*vReal[i] + vImag[i]*vImag[i]);
 
I use Teensy 4

What I need:
Get values in g or mV after FFT windowing
=
Read 16384 values (samples) from ADXL to buffer
FFT, hanning windowing
Calculate frequency for every bin from FFT
Calculate gRMS and mm/s from FFT
But the FFT has just Magnitude. And magnitude is for me dimensionless value - maybe I am wrong. Amplitude is in mV or in g and I will be able to calculate for example mm/s for 10-1000Hz - after band pass...
 
multiply FFT value with sensitivity and divide by i*omega to integrate accelerometer to get particle velocity
 
I try this code now, but I do not know what number is SUMA - I need mV or g...

and signal is not uniform.


Code:
void rms_mm()
{ 
  suma = 0;
  microseconds = micros();
  for(int i=0; i<BUFFER_SIZE; i++)
  {
      vReal[i] = analogRead(A1);
      vImag[i] = 0;
      while(micros() - microseconds < sampling_period_us)
      {
      }
      microseconds += sampling_period_us;
  }

  F_LOW_F();

 for(int i= 0;i < BUFFER_SIZE; i++)
   {  
    micros_0 = 0;  
    Filter_LOWPASS =  (vReal[i] *alpha_LOW) + (Filter_LOWPASS * (1.0 - alpha_LOW));
    Filter_HIGHPASS = (vReal[i] *alpha_HIGH) + (Filter_HIGHPASS * (1.0 - alpha_HIGH));
    Filter_BANDPASS =  Filter_HIGHPASS-Filter_LOWPASS;        
    sensorValueG = Filter_BANDPASS; 
    vReal[i] =    sensorValueG;    
    vImag[i] = period *i;
    while(micros_0 <=  period) 
    {      
    }     
   }  
    
  FFT.Windowing(vReal, BUFFER_SIZE, FFT_WIN_TYP_HANN, FFT_FORWARD);
  FFT.Compute(vReal, vImag, BUFFER_SIZE, FFT_FORWARD);
  FFT.ComplexToMagnitude(vReal, vImag, BUFFER_SIZE);

  for(int i= 4;i < BUFFER_SIZE / 20; i++)    //2 pro 10000Hz /20 pro 1000Hz
   {
    vReal[i] = pow(vReal[i], 2);
    frekvence = (i * 1.0 * 20000) / BUFFER_SIZE;    
    magnitude = vReal[i];  
    suma = suma + magnitude;
   }
   suma = sqrt(suma);
   Serial.println(suma);
  
}

sensor no movement.jpg
sensor motor started.jpg
 
Measuring bearings vibrations in mm/s.
Now I am able measure g RMS for 500-10000 Hz from time signal - confirmed with professional tools.

But for mm/s I have to use FFT.

As you know accelerometer give you length/sec^2 (acceleration). If you wanted velocity (mm/s in your case) you need to integrate your data
The so-called spectral method does not integrate the time-series but the the Fourier spectrum which is equivalent to divide the spectrum by i*omega.
So what you get is the spectral particle velocity (in your case vibration)
To get the spectrum correctly estimated, you may wanted to check that the sum of all spectra squared equals to sum of time series squared (energy conservation; Parseval theorem)
after that you divide spectrum by i*omega

Edit: again, where is FFT class defined?
 
Last edited:
I am just a simple mathematician, and I know that you can integrate, at least formally since Fourier techniques makes assumptions that are not always satisfied by measured signals, by dividing the Fourier transform by i omega. The units and how to convert them can be seen from the formulas for FFT and iFFT.

Personally I would start by simply making numerical integration, no FFT and calculate the RMS amplitude of this signal as you do with the acceleration. This of course assumes you simply want the RMS and not the spectral components.
 
I am just a simple mathematician, and I know that you can integrate, at least formally since Fourier techniques makes assumptions that are not always satisfied by measured signals, by dividing the Fourier transform by i omega.

implicit periodicity due to FFT is one of these assumptions that is not always satisfied.
But, some scientific areas (seismology, vibration testing, acoustics, etc) like to work in spectral domain.
 
As you know accelerometer give you length/sec^2 (acceleration). If you wanted velocity (mm/s in your case) you need to integrate your data
The so-called spectral method does not integrate the time-series but the the Fourier spectrum which is equivalent to divide the spectrum by i*omega.
So what you get is the spectral particle velocity (in your case vibration)
To get the spectrum correctly estimated, you may wanted to check that the sum of all spectra squared equals to sum of time series squared (energy conservation; Parseval theorem)
after that you divide spectrum by i*omega

Edit: again, where is FFT class defined?


Code:
#include <ADC.h>
#include <ADC_util.h>

#include "arduinoFFT.h"
arduinoFFT FFT = arduinoFFT();

IntervalTimer timer;
elapsedMicros micros_0;
float cas; 
float period; 
float alpha_LOW = 0;
float alpha_HIGH = 0;
double Filter_LOWPASS = 0;
double Filter_HIGHPASS = 0;
float Filter_BANDPASS = 0;
int i;
int pole = 0;
const int readPin = A1;
ADC *adc = new ADC();
int32_t adc_val = 0;
int startTimerValue = 0;
#define BUFFER_SIZE 16384 //16384 
uint16_t  BUFFER_ADC[BUFFER_SIZE];
double  sensorValueG;
double vReal[BUFFER_SIZE];
double vImag[BUFFER_SIZE];
float soucet = 0;
float g_rms = 0;
double frekvence;
float magnitude;
float suma;
float prumer;


unsigned int sampling_period_us;
unsigned long microseconds;
const double samplingFrequency = 20000;

void setup()
{    
 
  Serial.begin(115200); 
  pinMode(readPin, INPUT);    
  adc->adc0->setResolution(12); // set bits of resolution
  adc->adc0->setConversionSpeed(ADC_CONVERSION_SPEED::VERY_LOW_SPEED); // change the conversion speed
  adc->adc0->setSamplingSpeed(ADC_SAMPLING_SPEED::VERY_LOW_SPEED); // change the sampling speed   
  delayMicroseconds(25);  

   sampling_period_us = round(1000000*(1.0/samplingFrequency));
    
}


void loop()
{   
   // Serial.clear();


  //  F_HIGH_F();  
  //  SAMPLE();
  //  gRMS();
    rms_mm();
   
  }

void READ_ADC()
{
   adc_val = 0;
   pole=0;
   startTimerValue = timer.begin(TIMER, period); 
   timer.priority (0);  
   ADC_ISR();
}

void TIMER()   
 { 
  adc->adc0->startSingleRead(readPin);  
  adc->adc0->enableInterrupts(ADC_ISR);
 }
void ADC_ISR(void)
{   
   adc_val = adc->adc0->readSingle(); 
   if (pole < BUFFER_SIZE) BUFFER_ADC[pole++] = adc_val;
   asm volatile("dsb");         
 }


and


Code:
void F_LOW_F()
{
   alpha_LOW = 0.003; //  10hz
   alpha_HIGH =0.25; // 1000hz
}
void F_HIGH_F()
{
   alpha_LOW = 0.145; //  500hz
   alpha_HIGH =0.8; // 7500hz
}
void SAMPLE()
{ 
   period = 50; //20kHz 
   READ_ADC(); 
   delay(100);
   cas = 0;
   for(int i= 0;i < BUFFER_SIZE; i++)
   {  
    micros_0 = 0;  
    Filter_LOWPASS =  (BUFFER_ADC[i] *alpha_LOW) + (Filter_LOWPASS * (1.0 - alpha_LOW));
    Filter_HIGHPASS = (BUFFER_ADC[i] *alpha_HIGH) + (Filter_HIGHPASS * (1.0 - alpha_HIGH));
    Filter_BANDPASS =  Filter_HIGHPASS-Filter_LOWPASS;        
    sensorValueG = (((Filter_BANDPASS * 3.0) / 4095) - (3.0 / 2)) / 0.020; 
    sensorValueG =sensorValueG + 75; // offset  
    vReal[i] =    sensorValueG;    
    vImag[i] = period *i;
    while(micros_0 <=  period) 
    {      
    }     
    //Serial.println(vReal[i]);   
    //delay(2000);
    }  
 }  

void gRMS()
{ 
  for(int i= 0;i < BUFFER_SIZE; i++)
   {
    vReal[i] = pow(vReal[i],2);
    soucet = soucet + vReal[i];
   }
    soucet = soucet / BUFFER_SIZE;
    g_rms = sqrt(soucet);
    Serial.print("g RMS: ");
    Serial.println(g_rms);
    Serial.print("g Peak: ");
    Serial.println(g_rms * 1.414);
} 


void rms_mm()
{ 
  suma = 0;
  microseconds = micros();
  for(int i=0; i<BUFFER_SIZE; i++)
  {
      vReal[i] = analogRead(A1);
      vImag[i] = 0;
      while(micros() - microseconds < sampling_period_us)
      {
      }
      microseconds += sampling_period_us;
  }

  F_LOW_F();

 for(int i= 0;i < BUFFER_SIZE; i++)
   {  
    micros_0 = 0;  
    Filter_LOWPASS =  (vReal[i] *alpha_LOW) + (Filter_LOWPASS * (1.0 - alpha_LOW));
    Filter_HIGHPASS = (vReal[i] *alpha_HIGH) + (Filter_HIGHPASS * (1.0 - alpha_HIGH));
    Filter_BANDPASS =  Filter_HIGHPASS-Filter_LOWPASS;        
    sensorValueG = Filter_BANDPASS; 
    vReal[i] =    sensorValueG;    
    vImag[i] = period *i;
    while(micros_0 <=  period) 
    {      
    }     
   }  
    
  FFT.Windowing(vReal, BUFFER_SIZE, FFT_WIN_TYP_HANN, FFT_FORWARD);
  FFT.Compute(vReal, vImag, BUFFER_SIZE, FFT_FORWARD);
  FFT.ComplexToMagnitude(vReal, vImag, BUFFER_SIZE);

  for(int i= 0;i < BUFFER_SIZE / 2; i++)    //2 pro 10000Hz /20 pro 1000Hz
   {
    vReal[i] = pow(vReal[i], 2);
    frekvence = (i * 1.0 * 20000) / BUFFER_SIZE;    
    magnitude = vReal[i];  
    suma = suma + magnitude;
   }
   suma = sqrt(suma);
   Serial.println(suma);
  
}
 
Teensy 4.0

The Arduino FFT implementation is using double, which is OK for Teensy 4.x

Now we see also the whole code, which is not clear and most likely will not work as intended.
e.g. you call from SAMPLE()
Code:
void READ_ADC()
{
   adc_val = 0;
   pole=0;
   startTimerValue = timer.begin(TIMER, period); 
   timer.priority (0);  
   ADC_ISR();
}

You start a timer which triggers the ISR and then you call the ISR directly.
Later you process the data without waiting that buffer is filled

Suggest you get first a acquisition system ready
Typically start a periodic timer and let the ISR fill the bufffer.
in main loop wait for buffer to be filled the process buffer
In order to work on available buffer use ping-pong method where you fill two buffers and only work on buffer that is not filled.
make sure that you do not fill buffer while you work on it.

After that you
convert samples to volt or directly to m/s^2 (get conversion formula from ADXL documentation)
then you can integrate to get particle velocity (in time or spectrum)
but first things first
 
Suggest you get first a acquisition system ready
Typically start a periodic timer and let the ISR fill the bufffer.
in main loop wait for buffer to be filled the process buffer
In order to work on available buffer use ping-pong method where you fill two buffers and only work on buffer that is not filled.
make sure that you do not fill buffer while you work on it.

Could you help me with it? I have no idea how to do it...
 
Could you help me with it? I have no idea how to do it...
You have two options:
- start with examples of ADC libraries.
- use Audio library
in both cases try to read first a single channel,
if that works increase numbers of channels to three
if you get reasonable time series (you see response if you hit accelerometer etc) then thing about conversion to particle velocity.
 
I think Setup is OK:

Code:
Serial.begin(115200); 
  pinMode(readPin, INPUT);    
  adc->adc0->setResolution(12); // set bits of resolution
  adc->adc0->setConversionSpeed(ADC_CONVERSION_SPEED::VERY_LOW_SPEED); // change the conversion speed
  adc->adc0->setSamplingSpeed(ADC_SAMPLING_SPEED::VERY_LOW_SPEED); // change the sampling speed   
  delayMicroseconds(25);  
  sampling_period_us = round(1000000*(1.0/samplingFrequency));

Here is mistake? Where?

Code:
void READ_ADC()
{
   adc_val = 0;
   pole=0;
   startTimerValue = timer.begin(TIMER, period); 
   timer.priority (0);  
   ADC_ISR();
}

void TIMER()   
 { 
  adc->adc0->startSingleRead(readPin);  
  adc->adc0->enableInterrupts(ADC_ISR);
 }
void ADC_ISR(void)
{   
   adc_val = adc->adc0->readSingle(); 
   if (pole < BUFFER_SIZE) BUFFER_ADC[pole++] = adc_val;
   asm volatile("dsb");         
 }

This works fine I think:

Code:
void SAMPLE()
{ 
   period = 50; //20kHz 
   READ_ADC(); 
   delay(1000);
   cas = 0;
   for(int i= 0;i < BUFFER_SIZE; i++)
   {  
    micros_0 = 0;  
    Filter_LOWPASS =  (BUFFER_ADC[i] *alpha_LOW) + (Filter_LOWPASS * (1.0 - alpha_LOW));
    Filter_HIGHPASS = (BUFFER_ADC[i] *alpha_HIGH) + (Filter_HIGHPASS * (1.0 - alpha_HIGH));
    Filter_BANDPASS =  Filter_HIGHPASS-Filter_LOWPASS;        
    sensorValueG = (((Filter_BANDPASS * 3.0) / 4095) - (3.0 / 2)) / 0.020; 
    sensorValueG =sensorValueG + 75; //75 = calculated average noise of ADXL when no movement
    vReal[i] =    sensorValueG;    
    vImag[i] = period *i;
    while(micros_0 <=  period) 
    {      
    }     
    }  
 }


Calculation of gRMS is OK I mean - no problem to calculate it from time signal / wave - checked

Code:
void gRMS()
{ 
  for(int i= 0;i < BUFFER_SIZE; i++)
   {
    vReal[i] = pow(vReal[i],2);
    soucet = soucet + vReal[i];
   }
    soucet = soucet / BUFFER_SIZE;
    g_rms = sqrt(soucet);
    Serial.print("g RMS: ");
    Serial.println(g_rms);
}


But for mm/s I have to use FFT - I mean it is not possible to calculate mm/s from time signal / wave.
I need FFT with correct amplitude in mV or g and NOT just magnitude.
If I will have for example 8196 values values from FFT after Hann windowing, then I am able assign to every frequency right amplitude and then is easy to calculate velocity.
 
Back
Top