3

I have been stuck for a long time on a piece of Arduino code that performs an FFT on a Photoplethysmogram signal and decides signal quality after that.

#include <math.h>
#include "arduinoFFT.h"
#include <EEPROM.h>

//N=10 pt moving average
//For bidmc data of first 10 patients, MA = 0.4613849
double baseline=0.4613849;
float lm=1.0; //amplitude threshold
float ref_voltage=1.0; //saturation reference voltage
float cutoff=4.0;
float fs=0.02;
float zero_amp=0.015;

float pi=3.14159265;  //Universal Constants
float e=2.71828;
const int flen=61;

#define SAMPLES 128           //Must be a power of 2
#define SAMPLING_FREQUENCY 50 //Hz, must be less than 10000 due to ADC

arduinoFFT FFT = arduinoFFT(); 
unsigned int sampling_period_us;
unsigned long microseconds;

double vReal[SAMPLES];
double vImag[SAMPLES];

void setup() {
  Serial.begin(9600);
  sampling_period_us = round(1000000*(1.0/SAMPLING_FREQUENCY));
}

void loop() {
  //float arr[]={0.43597,0.43206,0.42815,0.42424,0.42131,0.41838,0.4174,0.41642,0.41544,0.41447,0.41251,0.41153,0.40958,0.40762,0.40665,0.40665,0.40665,0.40762,0.4086,0.4086,0.41056,0.41153,0.41153,0.41056,0.4086,0.40665,0.40469,0.40469,0.40567,0.4086,0.41251,0.42033,0.43206,0.45259,0.4741,0.4956,0.51808,0.5347,0.55034,0.56403,0.57576,0.58456,0.59238,0.59629,0.59922,0.60117,0.60215,0.60313,0.60313,0.60117,0.59824,0.59335,0.58553,0.57771,0.57185,0.56207,0.54741,0.53372,0.52102,0.50929,0.50244,0.49756,0.4956,0.49462,0.49365,0.49365,0.49365,0.49169,0.48876,0.4868,0.48387,0.48094,0.47801,0.47507,0.47116,0.46823,0.46334,0.45846,0.45357,0.44868,0.44379,0.43891,0.435,0.43109,0.4262,0.42326,0.42033,0.41838,0.41642,0.41447,0.41447,0.41349,0.41251,0.41251,0.41251,0.41251,0.41251,0.41349,0.41544,0.4174,0.41935,0.42033,0.42131,0.42131,0.42033,0.41935,0.41642,0.41447,0.41349,0.41349,0.41642,0.42131,0.42815,0.43695,0.45064,0.46823,0.4868,0.50538,0.52395,0.54154,0.55816,0.57185,0.58162,0.58749,0.59042,0.59042,0.5914,0.59042,0.58944,0.58651,0.5826,0.57771,0.57087,0.56305,0.5523,0.54057,0.52884,0.51711,0.50733,0.49756,0.48778,0.47996,0.4741,0.47019,0.46823,0.46725,0.46628,0.46432,0.46041,0.4565}; 
  float arr[]={0.30303,0.29814,0.29326,0.28837,0.28446,0.28152,0.27859,0.27468,0.26979,0.26393,0.25806,0.25611,0.25904,0.27175,0.29521,0.32942,0.37341,0.42522,0.47116,0.51417,0.55914,0.6002,0.63539,0.66373,0.68622,0.70381,0.71652,0.72532,0.73118,0.73314,0.73314,0.73021,0.7263,0.72043,0.71359,0.70479,0.69501,0.68524,0.67351,0.6608,0.64614,0.6305,0.61388,0.59726,0.58065,0.56207,0.54252,0.52395,0.50538,0.48778,0.47312,0.46041,0.44868,0.43988,0.43304,0.42815,0.42522,0.42229,0.41935,0.4174,0.41544,0.41349,0.41056,0.40762,0.40469,0.40176,0.39785,0.39394,0.39003,0.38612,0.38221,0.37634,0.37146,0.36559,0.35973,0.35386,0.348,0.34213,0.33627,0.3304,0.32551,0.32258,0.31769,0.31281,0.30792,0.30303,0.29814,0.29423,0.28935,0.28446,0.27957,0.2737,0.26784,0.26393,0.261,0.26197,0.26882,0.28446,0.31183,0.34897,0.39589,0.44868,0.50147,0.55132,0.59629,0.6305,0.65591,0.67937,0.69795,0.71261,0.72434,0.73216,0.73607,0.73803,0.73803,0.73509,0.73021,0.72434,0.7175,0.70968,0.6999,0.69013,0.6784,0.66569,0.65103,0.63441,0.61681,0.59824,0.57771,0.55718,0.53959,0.52297,0.50635,0.48974,0.47507,0.46334,0.45259,0.44379,0.43695,0.43206,0.42815,0.42522,0.42229,0.42033,0.41838,0.41642,0.41447,0.41153,0.4086,0.40469,0.40078,0.39687,0.39198,0.3871,0.38221,0.3783,0.37439,0.3695,0.36364,0.35875,0.35288,0.348,0.34213,0.33724,0.33236,0.32747,0.32258,0.31769,0.31378,0.31085,0.30792,0.30401,0.3001,0.29619,0.29032,0.28446,0.27859,0.2737,0.2737,0.28055,0.29423,0.31672,0.35191,0.39785,0.45161,0.50733,0.56109,0.60899,0.65005,0.68328,0.7087,0.72727,0.74096,0.74878,0.75367,0.75562,0.7566,0.75464,0.75171,0.74585};
  //float arr[]={0.30303,0.28446,0.26979,0.25904,0.37341,0.55914,0.68622,0.73118};
  int len=sizeof(arr)/sizeof(arr[0]);
  
  /*for(int i=0;i<sizeof(arr)/sizeof(arr[0]);i++)  //ppg data Array display
  {
    Serial.println(arr[i]);  
    delay(50);
  }*/
  
  //SAMPLES=len;
  for(int i=0; i<SAMPLES; i++)
    {
        microseconds = micros();    //Overflows after around 70 minutes!
        vReal[i] = arr[SAMPLES]; 
        vImag[i] = 0;
        
        while(micros() < (microseconds + sampling_period_us)){
        }
    }
    
  /*FFT*/
  FFT.Windowing(vReal, SAMPLES, FFT_WIN_TYP_HAMMING, FFT_FORWARD);
  FFT.Compute(vReal, vImag, SAMPLES, FFT_FORWARD);
  FFT.ComplexToMagnitude(vReal, vImag, SAMPLES);
  double peak = FFT.MajorPeak(vReal, SAMPLES, SAMPLING_FREQUENCY);
   
  /*PRINT RESULTS*/
  //Serial.println(peak);     //Print out what frequency is the most dominant.
  
  for(int i=0; i<(SAMPLES/2); i++)
  {
       /*View all these three lines in serial terminal to see which frequencies has which amplitudes*/ 
       Serial.print((i * 1.0 * SAMPLING_FREQUENCY) / SAMPLES, 2);
       Serial.print(" ");
       Serial.println(vReal[i], 1);  //View only this line in serial plotter to visualize the bins
       delay(50);
  }
   
  //delay(1000);  //Repeat the process every second OR:
  //while(1);       //Run code once
  Serial.println("Hello");
  //int freq_N=61;
  for (int i = 0 ; i < EEPROM.length() ; i++) {
    EEPROM.write(i, 0);
  }
  
  /*double D[3][freq_N]; //1st row for frequencies, 2nd row for magnitude, 3rd row for phase
  for(int i=0;i<freq_N;i++)
  {
    D[0][i]=-(pi/2)+i*0.1;
    double *d=fourier_tran(arr,len,D[0][i]);
    D[1][i]=*(d+0);
    D[2][i]=*(d+1); 
    Serial.println(i);
    delay(300);
  }*/
  
  int z_crsg=zero_cross(arr,len);
  //Serial.println(z_crsg);
  //delay(1000);
  float diff_arr[len-1];   //Differenced PPG signal
  for(int i=1;i<=len-1;i++)
  {
    diff_arr[i-1]=arr[i]-arr[i-1];
  }
  
  int rise=rise_count(arr,len);
  int fall=fall_count(arr,len);
  float fr=(float)fall/rise;      //fall/rise count ratio 
  uint16_t m=max_amplitude(arr,len);   //maximum amplitude detection
  int sc[]={0,0,0};  //sensor connectivity vector
  if(fr>1.65 && fr<2.35)
  { 
    sc[0]=1;
  }
  if(m<lm)
  {
    sc[1]=1;
  }
  if(m>zero_amp)
  {
    sc[2]=1;;
  }
  int sd=sig_saturation(arr,len);
  //Serial.println(sd);
  //Serial.print(sc[0]);
  //Serial.print(sc[1]);
  //Serial.print(sc[2]);
 
  if(sum(sc,3)>=2)
  {
    Serial.println("Sensor well connected");
  }
  else
  {
    Serial.println("Sensor not connected properly");
  }
  delay(500);
  if(base_drift(arr,len)>=0.4612 && base_drift(arr,len)<=0.4614)
  {
    Serial.println("Motion detected");
  }
  else
  {
    Serial.println("No motion artefacts problems");
  }
  delay(500);
  if(sd==0)
  {
    Serial.println("Saturated signal");
  }
  else
  {
    Serial.println("Non-saturated signal");
  }
  delay(500); 
  
  int d_z_crsg=zero_cross(diff_arr,len-1);
  
}

//Sensor connectivity detection 
int fall_count(float arr[],int l)
{
  int c=0;
  for(int i=0;i<l-1;i++)
  {
    if(arr[i]>arr[i+1])
    {
      c=c+1;
    }
  }
  return c;
}
int rise_count(float arr[],int l)
{
  int c=0;
  for(int i=0;i<l-1;i++)
  {
    if(arr[i]<arr[i+1])
    {
      c=c+1;
    }
  }
  return c;
}
int max_amplitude(float arr[],int l)
{
  float h=arr[0];
  for(int i=1;i<l;i=i+1)
  {
    if(arr[i]>h)
    {
      h=arr[i];
    }
  }
  if(h>lm)
  {
    return 0;
  }
  return 1;
  
}

//Noise segment detection
int zero_cross(float arr[],int l)
{
  int c=0;
  for(int i=0;i<l-1;i++)
  {
    if((arr[i]==0.0) && arr[i+1]!=0.0)
    {
      c=c+1; 
    }
  }
  return c;
}
int local_maxima(float arr[], int l)
{
  int frame_len=15;
  int a=0;
  int frames=l/frame_len;
  int maxima[frames];
  for(int k=0;k<frames;k++)
  {
    uint16_t h=0;
    for(int i=a;i<a+frame_len;i++)
    {
      if(arr[i]>h)
      {
        h=arr[i];
      }
    }
    maxima[k]=h;
    a=a+frame_len;
  }
  for(int i=0;i<frames-1;i++)
  {
    if(maxima[i]>lm && maxima[i+1]>lm)
    {
      return 1;
    }
    if(i<=frames-4)
    {
      if(maxima[i]<lm && maxima[i+1]<lm && maxima[i+2]<lm && maxima[i+3]<lm)
      {
        return 1;
      }
    }
    if(i<=frames-3)
    {
      if(maxima[i]<lm && maxima[i+2]<lm)
      {
        return 1;
      }
    }
  }
  return 0;
}


//Saturation detection
int sig_saturation(float arr[],int l)
{
  int c=0;
  //uint16_t lev=200000;
  for(int i=0;i<l-1;i++)
  {
    if(arr[i]<10 && arr[i+1]<10)
    {
      c++;
    }
    
    if(arr[i]>lm && arr[i+1]>lm)
    {
      c++;
    }
  }
  if(c>=20)
  {
    return 0;
  }
  return 1;
}

//motion detection module
double base_drift(float arr[],int l)
{
  int N=10; //N pt moving average
  int frame=l-N+1; //number of frames
  float s=0;
  double sum=0;
  for(int i=0;i<frame;i++)
  {
   s=0;
   for(int j=i;j<N+i;i++)
   {
      s=s+arr[j];
   }
   sum=sum+(s/N);
  }
  double avg=sum/frame;
  return avg;
}

int sum(int n[],int l)
{
  float s=0;
  for(int i=0;i<l;i++)
  {
    s=s+n[i];
  }
  return s;
}

float butterworth(float f, float fc,int n)
{
  float h;
  h=1/sqrt(1+pow((fc/f),2*n)); 
  return h;
}

double *fourier_tran(float arr[],int n,float omega)
{
   double Xr=0.00;
   double Xi=0.00;
   double DTFT[2];
   for(int j=0;j<n;j++)
   {
     Xr=Xr+arr[j]*cos(omega*j);
     Xi=Xi+arr[j]*sin(omega*j);
   }
   float X=sqrt(pow(Xr,2)+pow(Xi,2));  //Fourier transform magnitude
   float ph=-atan(Xi/Xr);  //Fourier transform phase
   
   DTFT[0]=X;
   DTFT[1]=ph;
   double *p=DTFT;
   Serial.println(*p);
   delay(100);
   Serial.println(*(p+1));
   delay(100);
   return p;
}
float *inv_fourier(float DTFT[][flen],int l,int flen,int n) //n is the length of the time array 
{
  float arr[n];
  for(int i=0;i<n;i++)
  {
    float X=0.000;
    for(int j=0;j<l;j++)
    {
      X=X+DTFT[1][j]*cos(DTFT[2][j]+DTFT[0][j]*n);
    }
    arr[i]=1/(2*pi) * X;
  }
  float *p=arr;
  return p;
}

When I upload the code, Right after performing the FFT, the code does not display anything more, implying it ceases to execute and shows an error message - Not enough memory; see http://www.arduino.cc/en/Guide/Troubleshooting#size for tips on reducing your footprint. data section exceeds available space in board

I cannot remove either section of the code, so what do I need to change with this?

9
  • This error is normally shown while attempting to upload a sketch, not during execution. Are you able to successfully load the sketch and the error only happens after the sensor has been read? – PMF Feb 1 at 6:28
  • It happens during uploading(sorry I didn;t mention, edited). – Divyayan Dey Feb 1 at 7:02
  • What program and data sizes does the compiler show right after your sketch has been compiled? – StarCat Feb 1 at 7:13
  • Your code as posted doesn't compile. Please post the entire sketch. We cannot debug snippets. – Nick Gammon Feb 1 at 7:21
  • 1
    With all those floats, have you considered that an 8-bit micro with no FPU is not the optimal MCU for this. – r_ahlskog Feb 1 at 11:08
2

Probably you use too much SRAM memory:

float arr[] seems to take 4 bytes (for float) * 128 = 512 bytes.

Then you have vReal and `vImage which takes 4 (double, float/doubles are the same in Arduino) * 128 * 2 (vReal/vImage) = 1024 bytes.

double vReal[SAMPLES];
double vImag[SAMPLES];

This already is 1536 bytes. The Arduino only has 2048 bytes so high likely you exceed it (including arguments passed and other globals).

arr seems to be a readonly array, so check into the use of F/PROGMEM to store them in Flash instead of SRAM (see StarCat's comment). This might get you well under the 2 KB Arduino Uno memory. Other possibilities are using a Mega which has 8 KB or use external SRAM memory but this decreases performance and takes more programming effort.

See also my anser at Optimizing code to use less Flash Memory and SRAM

1
  • 1
    See my comment: string constants (such as those used in Serial.print("");) also take up valuable RAM. You can use the F("...") macro to reduce their footprint. – StarCat Feb 1 at 10:59
0

Consider using one of the Teensy family. The latest Teensy 4.0 has an onboard FPU and lots of other goodies, and runs at 600MHz. Why drive a VW when you can jump into a Ferrari for very little money?

I've been using the Teensy 3.2/3.5 modules for years, and have found them to be very stable and easy to use.

2
  • Wouldn't comparing the memory available be the most helpful selling point for this question? – chicks Feb 3 at 17:28
  • Too lazy to click on the link? – user3765883 Feb 3 at 21:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.