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)+i0.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(omegaj);
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/(2pi) 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?
F("")macro on theSerial.println("");strings (storing them in Flash instead of copying them to RAM), I've been able to reduce the Global Variables to 2038 bytes, but that still leaves too little for running the code. – StarCat Feb 01 '21 at 10:56