ApproxFFT: Fastest FFT Function for Arduino

This function performs FFT with very good speed while maintaining accuracy.

Dec 23, 2020

•

25026 views

•

1 respects

fft

frequency

Components and supplies

Arduino Mega 2560

Apps and platforms

Arduino IDE

Project description

Code

ApproxFFT

arduino

This code performs FFT with good speed

1//---------------------------------lookup data------------------------------------//
2byte  isin_data[128]=
3{0,  1,   3,   4,   5,   6,   8,   9,   10,  11,  13,  14,  15,  17,  18,  19,  20, 
422,  23,  24,  26,  27,  28,  29,  31,  32,  33,  35,  36,  37,  39,  40,  41,  42, 
544,  45,  46,  48,  49,  50,  52,  53,  54,  56,  57,  59,  60,  61,  63,  64,  65, 
667,  68,  70,  71,  72,  74,  75,  77,  78,  80,  81,  82,  84,  85,  87,  88,  90, 
791,  93,  94,  96,  97,  99,  100, 102,  104, 105, 107, 108, 110, 112, 113, 115, 117, 
8118, 120, 122, 124, 125, 127, 129,  131, 133, 134, 136, 138, 140, 142, 144, 146, 148, 
9150, 152, 155, 157, 159, 161,  164, 166, 169, 171, 174, 176, 179, 182, 185, 188, 191, 
10195, 198, 202, 206, 210,  215, 221, 227, 236};
11unsigned int Pow2[14]={1,2,4,8,16,32,64,128,256,512,1024,2048,4096};
12byte  RSSdata[20]={7,6,6,5,5,5,4,4,4,4,3,3,3,3,3,3,3,2,2,2};
13//---------------------------------------------------------------------------------//
14
15
16//int  data[256]={};
17
18
19void setup() 
20        {
21     Serial.begin(250000);  
22        }
23
24
25void loop() {
26
27//float f=Approx_FFT(data,512,100);
28//Serial.println(f);
29            }
30
31
32//-----------------------------FFT Function----------------------------------------------//
33/*
34Code  to perform High speed and Accurate FFT on arduino,
35setup:
36
371. in[]     :  Data array, 
382. N        : Number of sample (recommended sample size 2,4,8,16,32,64,128,256,512...)
393.  Frequency: sampling frequency required as input (Hz)
40
41It will by default return  frequency with max aplitude,
42if you need complex output or magnitudes uncomment  required sections
43
44If sample size is not in power of 2 it will be clipped  to lower side of number. 
45i.e, for 150 number of samples, code will consider  first 128 sample, remaining sample  will be omitted.
46For Arduino nano, FFT of  more than 256 sample not possible due to mamory limitation 
47Code by ABHILASH
48Contact:  abhilashpatel121@gmail.com
49Documentation & details: https://www.instructables.com/member/abhilash_patel/instructables/
50*/
51
52float  Approx_FFT(int in[],int N,float Frequency)
53{ 
54int a,c1,f,o,x,data_max,data_min=0;
55long  data_avg,data_mag,temp11;         
56byte scale,check=0;
57
58data_max=0;
59data_avg=0;
60data_min=0;
61
62      for(int i=0;i<12;i++)                 //calculating the levels
63         {  if(Pow2[i]<=N){o=i;} }
64     a=Pow2[o];  
65int out_r[a];   //real part of transform
66int  out_im[a];  //imaginory part of transform
67
68      for(int i=0;i<a;i++)                //getting  min max and average for scalling
69          { out_r[i]=0; out_r[i]=0;
70            data_avg=data_avg+in[i];
71            if(in[i]>data_max){data_max=in[i];}
72            if(in[i]<data_min){data_min=in[i];}
73          }
74
75data_avg=data_avg>>o;
76scale=0;
77data_mag=data_max-data_min;
78temp11=data_mag;
79
80for(int  i;i<128;i++)              //scalling data  from +512 to -512
81
82      if(data_mag>1024)
83          {while(temp11>1024)
84                {temp11=temp11>>1;
85                scale=scale+1;
86                }   
87          }
88 
89      if(data_mag<1024)
90          {while(temp11<1024)
91                {temp11=temp11<<1;
92                scale=scale+1;
93                }
94          }
95
96
97      if(data_mag>1024)
98          {
99              for(int  i=0;i<a;i++)
100                    { in[i]=in[i]-data_avg;
101                      in[i]=in[i]>>scale;
102                    }
103                    scale=128-scale;
104          }
105
106      if(data_mag<1024)
107          { scale=scale-1;
108            for(int i=0;i<a;i++)
109                    {
110                      in[i]=in[i]-data_avg;
111                      in[i]=in[i]<<scale;
112                    }
113
114                    scale=128+scale;
115          }
116
117       
118x=0;  
119      for(int b=0;b<o;b++)                     // bit reversal  order stored in im_out array
120         {
121          c1=Pow2[b];
122          f=Pow2[o]/(c1+c1);
123                for(int j=0;j<c1;j++)
124                    { 
125                     x=x+1;
126                     out_im[x]=out_im[j]+f;
127                    }
128         }
129
130  
131      for(int i=0;i<a;i++)            // update input array as per bit reverse  order
132         {
133          out_r[i]=in[out_im[i]]; 
134          out_im[i]=0;
135         }
136
137
138int i10,i11,n1,tr,ti;
139float e;
140int c,s,temp4;
141    for(int  i=0;i<o;i++)                                    //fft
142    {
143     i10=Pow2[i];              // overall values of sine/cosine  
144     i11=Pow2[o]/Pow2[i+1];    // loop with similar sine cosine
145     e=1024/Pow2[i+1];  //1024 is equivalent  to 360 deg
146     e=0-e;
147     n1=0;
148
149          for(int j=0;j<i10;j++)
150          {
151            c=e*j;    //c is angle as where 1024 unit is 360 deg
152  while(c<0){c=c+1024;}
153  while(c>1024){c=c-1024;}
154
155          n1=j;
156          
157          for(int k=0;k<i11;k++)
158                 {
159                   temp4=i10+n1;
160       if(c==0)   {tr=out_r[temp4];
161                   ti=out_im[temp4];}
162  else if(c==256) {tr= -out_im[temp4];
163                   ti=out_r[temp4];}
164  else if(c==512) {tr=-out_r[temp4];
165                  ti=-out_im[temp4];}
166  else if(c==768) {tr=out_im[temp4];
167                   ti=-out_r[temp4];}
168  else if(c==1024){tr=out_r[temp4];
169                   ti=out_im[temp4];}
170  else{
171    tr=fast_cosine(out_r[temp4],c)-fast_sine(out_im[temp4],c);            //the  fast sine/cosine function gives direct (approx) output for A*sinx
172    ti=fast_sine(out_r[temp4],c)+fast_cosine(out_im[temp4],c);            
173      }
174          
175                 out_r[n1+i10]=out_r[n1]-tr;
176                 out_r[n1]=out_r[n1]+tr;
177                 if(out_r[n1]>15000  || out_r[n1]<-15000){check=1;}   //check for int size, it can handle only +31000  to -31000,
178          
179                 out_im[n1+i10]=out_im[n1]-ti;
180                 out_im[n1]=out_im[n1]+ti;
181                 if(out_im[n1]>15000 || out_im[n1]<-15000){check=1;}          
182          
183                 n1=n1+i10+i10;
184                  }       
185             }
186
187    if(check==1){                                             //  scalling the matrics if value higher than 15000 to prevent varible from overflowing
188                for(int i=0;i<a;i++)
189                    {
190                     out_r[i]=out_r[i]>>1;           
191                     out_im[i]=out_im[i]>>1; 
192                    }
193                     check=0; 
194                     scale=scale-1;                 //  tracking overall scalling of input data
195                }           
196
197     }
198
199
200if(scale>128)
201    {scale=scale-128;
202     for(int i=0;i<a;i++)
203      {out_r[i]=out_r[i]>>scale;
204       out_im[i]=out_im[i]>>scale;
205      }
206      scale=0;
207    }                                                   // revers  all scalling we done till here,
208else{scale=128-scale;}                             //  in case of nnumber getting higher than 32000, we will represent in as multiple of  2^scale
209
210/*
211for(int i=0;i<a;i++)
212{
213Serial.print(out_r[i]);Serial.print("\	");                    // un comment to print RAW o/p    
214Serial.print(out_im[i]);  
215Serial.print("i");Serial.print("\	"); 
216Serial.print("*2^");Serial.println(scale);  
217}
218*/
219
220//---> here onward out_r contains amplitude and our_in conntains  frequency (Hz)
221int fout,fm,fstp;
222float fstep;
223fstep=Frequency/N;
224fstp=fstep;
225fout=0;fm=0;
226
227    for(int i=1;i<Pow2[o-1];i++)               // getting amplitude from compex  number
228        { 
229              out_r[i]=fastRSS(out_r[i],out_im[i]);
230   // Approx RSS function used to calculated magnitude quickly
231        
232out_im[i]=out_im[i-1]+fstp;
233if  (fout<out_r[i]){fm=i; fout=out_r[i];}
234         /*
235         // un comment to  print Amplitudes (1st value (offset) is not printed)
236         Serial.print(out_r[i]);  Serial.print("\	"); 
237         Serial.print("*2^");Serial.println(scale);  
238         */ 
239        }
240
241
242float fa,fb,fc;
243fa=out_r[fm-1];
244fb=out_r[fm];  
245fc=out_r[fm+1];
246fstep=(fa*(fm-1)+fb*fm+fc*(fm+1))/(fa+fb+fc);
247
248return(fstep*Frequency/N);
249}
250
251//---------------------------------fast  sine/cosine---------------------------------------//
252
253int fast_sine(int Amp,  int th)
254{
255int temp3,m1,m2;
256byte temp1,temp2, test,quad,accuracy;
257accuracy=5;    // set it value from 1 to 7, where 7 being most accurate but slowest
258               //  accuracy value of 5 recommended for typical applicaiton
259while(th>1024){th=th-1024;}   // here 1024 = 2*pi or 360 deg
260while(th<0){th=th+1024;}
261quad=th>>8;
262
263       if(quad==1){th= 512-th;}
264  else if(quad==2){th= th-512;}
265  else if(quad==3){th=  1024-th;}
266
267temp1= 0;
268temp2= 128;     //2 multiple
269m1=0;
270m2=Amp;
271
272    temp3=(m1+m2)>>1;
273    Amp=temp3;
274      for(int i=0;i<accuracy;i++)
275        { test=(temp1+temp2)>>1;
276          temp3=temp3>>1; 
277          if(th>isin_data[test]){temp1=test;  Amp=Amp+temp3; m1=Amp;}
278          else if(th<isin_data[test]){temp2=test; Amp=Amp-temp3;  m2=Amp;}
279        }
280
281         if(quad==2){Amp= 0-Amp;}
282    else if(quad==3){Amp=  0-Amp;}
283return(Amp);
284}
285
286int fast_cosine(int Amp, int th)
287  {  
288  th=256-th;  //cos th = sin (90-th) formula
289  return(fast_sine(Amp,th)); 
290  }
291
292//--------------------------------------------------------------------------------//
293
294
295//--------------------------------Fast  RSS----------------------------------------//
296int fastRSS(int a, int b)
297{  if(a==0 && b==0){return(0);}
298  int min,max,temp1,temp2;
299  byte clevel;
300  if(a<0){a=-a;}
301  if(b<0){b=-b;}
302clevel=0;
303if(a>b){max=a;min=b;} else{max=b;min=a;}
304
305  if(max>(min+min+min))
306    {return max;}
307  else
308    {
309     temp1=min>>3;  if(temp1==0){temp1=1;}
310     temp2=min;
311     while(temp2<max){temp2=temp2+temp1;clevel=clevel+1;}
312     temp2=RSSdata[clevel];temp1=temp1>>1;  
313     for(int i=0;i<temp2;i++){max=max+temp1;}
314     return(max);
315    }
316}
317//--------------------------------------------------------------------------------//

//---------------------------------lookup data------------------------------------//
byte  isin_data[128]=
{0,  1,   3,   4,   5,   6,   8,   9,   10,  11,  13,  14,  15,  17,  18,  19,  20, 
22,  23,  24,  26,  27,  28,  29,  31,  32,  33,  35,  36,  37,  39,  40,  41,  42, 
44,  45,  46,  48,  49,  50,  52,  53,  54,  56,  57,  59,  60,  61,  63,  64,  65, 
67,  68,  70,  71,  72,  74,  75,  77,  78,  80,  81,  82,  84,  85,  87,  88,  90, 
91,  93,  94,  96,  97,  99,  100, 102,  104, 105, 107, 108, 110, 112, 113, 115, 117, 
118, 120, 122, 124, 125, 127, 129,  131, 133, 134, 136, 138, 140, 142, 144, 146, 148, 
150, 152, 155, 157, 159, 161,  164, 166, 169, 171, 174, 176, 179, 182, 185, 188, 191, 
195, 198, 202, 206, 210,  215, 221, 227, 236};
unsigned int Pow2[14]={1,2,4,8,16,32,64,128,256,512,1024,2048,4096};
byte  RSSdata[20]={7,6,6,5,5,5,4,4,4,4,3,3,3,3,3,3,3,2,2,2};
//---------------------------------------------------------------------------------//


//int  data[256]={};


void setup() 
        {
     Serial.begin(250000);  
        }


void loop() {

//float f=Approx_FFT(data,512,100);
//Serial.println(f);
            }


//-----------------------------FFT Function----------------------------------------------//
/*
Code  to perform High speed and Accurate FFT on arduino,
setup:

1. in[]     :  Data array, 
2. N        : Number of sample (recommended sample size 2,4,8,16,32,64,128,256,512...)
3.  Frequency: sampling frequency required as input (Hz)

It will by default return  frequency with max aplitude,
if you need complex output or magnitudes uncomment  required sections

If sample size is not in power of 2 it will be clipped  to lower side of number. 
i.e, for 150 number of samples, code will consider  first 128 sample, remaining sample  will be omitted.
For Arduino nano, FFT of  more than 256 sample not possible due to mamory limitation 
Code by ABHILASH
Contact:  abhilashpatel121@gmail.com
Documentation & details: https://www.instructables.com/member/abhilash_patel/instructables/
*/

float  Approx_FFT(int in[],int N,float Frequency)
{ 
int a,c1,f,o,x,data_max,data_min=0;
long  data_avg,data_mag,temp11;         
byte scale,check=0;

data_max=0;
data_avg=0;
data_min=0;

      for(int i=0;i<12;i++)                 //calculating the levels
         {  if(Pow2[i]<=N){o=i;} }
     a=Pow2[o];  
int out_r[a];   //real part of transform
int  out_im[a];  //imaginory part of transform

      for(int i=0;i<a;i++)                //getting  min max and average for scalling
          { out_r[i]=0; out_r[i]=0;
            data_avg=data_avg+in[i];
            if(in[i]>data_max){data_max=in[i];}
            if(in[i]<data_min){data_min=in[i];}
          }

data_avg=data_avg>>o;
scale=0;
data_mag=data_max-data_min;
temp11=data_mag;

for(int  i;i<128;i++)              //scalling data  from +512 to -512

      if(data_mag>1024)
          {while(temp11>1024)
                {temp11=temp11>>1;
                scale=scale+1;
                }   
          }
 
      if(data_mag<1024)
          {while(temp11<1024)
                {temp11=temp11<<1;
                scale=scale+1;
                }
          }


      if(data_mag>1024)
          {
              for(int  i=0;i<a;i++)
                    { in[i]=in[i]-data_avg;
                      in[i]=in[i]>>scale;
                    }
                    scale=128-scale;
          }

      if(data_mag<1024)
          { scale=scale-1;
            for(int i=0;i<a;i++)
                    {
                      in[i]=in[i]-data_avg;
                      in[i]=in[i]<<scale;
                    }

                    scale=128+scale;
          }

       
x=0;  
      for(int b=0;b<o;b++)                     // bit reversal  order stored in im_out array
         {
          c1=Pow2[b];
          f=Pow2[o]/(c1+c1);
                for(int j=0;j<c1;j++)
                    { 
                     x=x+1;
                     out_im[x]=out_im[j]+f;
                    }
         }

  
      for(int i=0;i<a;i++)            // update input array as per bit reverse  order
         {
          out_r[i]=in[out_im[i]]; 
          out_im[i]=0;
         }


int i10,i11,n1,tr,ti;
float e;
int c,s,temp4;
    for(int  i=0;i<o;i++)                                    //fft
    {
     i10=Pow2[i];              // overall values of sine/cosine  
     i11=Pow2[o]/Pow2[i+1];    // loop with similar sine cosine
     e=1024/Pow2[i+1];  //1024 is equivalent  to 360 deg
     e=0-e;
     n1=0;

          for(int j=0;j<i10;j++)
          {
            c=e*j;    //c is angle as where 1024 unit is 360 deg
  while(c<0){c=c+1024;}
  while(c>1024){c=c-1024;}

          n1=j;
          
          for(int k=0;k<i11;k++)
                 {
                   temp4=i10+n1;
       if(c==0)   {tr=out_r[temp4];
                   ti=out_im[temp4];}
  else if(c==256) {tr= -out_im[temp4];
                   ti=out_r[temp4];}
  else if(c==512) {tr=-out_r[temp4];
                  ti=-out_im[temp4];}
  else if(c==768) {tr=out_im[temp4];
                   ti=-out_r[temp4];}
  else if(c==1024){tr=out_r[temp4];
                   ti=out_im[temp4];}
  else{
    tr=fast_cosine(out_r[temp4],c)-fast_sine(out_im[temp4],c);            //the  fast sine/cosine function gives direct (approx) output for A*sinx
    ti=fast_sine(out_r[temp4],c)+fast_cosine(out_im[temp4],c);            
      }
          
                 out_r[n1+i10]=out_r[n1]-tr;
                 out_r[n1]=out_r[n1]+tr;
                 if(out_r[n1]>15000  || out_r[n1]<-15000){check=1;}   //check for int size, it can handle only +31000  to -31000,
          
                 out_im[n1+i10]=out_im[n1]-ti;
                 out_im[n1]=out_im[n1]+ti;
                 if(out_im[n1]>15000 || out_im[n1]<-15000){check=1;}          
          
                 n1=n1+i10+i10;
                  }       
             }

    if(check==1){                                             //  scalling the matrics if value higher than 15000 to prevent varible from overflowing
                for(int i=0;i<a;i++)
                    {
                     out_r[i]=out_r[i]>>1;           
                     out_im[i]=out_im[i]>>1; 
                    }
                     check=0; 
                     scale=scale-1;                 //  tracking overall scalling of input data
                }           

     }


if(scale>128)
    {scale=scale-128;
     for(int i=0;i<a;i++)
      {out_r[i]=out_r[i]>>scale;
       out_im[i]=out_im[i]>>scale;
      }
      scale=0;
    }                                                   // revers  all scalling we done till here,
else{scale=128-scale;}                             //  in case of nnumber getting higher than 32000, we will represent in as multiple of  2^scale

/*
for(int i=0;i<a;i++)
{
Serial.print(out_r[i]);Serial.print("\	");                    // un comment to print RAW o/p    
Serial.print(out_im[i]);  
Serial.print("i");Serial.print("\	"); 
Serial.print("*2^");Serial.println(scale);  
}
*/

//---> here onward out_r contains amplitude and our_in conntains  frequency (Hz)
int fout,fm,fstp;
float fstep;
fstep=Frequency/N;
fstp=fstep;
fout=0;fm=0;

    for(int i=1;i<Pow2[o-1];i++)               // getting amplitude from compex  number
        { 
              out_r[i]=fastRSS(out_r[i],out_im[i]);
   // Approx RSS function used to calculated magnitude quickly
        
out_im[i]=out_im[i-1]+fstp;
if  (fout<out_r[i]){fm=i; fout=out_r[i];}
         /*
         // un comment to  print Amplitudes (1st value (offset) is not printed)
         Serial.print(out_r[i]);  Serial.print("\	"); 
         Serial.print("*2^");Serial.println(scale);  
         */ 
        }


float fa,fb,fc;
fa=out_r[fm-1];
fb=out_r[fm];  
fc=out_r[fm+1];
fstep=(fa*(fm-1)+fb*fm+fc*(fm+1))/(fa+fb+fc);

return(fstep*Frequency/N);
}

//---------------------------------fast  sine/cosine---------------------------------------//

int fast_sine(int Amp,  int th)
{
int temp3,m1,m2;
byte temp1,temp2, test,quad,accuracy;
accuracy=5;    // set it value from 1 to 7, where 7 being most accurate but slowest
               //  accuracy value of 5 recommended for typical applicaiton
while(th>1024){th=th-1024;}   // here 1024 = 2*pi or 360 deg
while(th<0){th=th+1024;}
quad=th>>8;

       if(quad==1){th= 512-th;}
  else if(quad==2){th= th-512;}
  else if(quad==3){th=  1024-th;}

temp1= 0;
temp2= 128;     //2 multiple
m1=0;
m2=Amp;

    temp3=(m1+m2)>>1;
    Amp=temp3;
      for(int i=0;i<accuracy;i++)
        { test=(temp1+temp2)>>1;
          temp3=temp3>>1; 
          if(th>isin_data[test]){temp1=test;  Amp=Amp+temp3; m1=Amp;}
          else if(th<isin_data[test]){temp2=test; Amp=Amp-temp3;  m2=Amp;}
        }

         if(quad==2){Amp= 0-Amp;}
    else if(quad==3){Amp=  0-Amp;}
return(Amp);
}

int fast_cosine(int Amp, int th)
  {  
  th=256-th;  //cos th = sin (90-th) formula
  return(fast_sine(Amp,th)); 
  }

//--------------------------------------------------------------------------------//


//--------------------------------Fast  RSS----------------------------------------//
int fastRSS(int a, int b)
{  if(a==0 && b==0){return(0);}
  int min,max,temp1,temp2;
  byte clevel;
  if(a<0){a=-a;}
  if(b<0){b=-b;}
clevel=0;
if(a>b){max=a;min=b;} else{max=b;min=a;}

  if(max>(min+min+min))
    {return max;}
  else
    {
     temp1=min>>3;  if(temp1==0){temp1=1;}
     temp2=min;
     while(temp2<max){temp2=temp2+temp1;clevel=clevel+1;}
     temp2=RSSdata[clevel];temp1=temp1>>1;  
     for(int i=0;i<temp2;i++){max=max+temp1;}
     return(max);
    }
}
//--------------------------------------------------------------------------------//

ApproxFFT.ino

arduino

ApproxFFT

arduino

This code performs FFT with good speed

1//---------------------------------lookup data------------------------------------//
2byte
3  isin_data[128]=
4{0,  1,   3,   4,   5,   6,   8,   9,   10,  11,  13,  14,  15,
5  17,  18,  19,  20, 
622,  23,  24,  26,  27,  28,  29,  31,  32,  33,  35,  36,
7  37,  39,  40,  41,  42, 
844,  45,  46,  48,  49,  50,  52,  53,  54,  56,  57,
9  59,  60,  61,  63,  64,  65, 
1067,  68,  70,  71,  72,  74,  75,  77,  78,  80,
11  81,  82,  84,  85,  87,  88,  90, 
1291,  93,  94,  96,  97,  99,  100, 102,
13  104, 105, 107, 108, 110, 112, 113, 115, 117, 
14118, 120, 122, 124, 125, 127, 129,
15  131, 133, 134, 136, 138, 140, 142, 144, 146, 148, 
16150, 152, 155, 157, 159, 161,
17  164, 166, 169, 171, 174, 176, 179, 182, 185, 188, 191, 
18195, 198, 202, 206, 210,
19  215, 221, 227, 236};
20unsigned int Pow2[14]={1,2,4,8,16,32,64,128,256,512,1024,2048,4096};
21byte
22  RSSdata[20]={7,6,6,5,5,5,4,4,4,4,3,3,3,3,3,3,3,2,2,2};
23//---------------------------------------------------------------------------------//
24
25
26//int
27  data[256]={};
28
29
30void setup() 
31        {
32     Serial.begin(250000);
33  
34        }
35
36
37void loop() {
38
39//float f=Approx_FFT(data,512,100);
40//Serial.println(f);
41
42            }
43
44
45//-----------------------------FFT Function----------------------------------------------//
46/*
47Code
48  to perform High speed and Accurate FFT on arduino,
49setup:
50
511. in[]     :
52  Data array, 
532. N        : Number of sample (recommended sample size 2,4,8,16,32,64,128,256,512...)
543.
55  Frequency: sampling frequency required as input (Hz)
56
57It will by default return
58  frequency with max aplitude,
59if you need complex output or magnitudes uncomment
60  required sections
61
62If sample size is not in power of 2 it will be clipped
63  to lower side of number. 
64i.e, for 150 number of samples, code will consider
65  first 128 sample, remaining sample  will be omitted.
66For Arduino nano, FFT of
67  more than 256 sample not possible due to mamory limitation 
68Code by ABHILASH
69Contact:
70  abhilashpatel121@gmail.com
71Documentation & details: https://www.instructables.com/member/abhilash_patel/instructables/
72*/
73
74float
75  Approx_FFT(int in[],int N,float Frequency)
76{ 
77int a,c1,f,o,x,data_max,data_min=0;
78long
79  data_avg,data_mag,temp11;         
80byte scale,check=0;
81
82data_max=0;
83data_avg=0;
84data_min=0;
85
86
87      for(int i=0;i<12;i++)                 //calculating the levels
88         {
89  if(Pow2[i]<=N){o=i;} }
90     a=Pow2[o];  
91int out_r[a];   //real part of transform
92int
93  out_im[a];  //imaginory part of transform
94
95      for(int i=0;i<a;i++)                //getting
96  min max and average for scalling
97          { out_r[i]=0; out_r[i]=0;
98            data_avg=data_avg+in[i];
99
100            if(in[i]>data_max){data_max=in[i];}
101            if(in[i]<data_min){data_min=in[i];}
102
103          }
104
105data_avg=data_avg>>o;
106scale=0;
107data_mag=data_max-data_min;
108temp11=data_mag;
109
110for(int
111  i;i<128;i++)              //scalling data  from +512 to -512
112
113      if(data_mag>1024)
114
115          {while(temp11>1024)
116                {temp11=temp11>>1;
117                scale=scale+1;
118
119                }   
120          }
121 
122      if(data_mag<1024)
123          {while(temp11<1024)
124
125                {temp11=temp11<<1;
126                scale=scale+1;
127                }
128
129          }
130
131
132      if(data_mag>1024)
133          {
134              for(int
135  i=0;i<a;i++)
136                    { in[i]=in[i]-data_avg;
137                      in[i]=in[i]>>scale;
138
139                    }
140                    scale=128-scale;
141          }
142
143
144      if(data_mag<1024)
145          { scale=scale-1;
146            for(int i=0;i<a;i++)
147
148                    {
149                      in[i]=in[i]-data_avg;
150                      in[i]=in[i]<<scale;
151
152                    }
153
154                    scale=128+scale;
155          }
156
157
158       
159x=0;  
160      for(int b=0;b<o;b++)                     // bit reversal
161  order stored in im_out array
162         {
163          c1=Pow2[b];
164          f=Pow2[o]/(c1+c1);
165
166                for(int j=0;j<c1;j++)
167                    { 
168                     x=x+1;
169
170                     out_im[x]=out_im[j]+f;
171                    }
172         }
173
174
175  
176      for(int i=0;i<a;i++)            // update input array as per bit reverse
177  order
178         {
179          out_r[i]=in[out_im[i]]; 
180          out_im[i]=0;
181
182         }
183
184
185int i10,i11,n1,tr,ti;
186float e;
187int c,s,temp4;
188    for(int
189  i=0;i<o;i++)                                    //fft
190    {
191     i10=Pow2[i];
192              // overall values of sine/cosine  
193     i11=Pow2[o]/Pow2[i+1];
194    // loop with similar sine cosine
195     e=1024/Pow2[i+1];  //1024 is equivalent
196  to 360 deg
197     e=0-e;
198     n1=0;
199
200          for(int j=0;j<i10;j++)
201
202          {
203            c=e*j;    //c is angle as where 1024 unit is 360 deg
204
205  while(c<0){c=c+1024;}
206  while(c>1024){c=c-1024;}
207
208          n1=j;
209
210          
211          for(int k=0;k<i11;k++)
212                 {
213                   temp4=i10+n1;
214
215       if(c==0)   {tr=out_r[temp4];
216                   ti=out_im[temp4];}
217
218  else if(c==256) {tr= -out_im[temp4];
219                   ti=out_r[temp4];}
220
221  else if(c==512) {tr=-out_r[temp4];
222                  ti=-out_im[temp4];}
223
224  else if(c==768) {tr=out_im[temp4];
225                   ti=-out_r[temp4];}
226
227  else if(c==1024){tr=out_r[temp4];
228                   ti=out_im[temp4];}
229
230  else{
231    tr=fast_cosine(out_r[temp4],c)-fast_sine(out_im[temp4],c);            //the
232  fast sine/cosine function gives direct (approx) output for A*sinx
233    ti=fast_sine(out_r[temp4],c)+fast_cosine(out_im[temp4],c);
234            
235      }
236          
237                 out_r[n1+i10]=out_r[n1]-tr;
238
239                 out_r[n1]=out_r[n1]+tr;
240                 if(out_r[n1]>15000
241  || out_r[n1]<-15000){check=1;}   //check for int size, it can handle only +31000
242  to -31000,
243          
244                 out_im[n1+i10]=out_im[n1]-ti;
245                 out_im[n1]=out_im[n1]+ti;
246
247                 if(out_im[n1]>15000 || out_im[n1]<-15000){check=1;}          
248
249          
250                 n1=n1+i10+i10;
251                  }       
252
253             }
254
255    if(check==1){                                             //
256  scalling the matrics if value higher than 15000 to prevent varible from overflowing
257
258                for(int i=0;i<a;i++)
259                    {
260                     out_r[i]=out_r[i]>>1;
261           
262                     out_im[i]=out_im[i]>>1; 
263                    }
264
265                     check=0; 
266                     scale=scale-1;                 //
267  tracking overall scalling of input data
268                }           
269
270
271     }
272
273
274if(scale>128)
275    {scale=scale-128;
276     for(int i=0;i<a;i++)
277
278      {out_r[i]=out_r[i]>>scale;
279       out_im[i]=out_im[i]>>scale;
280      }
281
282      scale=0;
283    }                                                   // revers
284  all scalling we done till here,
285else{scale=128-scale;}                             //
286  in case of nnumber getting higher than 32000, we will represent in as multiple of
287  2^scale
288
289/*
290for(int i=0;i<a;i++)
291{
292Serial.print(out_r[i]);Serial.print("\	");
293                    // un comment to print RAW o/p    
294Serial.print(out_im[i]);
295  
296Serial.print("i");Serial.print("\	"); 
297Serial.print("*2^");Serial.println(scale);
298  
299}
300*/
301
302//---> here onward out_r contains amplitude and our_in conntains
303  frequency (Hz)
304int fout,fm,fstp;
305float fstep;
306fstep=Frequency/N;
307fstp=fstep;
308fout=0;fm=0;
309
310
311    for(int i=1;i<Pow2[o-1];i++)               // getting amplitude from compex
312  number
313        { 
314              out_r[i]=fastRSS(out_r[i],out_im[i]);
315
316   // Approx RSS function used to calculated magnitude quickly
317        
318out_im[i]=out_im[i-1]+fstp;
319if
320  (fout<out_r[i]){fm=i; fout=out_r[i];}
321         /*
322         // un comment to
323  print Amplitudes (1st value (offset) is not printed)
324         Serial.print(out_r[i]);
325  Serial.print("\	"); 
326         Serial.print("*2^");Serial.println(scale);
327  
328         */ 
329        }
330
331
332float fa,fb,fc;
333fa=out_r[fm-1];
334fb=out_r[fm];
335  
336fc=out_r[fm+1];
337fstep=(fa*(fm-1)+fb*fm+fc*(fm+1))/(fa+fb+fc);
338
339return(fstep*Frequency/N);
340}
341
342//---------------------------------fast
343  sine/cosine---------------------------------------//
344
345int fast_sine(int Amp,
346  int th)
347{
348int temp3,m1,m2;
349byte temp1,temp2, test,quad,accuracy;
350accuracy=5;
351    // set it value from 1 to 7, where 7 being most accurate but slowest
352               //
353  accuracy value of 5 recommended for typical applicaiton
354while(th>1024){th=th-1024;}
355   // here 1024 = 2*pi or 360 deg
356while(th<0){th=th+1024;}
357quad=th>>8;
358
359
360       if(quad==1){th= 512-th;}
361  else if(quad==2){th= th-512;}
362  else if(quad==3){th=
363  1024-th;}
364
365temp1= 0;
366temp2= 128;     //2 multiple
367m1=0;
368m2=Amp;
369
370
371    temp3=(m1+m2)>>1;
372    Amp=temp3;
373      for(int i=0;i<accuracy;i++)
374
375        { test=(temp1+temp2)>>1;
376          temp3=temp3>>1; 
377          if(th>isin_data[test]){temp1=test;
378  Amp=Amp+temp3; m1=Amp;}
379          else if(th<isin_data[test]){temp2=test; Amp=Amp-temp3;
380  m2=Amp;}
381        }
382
383         if(quad==2){Amp= 0-Amp;}
384    else if(quad==3){Amp=
385  0-Amp;}
386return(Amp);
387}
388
389int fast_cosine(int Amp, int th)
390  {  
391
392  th=256-th;  //cos th = sin (90-th) formula
393  return(fast_sine(Amp,th)); 
394
395  }
396
397//--------------------------------------------------------------------------------//
398
399
400//--------------------------------Fast
401  RSS----------------------------------------//
402int fastRSS(int a, int b)
403{
404  if(a==0 && b==0){return(0);}
405  int min,max,temp1,temp2;
406  byte clevel;
407
408  if(a<0){a=-a;}
409  if(b<0){b=-b;}
410clevel=0;
411if(a>b){max=a;min=b;} else{max=b;min=a;}
412
413
414  if(max>(min+min+min))
415    {return max;}
416  else
417    {
418     temp1=min>>3;
419  if(temp1==0){temp1=1;}
420     temp2=min;
421     while(temp2<max){temp2=temp2+temp1;clevel=clevel+1;}
422
423     temp2=RSSdata[clevel];temp1=temp1>>1;  
424     for(int i=0;i<temp2;i++){max=max+temp1;}
425
426     return(max);
427    }
428}
429//--------------------------------------------------------------------------------//

//---------------------------------lookup data------------------------------------//
byte
  isin_data[128]=
{0,  1,   3,   4,   5,   6,   8,   9,   10,  11,  13,  14,  15,
  17,  18,  19,  20, 
22,  23,  24,  26,  27,  28,  29,  31,  32,  33,  35,  36,
  37,  39,  40,  41,  42, 
44,  45,  46,  48,  49,  50,  52,  53,  54,  56,  57,
  59,  60,  61,  63,  64,  65, 
67,  68,  70,  71,  72,  74,  75,  77,  78,  80,
  81,  82,  84,  85,  87,  88,  90, 
91,  93,  94,  96,  97,  99,  100, 102,
  104, 105, 107, 108, 110, 112, 113, 115, 117, 
118, 120, 122, 124, 125, 127, 129,
  131, 133, 134, 136, 138, 140, 142, 144, 146, 148, 
150, 152, 155, 157, 159, 161,
  164, 166, 169, 171, 174, 176, 179, 182, 185, 188, 191, 
195, 198, 202, 206, 210,
  215, 221, 227, 236};
unsigned int Pow2[14]={1,2,4,8,16,32,64,128,256,512,1024,2048,4096};
byte
  RSSdata[20]={7,6,6,5,5,5,4,4,4,4,3,3,3,3,3,3,3,2,2,2};
//---------------------------------------------------------------------------------//


//int
  data[256]={};


void setup() 
        {
     Serial.begin(250000);
  
        }


void loop() {

//float f=Approx_FFT(data,512,100);
//Serial.println(f);

            }


//-----------------------------FFT Function----------------------------------------------//
/*
Code
  to perform High speed and Accurate FFT on arduino,
setup:

1. in[]     :
  Data array, 
2. N        : Number of sample (recommended sample size 2,4,8,16,32,64,128,256,512...)
3.
  Frequency: sampling frequency required as input (Hz)

It will by default return
  frequency with max aplitude,
if you need complex output or magnitudes uncomment
  required sections

If sample size is not in power of 2 it will be clipped
  to lower side of number. 
i.e, for 150 number of samples, code will consider
  first 128 sample, remaining sample  will be omitted.
For Arduino nano, FFT of
  more than 256 sample not possible due to mamory limitation 
Code by ABHILASH
Contact:
  abhilashpatel121@gmail.com
Documentation & details: https://www.instructables.com/member/abhilash_patel/instructables/
*/

float
  Approx_FFT(int in[],int N,float Frequency)
{ 
int a,c1,f,o,x,data_max,data_min=0;
long
  data_avg,data_mag,temp11;         
byte scale,check=0;

data_max=0;
data_avg=0;
data_min=0;


      for(int i=0;i<12;i++)                 //calculating the levels
         {
  if(Pow2[i]<=N){o=i;} }
     a=Pow2[o];  
int out_r[a];   //real part of transform
int
  out_im[a];  //imaginory part of transform

      for(int i=0;i<a;i++)                //getting
  min max and average for scalling
          { out_r[i]=0; out_r[i]=0;
            data_avg=data_avg+in[i];

            if(in[i]>data_max){data_max=in[i];}
            if(in[i]<data_min){data_min=in[i];}

          }

data_avg=data_avg>>o;
scale=0;
data_mag=data_max-data_min;
temp11=data_mag;

for(int
  i;i<128;i++)              //scalling data  from +512 to -512

      if(data_mag>1024)

          {while(temp11>1024)
                {temp11=temp11>>1;
                scale=scale+1;

                }   
          }
 
      if(data_mag<1024)
          {while(temp11<1024)

                {temp11=temp11<<1;
                scale=scale+1;
                }

          }


      if(data_mag>1024)
          {
              for(int
  i=0;i<a;i++)
                    { in[i]=in[i]-data_avg;
                      in[i]=in[i]>>scale;

                    }
                    scale=128-scale;
          }


      if(data_mag<1024)
          { scale=scale-1;
            for(int i=0;i<a;i++)

                    {
                      in[i]=in[i]-data_avg;
                      in[i]=in[i]<<scale;

                    }

                    scale=128+scale;
          }


       
x=0;  
      for(int b=0;b<o;b++)                     // bit reversal
  order stored in im_out array
         {
          c1=Pow2[b];
          f=Pow2[o]/(c1+c1);

                for(int j=0;j<c1;j++)
                    { 
                     x=x+1;

                     out_im[x]=out_im[j]+f;
                    }
         }


  
      for(int i=0;i<a;i++)            // update input array as per bit reverse
  order
         {
          out_r[i]=in[out_im[i]]; 
          out_im[i]=0;

         }


int i10,i11,n1,tr,ti;
float e;
int c,s,temp4;
    for(int
  i=0;i<o;i++)                                    //fft
    {
     i10=Pow2[i];
              // overall values of sine/cosine  
     i11=Pow2[o]/Pow2[i+1];
    // loop with similar sine cosine
     e=1024/Pow2[i+1];  //1024 is equivalent
  to 360 deg
     e=0-e;
     n1=0;

          for(int j=0;j<i10;j++)

          {
            c=e*j;    //c is angle as where 1024 unit is 360 deg

  while(c<0){c=c+1024;}
  while(c>1024){c=c-1024;}

          n1=j;

          
          for(int k=0;k<i11;k++)
                 {
                   temp4=i10+n1;

       if(c==0)   {tr=out_r[temp4];
                   ti=out_im[temp4];}

  else if(c==256) {tr= -out_im[temp4];
                   ti=out_r[temp4];}

  else if(c==512) {tr=-out_r[temp4];
                  ti=-out_im[temp4];}

  else if(c==768) {tr=out_im[temp4];
                   ti=-out_r[temp4];}

  else if(c==1024){tr=out_r[temp4];
                   ti=out_im[temp4];}

  else{
    tr=fast_cosine(out_r[temp4],c)-fast_sine(out_im[temp4],c);            //the
  fast sine/cosine function gives direct (approx) output for A*sinx
    ti=fast_sine(out_r[temp4],c)+fast_cosine(out_im[temp4],c);
            
      }
          
                 out_r[n1+i10]=out_r[n1]-tr;

                 out_r[n1]=out_r[n1]+tr;
                 if(out_r[n1]>15000
  || out_r[n1]<-15000){check=1;}   //check for int size, it can handle only +31000
  to -31000,
          
                 out_im[n1+i10]=out_im[n1]-ti;
                 out_im[n1]=out_im[n1]+ti;

                 if(out_im[n1]>15000 || out_im[n1]<-15000){check=1;}          

          
                 n1=n1+i10+i10;
                  }       

             }

    if(check==1){                                             //
  scalling the matrics if value higher than 15000 to prevent varible from overflowing

                for(int i=0;i<a;i++)
                    {
                     out_r[i]=out_r[i]>>1;
           
                     out_im[i]=out_im[i]>>1; 
                    }

                     check=0; 
                     scale=scale-1;                 //
  tracking overall scalling of input data
                }           


     }


if(scale>128)
    {scale=scale-128;
     for(int i=0;i<a;i++)

      {out_r[i]=out_r[i]>>scale;
       out_im[i]=out_im[i]>>scale;
      }

      scale=0;
    }                                                   // revers
  all scalling we done till here,
else{scale=128-scale;}                             //
  in case of nnumber getting higher than 32000, we will represent in as multiple of
  2^scale

/*
for(int i=0;i<a;i++)
{
Serial.print(out_r[i]);Serial.print("\	");
                    // un comment to print RAW o/p    
Serial.print(out_im[i]);
  
Serial.print("i");Serial.print("\	"); 
Serial.print("*2^");Serial.println(scale);
  
}
*/

//---> here onward out_r contains amplitude and our_in conntains
  frequency (Hz)
int fout,fm,fstp;
float fstep;
fstep=Frequency/N;
fstp=fstep;
fout=0;fm=0;


    for(int i=1;i<Pow2[o-1];i++)               // getting amplitude from compex
  number
        { 
              out_r[i]=fastRSS(out_r[i],out_im[i]);

   // Approx RSS function used to calculated magnitude quickly
        
out_im[i]=out_im[i-1]+fstp;
if
  (fout<out_r[i]){fm=i; fout=out_r[i];}
         /*
         // un comment to
  print Amplitudes (1st value (offset) is not printed)
         Serial.print(out_r[i]);
  Serial.print("\	"); 
         Serial.print("*2^");Serial.println(scale);
  
         */ 
        }


float fa,fb,fc;
fa=out_r[fm-1];
fb=out_r[fm];
  
fc=out_r[fm+1];
fstep=(fa*(fm-1)+fb*fm+fc*(fm+1))/(fa+fb+fc);

return(fstep*Frequency/N);
}

//---------------------------------fast
  sine/cosine---------------------------------------//

int fast_sine(int Amp,
  int th)
{
int temp3,m1,m2;
byte temp1,temp2, test,quad,accuracy;
accuracy=5;
    // set it value from 1 to 7, where 7 being most accurate but slowest
               //
  accuracy value of 5 recommended for typical applicaiton
while(th>1024){th=th-1024;}
   // here 1024 = 2*pi or 360 deg
while(th<0){th=th+1024;}
quad=th>>8;


       if(quad==1){th= 512-th;}
  else if(quad==2){th= th-512;}
  else if(quad==3){th=
  1024-th;}

temp1= 0;
temp2= 128;     //2 multiple
m1=0;
m2=Amp;


    temp3=(m1+m2)>>1;
    Amp=temp3;
      for(int i=0;i<accuracy;i++)

        { test=(temp1+temp2)>>1;
          temp3=temp3>>1; 
          if(th>isin_data[test]){temp1=test;
  Amp=Amp+temp3; m1=Amp;}
          else if(th<isin_data[test]){temp2=test; Amp=Amp-temp3;
  m2=Amp;}
        }

         if(quad==2){Amp= 0-Amp;}
    else if(quad==3){Amp=
  0-Amp;}
return(Amp);
}

int fast_cosine(int Amp, int th)
  {  

  th=256-th;  //cos th = sin (90-th) formula
  return(fast_sine(Amp,th)); 

  }

//--------------------------------------------------------------------------------//


//--------------------------------Fast
  RSS----------------------------------------//
int fastRSS(int a, int b)
{
  if(a==0 && b==0){return(0);}
  int min,max,temp1,temp2;
  byte clevel;

  if(a<0){a=-a;}
  if(b<0){b=-b;}
clevel=0;
if(a>b){max=a;min=b;} else{max=b;min=a;}


  if(max>(min+min+min))
    {return max;}
  else
    {
     temp1=min>>3;
  if(temp1==0){temp1=1;}
     temp2=min;
     while(temp2<max){temp2=temp2+temp1;clevel=clevel+1;}

     temp2=RSSdata[clevel];temp1=temp1>>1;  
     for(int i=0;i<temp2;i++){max=max+temp1;}

     return(max);
    }
}
//--------------------------------------------------------------------------------//

ApproxFFT.ino

arduino

ApproxFFT_V2_060521

arduino

1
2//---------------------------------lookup data------------------------------------//
3byte  isin_data[128]=
4{0,  1,   3,   4,   5,   6,   8,   9,   10,  11,  13,  14,  15,  17,  18,  19,  20, 
522,  23,  24,  26,  27,  28,  29,  31,  32,  33,  35,  36,  37,  39,  40,  41,  42, 
644,  45,  46,  48,  49,  50,  52,  53,  54,  56,  57,  59,  60,  61,  63,  64,  65, 
767,  68,  70,  71,  72,  74,  75,  77,  78,  80,  81,  82,  84,  85,  87,  88,  90, 
891,  93,  94,  96,  97,  99,  100, 102,  104, 105, 107, 108, 110, 112, 113, 115, 117, 
9118, 120, 122, 124, 125, 127, 129,  131, 133, 134, 136, 138, 140, 142, 144, 146, 148, 
10150, 152, 155, 157, 159, 161,  164, 166, 169, 171, 174, 176, 179, 182, 185, 188, 191, 
11195, 198, 202, 206, 210,  215, 221, 227, 236};
12unsigned int Pow2[14]={1,2,4,8,16,32,64,128,256,512,1024,2048,4096};
13byte  RSSdata[20]={7,6,6,5,5,5,4,4,4,4,3,3,3,3,3,3,3,2,2,2};
14//---------------------------------------------------------------------------------//
15
16
17
18void  setup() 
19        {
20     Serial.begin(250000);  
21        }
22
23void loop()  {
24
25//float f=Approx_FFT(data,256,100);
26//Serial.println(f);
27//delay(99999);
28            }
29
30
31
32
33//-----------------------------FFT Function----------------------------------------------//
34/*
35Code  to perform High speed and Accurate FFT on arduino,
36setup:
37
381. in[]     :  Data array, 
392. N        : Number of sample (recommended sample size 2,4,8,16,32,64,128,256,512...)
403.  Frequency: sampling frequency required as input (Hz)
41
42It will by default return  frequency with max aplitude,
43if you need complex output or magnitudes uncomment  required sections
44
45If sample size is not in power of 2 it will be clipped  to lower side of number. 
46i.e, for 150 number of samples, code will consider  first 128 sample, remaining sample  will be omitted.
47For Arduino nano, FFT of  more than 256 sample not possible due to mamory limitation 
48Code by ABHILASH
49Contact:  abhilashpatel121@gmail.com
50Documentation & details: https://www.instructables.com/member/abhilash_patel/instructables/
51
52Update(06/05/21):  Correction made for support on Arduino Due
53*/
54
55float Approx_FFT(int in[],int  N,float Frequency)
56{ 
57int a,c1,f,o,x,data_max,data_min=0;
58long data_avg,data_mag,temp11;         
59byte scale,check=0;
60
61data_max=0;
62data_avg=0;
63data_min=0;
64
65      for(int i=0;i<12;i++)                 //calculating the levels
66         {  if(Pow2[i]<=N){o=i;}}
67     a=Pow2[o];  
68int out_r[a];   //real part of transform
69int  out_im[a];  //imaginory part of transform
70
71      for(int i=0;i<a;i++)                //getting  min max and average for scalling
72          { out_r[i]=0; out_im[i]=0;
73            data_avg=data_avg+in[i];
74            if(in[i]>data_max){data_max=in[i];}
75            if(in[i]<data_min){data_min=in[i];}  
76          }
77
78data_avg=data_avg>>o;
79scale=0;
80data_mag=data_max-data_min;
81temp11=data_mag;
82
83  //scalling data  from +512 to -512
84
85      if(data_mag>1024)
86          {while(temp11>1024)
87                {temp11=temp11>>1;
88                scale=scale+1;
89                }   
90          }
91 
92      if(data_mag<1024)
93          {while(temp11<1024)
94                {temp11=temp11<<1;
95                scale=scale+1;
96                }
97          }
98
99
100      if(data_mag>1024)
101          {
102              for(int  i=0;i<a;i++)
103                    { in[i]=in[i]-data_avg;
104                      in[i]=in[i]>>scale;
105                    }
106                    scale=128-scale;
107          }
108
109      if(data_mag<1024)
110          { scale=scale-1;
111            for(int i=0;i<a;i++)
112                    {
113                      in[i]=in[i]-data_avg;
114                      in[i]=in[i]<<scale;
115                    }
116
117                    scale=128+scale;
118          }
119
120
121x=0;  
122      for(int b=0;b<o;b++)                     // bit reversal order stored  in im_out array
123         {
124          c1=Pow2[b];
125          f=Pow2[o]/(c1+c1);
126                for(int j=0;j<c1;j++)
127                    { 
128                     x=x+1;
129                     out_im[x]=out_im[j]+f;
130                    }
131         }
132
133      for(int i=0;i<a;i++)            // update input array as per bit reverse order
134         {
135          out_r[i]=in[out_im[i]]; 
136          out_im[i]=0;
137         }
138
139
140int  i10,i11,n1,tr,ti;
141float e;
142int c,s,temp4;
143    for(int i=0;i<o;i++)                                    //fft
144    {
145     i10=Pow2[i];              // overall values of sine/cosine  
146     i11=Pow2[o]/Pow2[i+1];    // loop with similar sine cosine
147     e=1024/Pow2[i+1];  //1024 is equivalent  to 360 deg
148     e=0-e;
149     n1=0;
150
151          for(int j=0;j<i10;j++)
152          {
153            c=e*j;    //c is angle as where 1024 unit is 360 deg
154  while(c<0){c=c+1024;}
155  while(c>1024){c=c-1024;}
156
157          n1=j;
158          
159          for(int k=0;k<i11;k++)
160                 {
161                   temp4=i10+n1;
162       if(c==0)   {tr=out_r[temp4];
163                   ti=out_im[temp4];}
164  else if(c==256) {tr= -out_im[temp4];
165                   ti=out_r[temp4];}
166  else if(c==512) {tr=-out_r[temp4];
167                  ti=-out_im[temp4];}
168  else if(c==768) {tr=out_im[temp4];
169                   ti=-out_r[temp4];}
170  else if(c==1024){tr=out_r[temp4];
171                   ti=out_im[temp4];}
172  else{
173    tr=fast_cosine(out_r[temp4],c)-fast_sine(out_im[temp4],c);            //the  fast sine/cosine function gives direct (approx) output for A*sinx
174    ti=fast_sine(out_r[temp4],c)+fast_cosine(out_im[temp4],c);            
175      }
176          
177                 out_r[n1+i10]=out_r[n1]-tr;
178                 out_r[n1]=out_r[n1]+tr;
179                 if(out_r[n1]>15000  || out_r[n1]<-15000){check=1;}   //check for int size, it can handle only +31000  to -31000,
180          
181                 out_im[n1+i10]=out_im[n1]-ti;
182                 out_im[n1]=out_im[n1]+ti;
183                 if(out_im[n1]>15000 || out_im[n1]<-15000){check=1;}          
184          
185                 n1=n1+i10+i10;
186                  }       
187             }
188
189    if(check==1){                                             //  scalling the matrics if value higher than 15000 to prevent varible from overflowing
190                for(int i=0;i<a;i++)
191                    {
192                     out_r[i]=out_r[i]>>1;           
193                     out_im[i]=out_im[i]>>1; 
194                    }
195                     check=0; 
196                     scale=scale-1;                 //  tracking overall scalling of input data
197                }           
198
199     }
200
201
202if(scale>128)
203    {scale=scale-128;
204     for(int i=0;i<a;i++)
205      {out_r[i]=out_r[i]>>scale;
206       out_im[i]=out_im[i]>>scale;
207      }
208      scale=0;
209    }                                                   // revers  all scalling we done till here,
210else{scale=128-scale;}                             //  in case of nnumber getting higher than 32000, we will represent in as multiple of  2^scale
211
212/*
213for(int i=0;i<a;i++)
214{
215Serial.print(out_r[i]);Serial.print("\	");                    // un comment to print RAW o/p    
216Serial.print(out_im[i]);  
217Serial.print("i");Serial.print("\	"); 
218Serial.print("*2^");Serial.println(scale);  
219}
220*/
221
222//---> here onward out_r contains amplitude and our_in conntains  frequency (Hz)
223int fout,fm,fstp;
224float fstep;
225fstep=Frequency/N;
226fstp=fstep;
227fout=0;fm=0;
228
229    for(int i=1;i<Pow2[o-1];i++)               // getting amplitude from compex  number
230        { 
231              out_r[i]=fastRSS(out_r[i],out_im[i]);
232   // Approx RSS function used to calculated magnitude quickly
233        
234out_im[i]=out_im[i-1]+fstp;
235if  (fout<out_r[i]){fm=i; fout=out_r[i];}
236         
237         // un comment to  print Amplitudes (1st value (offset) is not printed)
238         //Serial.print(out_r[i]);  Serial.print("\	"); 
239         //Serial.print("*2^");Serial.println(scale);  
240        }
241
242
243float fa,fb,fc;
244fa=out_r[fm-1];
245fb=out_r[fm]; 
246fc=out_r[fm+1];
247fstep=(fa*(fm-1)+fb*fm+fc*(fm+1))/(fa+fb+fc);
248
249return(fstep*Frequency/N);
250}
251
252//---------------------------------fast  sine/cosine---------------------------------------//
253
254int fast_sine(int Amp,  int th)
255{
256int temp3,m1,m2;
257byte temp1,temp2, test,quad,accuracy;
258accuracy=5;    // set it value from 1 to 7, where 7 being most accurate but slowest
259               //  accuracy value of 5 recommended for typical applicaiton
260while(th>1024){th=th-1024;}   // here 1024 = 2*pi or 360 deg
261while(th<0){th=th+1024;}
262quad=th>>8;
263
264       if(quad==1){th= 512-th;}
265  else if(quad==2){th= th-512;}
266  else if(quad==3){th=  1024-th;}
267
268temp1= 0;
269temp2= 128;     //2 multiple
270m1=0;
271m2=Amp;
272
273    temp3=(m1+m2)>>1;
274    Amp=temp3;
275      for(int i=0;i<accuracy;i++)
276        { test=(temp1+temp2)>>1;
277          temp3=temp3>>1; 
278          if(th>isin_data[test]){temp1=test;  Amp=Amp+temp3; m1=Amp;}
279          else if(th<isin_data[test]){temp2=test; Amp=Amp-temp3;  m2=Amp;}
280        }
281
282         if(quad==2){Amp= 0-Amp;}
283    else if(quad==3){Amp=  0-Amp;}
284return(Amp);
285}
286
287int fast_cosine(int Amp, int th)
288  {  
289  th=256-th;  //cos th = sin (90-th) formula
290  return(fast_sine(Amp,th)); 
291  }
292
293//--------------------------------------------------------------------------------//
294
295
296//--------------------------------Fast  RSS----------------------------------------//
297int fastRSS(int a, int b)
298{  if(a==0 && b==0){return(0);}
299  int min,max,temp1,temp2;
300  byte clevel;
301  if(a<0){a=-a;}
302  if(b<0){b=-b;}
303clevel=0;
304if(a>b){max=a;min=b;} else{max=b;min=a;}
305
306  if(max>(min+min+min))
307    {return max;}
308  else
309    {
310     temp1=min>>3;  if(temp1==0){temp1=1;}
311     temp2=min;
312     while(temp2<max){temp2=temp2+temp1;clevel=clevel+1;}
313     temp2=RSSdata[clevel];temp1=temp1>>1;  
314     for(int i=0;i<temp2;i++){max=max+temp1;}
315     return(max);
316    }
317}
318//--------------------------------------------------------------------------------//

//---------------------------------lookup data------------------------------------//
byte  isin_data[128]=
{0,  1,   3,   4,   5,   6,   8,   9,   10,  11,  13,  14,  15,  17,  18,  19,  20, 
22,  23,  24,  26,  27,  28,  29,  31,  32,  33,  35,  36,  37,  39,  40,  41,  42, 
44,  45,  46,  48,  49,  50,  52,  53,  54,  56,  57,  59,  60,  61,  63,  64,  65, 
67,  68,  70,  71,  72,  74,  75,  77,  78,  80,  81,  82,  84,  85,  87,  88,  90, 
91,  93,  94,  96,  97,  99,  100, 102,  104, 105, 107, 108, 110, 112, 113, 115, 117, 
118, 120, 122, 124, 125, 127, 129,  131, 133, 134, 136, 138, 140, 142, 144, 146, 148, 
150, 152, 155, 157, 159, 161,  164, 166, 169, 171, 174, 176, 179, 182, 185, 188, 191, 
195, 198, 202, 206, 210,  215, 221, 227, 236};
unsigned int Pow2[14]={1,2,4,8,16,32,64,128,256,512,1024,2048,4096};
byte  RSSdata[20]={7,6,6,5,5,5,4,4,4,4,3,3,3,3,3,3,3,2,2,2};
//---------------------------------------------------------------------------------//



void  setup() 
        {
     Serial.begin(250000);  
        }

void loop()  {

//float f=Approx_FFT(data,256,100);
//Serial.println(f);
//delay(99999);
            }




//-----------------------------FFT Function----------------------------------------------//
/*
Code  to perform High speed and Accurate FFT on arduino,
setup:

1. in[]     :  Data array, 
2. N        : Number of sample (recommended sample size 2,4,8,16,32,64,128,256,512...)
3.  Frequency: sampling frequency required as input (Hz)

It will by default return  frequency with max aplitude,
if you need complex output or magnitudes uncomment  required sections

If sample size is not in power of 2 it will be clipped  to lower side of number. 
i.e, for 150 number of samples, code will consider  first 128 sample, remaining sample  will be omitted.
For Arduino nano, FFT of  more than 256 sample not possible due to mamory limitation 
Code by ABHILASH
Contact:  abhilashpatel121@gmail.com
Documentation & details: https://www.instructables.com/member/abhilash_patel/instructables/

Update(06/05/21):  Correction made for support on Arduino Due
*/

float Approx_FFT(int in[],int  N,float Frequency)
{ 
int a,c1,f,o,x,data_max,data_min=0;
long data_avg,data_mag,temp11;         
byte scale,check=0;

data_max=0;
data_avg=0;
data_min=0;

      for(int i=0;i<12;i++)                 //calculating the levels
         {  if(Pow2[i]<=N){o=i;}}
     a=Pow2[o];  
int out_r[a];   //real part of transform
int  out_im[a];  //imaginory part of transform

      for(int i=0;i<a;i++)                //getting  min max and average for scalling
          { out_r[i]=0; out_im[i]=0;
            data_avg=data_avg+in[i];
            if(in[i]>data_max){data_max=in[i];}
            if(in[i]<data_min){data_min=in[i];}  
          }

data_avg=data_avg>>o;
scale=0;
data_mag=data_max-data_min;
temp11=data_mag;

  //scalling data  from +512 to -512

      if(data_mag>1024)
          {while(temp11>1024)
                {temp11=temp11>>1;
                scale=scale+1;
                }   
          }
 
      if(data_mag<1024)
          {while(temp11<1024)
                {temp11=temp11<<1;
                scale=scale+1;
                }
          }


      if(data_mag>1024)
          {
              for(int  i=0;i<a;i++)
                    { in[i]=in[i]-data_avg;
                      in[i]=in[i]>>scale;
                    }
                    scale=128-scale;
          }

      if(data_mag<1024)
          { scale=scale-1;
            for(int i=0;i<a;i++)
                    {
                      in[i]=in[i]-data_avg;
                      in[i]=in[i]<<scale;
                    }

                    scale=128+scale;
          }


x=0;  
      for(int b=0;b<o;b++)                     // bit reversal order stored  in im_out array
         {
          c1=Pow2[b];
          f=Pow2[o]/(c1+c1);
                for(int j=0;j<c1;j++)
                    { 
                     x=x+1;
                     out_im[x]=out_im[j]+f;
                    }
         }

      for(int i=0;i<a;i++)            // update input array as per bit reverse order
         {
          out_r[i]=in[out_im[i]]; 
          out_im[i]=0;
         }


int  i10,i11,n1,tr,ti;
float e;
int c,s,temp4;
    for(int i=0;i<o;i++)                                    //fft
    {
     i10=Pow2[i];              // overall values of sine/cosine  
     i11=Pow2[o]/Pow2[i+1];    // loop with similar sine cosine
     e=1024/Pow2[i+1];  //1024 is equivalent  to 360 deg
     e=0-e;
     n1=0;

          for(int j=0;j<i10;j++)
          {
            c=e*j;    //c is angle as where 1024 unit is 360 deg
  while(c<0){c=c+1024;}
  while(c>1024){c=c-1024;}

          n1=j;
          
          for(int k=0;k<i11;k++)
                 {
                   temp4=i10+n1;
       if(c==0)   {tr=out_r[temp4];
                   ti=out_im[temp4];}
  else if(c==256) {tr= -out_im[temp4];
                   ti=out_r[temp4];}
  else if(c==512) {tr=-out_r[temp4];
                  ti=-out_im[temp4];}
  else if(c==768) {tr=out_im[temp4];
                   ti=-out_r[temp4];}
  else if(c==1024){tr=out_r[temp4];
                   ti=out_im[temp4];}
  else{
    tr=fast_cosine(out_r[temp4],c)-fast_sine(out_im[temp4],c);            //the  fast sine/cosine function gives direct (approx) output for A*sinx
    ti=fast_sine(out_r[temp4],c)+fast_cosine(out_im[temp4],c);            
      }
          
                 out_r[n1+i10]=out_r[n1]-tr;
                 out_r[n1]=out_r[n1]+tr;
                 if(out_r[n1]>15000  || out_r[n1]<-15000){check=1;}   //check for int size, it can handle only +31000  to -31000,
          
                 out_im[n1+i10]=out_im[n1]-ti;
                 out_im[n1]=out_im[n1]+ti;
                 if(out_im[n1]>15000 || out_im[n1]<-15000){check=1;}          
          
                 n1=n1+i10+i10;
                  }       
             }

    if(check==1){                                             //  scalling the matrics if value higher than 15000 to prevent varible from overflowing
                for(int i=0;i<a;i++)
                    {
                     out_r[i]=out_r[i]>>1;           
                     out_im[i]=out_im[i]>>1; 
                    }
                     check=0; 
                     scale=scale-1;                 //  tracking overall scalling of input data
                }           

     }


if(scale>128)
    {scale=scale-128;
     for(int i=0;i<a;i++)
      {out_r[i]=out_r[i]>>scale;
       out_im[i]=out_im[i]>>scale;
      }
      scale=0;
    }                                                   // revers  all scalling we done till here,
else{scale=128-scale;}                             //  in case of nnumber getting higher than 32000, we will represent in as multiple of  2^scale

/*
for(int i=0;i<a;i++)
{
Serial.print(out_r[i]);Serial.print("\	");                    // un comment to print RAW o/p    
Serial.print(out_im[i]);  
Serial.print("i");Serial.print("\	"); 
Serial.print("*2^");Serial.println(scale);  
}
*/

//---> here onward out_r contains amplitude and our_in conntains  frequency (Hz)
int fout,fm,fstp;
float fstep;
fstep=Frequency/N;
fstp=fstep;
fout=0;fm=0;

    for(int i=1;i<Pow2[o-1];i++)               // getting amplitude from compex  number
        { 
              out_r[i]=fastRSS(out_r[i],out_im[i]);
   // Approx RSS function used to calculated magnitude quickly
        
out_im[i]=out_im[i-1]+fstp;
if  (fout<out_r[i]){fm=i; fout=out_r[i];}
         
         // un comment to  print Amplitudes (1st value (offset) is not printed)
         //Serial.print(out_r[i]);  Serial.print("\	"); 
         //Serial.print("*2^");Serial.println(scale);  
        }


float fa,fb,fc;
fa=out_r[fm-1];
fb=out_r[fm]; 
fc=out_r[fm+1];
fstep=(fa*(fm-1)+fb*fm+fc*(fm+1))/(fa+fb+fc);

return(fstep*Frequency/N);
}

//---------------------------------fast  sine/cosine---------------------------------------//

int fast_sine(int Amp,  int th)
{
int temp3,m1,m2;
byte temp1,temp2, test,quad,accuracy;
accuracy=5;    // set it value from 1 to 7, where 7 being most accurate but slowest
               //  accuracy value of 5 recommended for typical applicaiton
while(th>1024){th=th-1024;}   // here 1024 = 2*pi or 360 deg
while(th<0){th=th+1024;}
quad=th>>8;

       if(quad==1){th= 512-th;}
  else if(quad==2){th= th-512;}
  else if(quad==3){th=  1024-th;}

temp1= 0;
temp2= 128;     //2 multiple
m1=0;
m2=Amp;

    temp3=(m1+m2)>>1;
    Amp=temp3;
      for(int i=0;i<accuracy;i++)
        { test=(temp1+temp2)>>1;
          temp3=temp3>>1; 
          if(th>isin_data[test]){temp1=test;  Amp=Amp+temp3; m1=Amp;}
          else if(th<isin_data[test]){temp2=test; Amp=Amp-temp3;  m2=Amp;}
        }

         if(quad==2){Amp= 0-Amp;}
    else if(quad==3){Amp=  0-Amp;}
return(Amp);
}

int fast_cosine(int Amp, int th)
  {  
  th=256-th;  //cos th = sin (90-th) formula
  return(fast_sine(Amp,th)); 
  }

//--------------------------------------------------------------------------------//


//--------------------------------Fast  RSS----------------------------------------//
int fastRSS(int a, int b)
{  if(a==0 && b==0){return(0);}
  int min,max,temp1,temp2;
  byte clevel;
  if(a<0){a=-a;}
  if(b<0){b=-b;}
clevel=0;
if(a>b){max=a;min=b;} else{max=b;min=a;}

  if(max>(min+min+min))
    {return max;}
  else
    {
     temp1=min>>3;  if(temp1==0){temp1=1;}
     temp2=min;
     while(temp2<max){temp2=temp2+temp1;clevel=clevel+1;}
     temp2=RSSdata[clevel];temp1=temp1>>1;  
     for(int i=0;i<temp2;i++){max=max+temp1;}
     return(max);
    }
}
//--------------------------------------------------------------------------------//

Comments

Only logged in users can leave comments

abhilashpatel121

0 Followers

•

0 Projects

Intro