QuickFFT: High Speed (low accuracy) FFT for Arduino

1/*
2//example data
3int data[256]={10, 42, 47, 19, -20,  -39,  -26,   7,  30, 23, -10,  -42,  -47,  -19,  20,
439, 26, -7, -30,  -23,  10, 42, 47,   19, -19,  -39,  -26,  7,  30, 23, -9, -42,  -47,  -19,
519, 39, 26, -6, -30,  -24,   9,  42, 47, 20, -19,  -39,  -27,  6,  30, 24, -9, -42,  -47,
6-20,  19, 39,   27, -6, -30,  -24,  9,  41, 47, 20, -19,  -39,  -27,  6,  30, 24, -9, -41,
7-47,   -20,  18, 39, 27, -6, -30,  -24,  8,  41, 47, 20, -18,  -39,  -27,  6,  30, 24,   -8,
8-41,  -47,  -21,  18, 39, 27, -5, -30,  -24,  8,  41, 47, 21, -18,  -39,   -27,  5,  30, 24,
9-8, -41,  -47,  -21,  18, 39, 27, -5, -30,  -24,  8,  41,   47, 21, -18,  -39,  -28,  5,  30,
1025, -7, -41,  -47,  -21,  17, 39, 28, -5, -30,   -25,  7,  40, 47, 22, -17,  -39,  -28,  5,
1130, 25, -7, -40,  -47,  -22,  17,   39, 28, -4, -30,  -25,  7,  40, 48, 22, -17,  -39,  -28,
124,  29, 25, -7, -40,   -48,  -22,  17, 39, 28, -4, -29,  -25,  7,  40, 48, 22, -17,  -39,  
13-28,  4,   29, 25, -6, -40,  -48,  -22,  16, 39, 29, -4, -29,  -25,  6,  40, 48, 23, -16,   
14-39,  -29,  4,  29, 25, -6, -40,  -48,  -23,  16, 39, 29, -3, -29,  -26,  6,   39, 48, 23,
15-16,  -39,  -29,  3,  29, 26, -6, -39,  -48,  -23,  16, 39, 29,   -3, -29,  -26,  5,  39, 48,
1623, -15,  -39,  -29,  3,  29, 26, -5, -39,  -48,   -24,  15, 39};
17*/
18
19void setup() 
20        {
21    // Serial.begin(250000);
22      
23        }
24
25        
26void loop() {
27/*
28float f=Q_FFT(data,256,100);
29Serial.println(f);
30delay(1000000);
31*/
32        
33            }
34
35
36
37//-----------------------------FFT Function----------------------------------------------//
38/*
39Code   to perform High speed (5-7 times faster) and low accuracy FFT on arduino,
40This   code compromises accuracy for speed,
41setup:
42
431. in[]     : Data array,   
442. N        : Number of sample (recommended sample size 2,4,8,16,32,64,128...)
453.   Frequency: sampling frequency required as input (Hz)
46
47It will by default return   frequency with max aplitude,
48
49If sample size is not in power of 2 it will   be clipped to lower side of number. 
50i.e, for 150 number of samples, code will   consider first 128 sample, remaining sample  will be omitted.
51For Arduino nano,   FFT of more than 256 sample not possible due to mamory limitation 
52Code by ABHILASH
53Contact:   abhilashpatel121@gmail.com 
54Documentation & deatails: https://www.instructables.com/member/abhilash_patel/instructables/
55*/
56
57
58float   Q_FFT(int in[],int N,float Frequency)
59{ 
60
61unsigned int Pow2[13]={1,2,4,8,16,32,64,128,256,512,1024,2048};   // declaring this as global array will save 1-2 ms of time
62
63
64int a,c1,f,o,x;          
65byte check=0;
66a=N;  
67                                 
68      for(int   i=0;i<12;i++)                 //calculating the levels
69         { if(Pow2[i]<=a){o=i;}   }
70      
71int out_r[Pow2[o]]={};   //real part of transform
72int out_im[Pow2[o]]={};   //imaginory part of transform
73           
74x=0;  
75      for(int b=0;b<o;b++)                      // bit reversal
76         {
77          c1=Pow2[b];
78           f=Pow2[o]/(c1+c1);
79                for(int j=0;j<c1;j++)
80                    {   
81                     x=x+1;
82                     out_im[x]=out_im[j]+f;
83                     }
84         }
85
86 
87      for(int i=0;i<Pow2[o];i++)             // update input array as per bit reverse order
88         {
89          out_r[i]=in[out_im[i]];   
90          out_im[i]=0;
91         }
92
93
94int i10,i11,n1,tr,ti;
95float   e;
96int c,s;
97    for(int i=0;i<o;i++)                                    //fft
98     {
99     i10=Pow2[i];              // overall values of sine/cosine  
100     i11=Pow2[o]/Pow2[i+1];     // loop with similar sine cosine
101     e=360/Pow2[i+1];
102     e=0-e;
103      n1=0;
104
105          for(int j=0;j<i10;j++)
106          {
107            c=e*j;
108   while(c<0){c=c+360;}
109  while(c>360){c=c-360;}
110
111          n1=j;
112          
113           for(int k=0;k<i11;k++)
114                 {
115
116       if(c==0) {   tr=out_r[i10+n1];
117                  ti=out_im[i10+n1];}
118  else if(c==90){   tr= -out_im[i10+n1];
119                  ti=out_r[i10+n1];}
120  else if(c==180){tr=-out_r[i10+n1];
121                   ti=-out_im[i10+n1];}
122  else if(c==270){tr=out_im[i10+n1];
123                   ti=-out_r[i10+n1];}
124  else if(c==360){tr=out_r[i10+n1];
125                   ti=out_im[i10+n1];}
126  else if(c>0  && c<90)   {tr=out_r[i10+n1]-out_im[i10+n1];
127                            ti=out_im[i10+n1]+out_r[i10+n1];}
128  else if(c>90  &&   c<180) {tr=-out_r[i10+n1]-out_im[i10+n1];
129                           ti=-out_im[i10+n1]+out_r[i10+n1];}
130   else if(c>180 && c<270) {tr=-out_r[i10+n1]+out_im[i10+n1];
131                           ti=-out_im[i10+n1]-out_r[i10+n1];}
132   else if(c>270 && c<360) {tr=out_r[i10+n1]+out_im[i10+n1];
133                           ti=out_im[i10+n1]-out_r[i10+n1];}
134           
135                 out_r[n1+i10]=out_r[n1]-tr;
136                 out_r[n1]=out_r[n1]+tr;
137                  if(out_r[n1]>15000 || out_r[n1]<-15000){check=1;}
138          
139                  out_im[n1+i10]=out_im[n1]-ti;
140                 out_im[n1]=out_im[n1]+ti;
141                  if(out_im[n1]>15000 || out_im[n1]<-15000){check=1;}          
142           
143                 n1=n1+i10+i10;
144                  }       
145              }
146
147    if(check==1){                                             //   scale the matrics if value higher than 15000 to prevent varible from overloading
148                 for(int i=0;i<Pow2[o];i++)
149                    {
150                     out_r[i]=out_r[i]/100;
151                      out_im[i]=out_im[i]/100;    
152                    }
153                     check=0;   
154                }           
155
156     }
157
158/*
159for(int i=0;i<Pow2[o];i++)
160{
161Serial.print(out_r[i]);
162Serial.print("\	");                                      // un comment to print RAW o/p    
163Serial.print(out_im[i]);   Serial.println("i");      
164}
165*/
166
167//---> here onward out_r contains   amplitude and our_in conntains frequency (Hz)
168int fout,fm,fstp;
169float fstep;
170fstep=Frequency/N;
171fstp=fstep;
172fout=0;fm=0;
173
174     for(int i=1;i<Pow2[o-1];i++)               // getting amplitude from compex   number
175        {
176        if((out_r[i]>=0) && (out_im[i]>=0)){out_r[i]=out_r[i]+out_im[i];}
177    else if((out_r[i]<=0) && (out_im[i]<=0)){out_r[i]=-out_r[i]-out_im[i];}
178   else   if((out_r[i]>=0) && (out_im[i]<=0)){out_r[i]=out_r[i]-out_im[i];}
179   else if((out_r[i]<=0)   && (out_im[i]>=0)){out_r[i]=-out_r[i]+out_im[i];}
180   // to find peak sum of mod   of real and imaginery part are considered to increase speed
181        
182out_im[i]=out_im[i-1]+fstp;
183if   (fout<out_r[i]){fm=i; fout=out_r[i];}
184         /*
185         Serial.print(out_im[i]);Serial.print("Hz");
186          Serial.print("\	");                            // un comment to print   freuency bin    
187         Serial.println(out_r[i]); 
188          */
189        }
190
191
192float   fa,fb,fc;
193fa=out_r[fm-1];
194fb=out_r[fm]; 
195fc=out_r[fm+1];
196fstep=(fa*(fm-1)+fb*fm+fc*(fm+1))/(fa+fb+fc);
197
198return(fstep*Frequency/N);
199}
200     
201
202//------------------------------------------------------------------------------------//

1/*
2//example data
3int data[256]={10, 42, 47, 19, -20,  -39,  -26,  7,  30, 23, -10,  -42,  -47,  -19,  20,
439, 26, -7, -30,  -23,  10, 42, 47,  19, -19,  -39,  -26,  7,  30, 23, -9, -42,  -47,  -19,
519, 39, 26, -6, -30,  -24,  9,  42, 47, 20, -19,  -39,  -27,  6,  30, 24, -9, -42,  -47,
6-20,  19, 39,  27, -6, -30,  -24,  9,  41, 47, 20, -19,  -39,  -27,  6,  30, 24, -9, -41,
7-47,  -20,  18, 39, 27, -6, -30,  -24,  8,  41, 47, 20, -18,  -39,  -27,  6,  30, 24,  -8,
8-41,  -47,  -21,  18, 39, 27, -5, -30,  -24,  8,  41, 47, 21, -18,  -39,  -27,  5,  30, 24,
9-8, -41,  -47,  -21,  18, 39, 27, -5, -30,  -24,  8,  41,  47, 21, -18,  -39,  -28,  5,  30,
1025, -7, -41,  -47,  -21,  17, 39, 28, -5, -30,  -25,  7,  40, 47, 22, -17,  -39,  -28,  5,
1130, 25, -7, -40,  -47,  -22,  17,  39, 28, -4, -30,  -25,  7,  40, 48, 22, -17,  -39,  -28,
124,  29, 25, -7, -40,  -48,  -22,  17, 39, 28, -4, -29,  -25,  7,  40, 48, 22, -17,  -39,  
13-28,  4,  29, 25, -6, -40,  -48,  -22,  16, 39, 29, -4, -29,  -25,  6,  40, 48, 23, -16,  
14-39,  -29,  4,  29, 25, -6, -40,  -48,  -23,  16, 39, 29, -3, -29,  -26,  6,  39, 48, 23,
15-16,  -39,  -29,  3,  29, 26, -6, -39,  -48,  -23,  16, 39, 29,  -3, -29,  -26,  5,  39, 48,
1623, -15,  -39,  -29,  3,  29, 26, -5, -39,  -48,  -24,  15, 39};
17*/
18
19void setup() 
20        {
21    // Serial.begin(250000);
22     
23        }
24
25        
26void loop() {
27/*
28float f=Q_FFT(data,256,100);
29Serial.println(f);
30delay(1000000);
31*/
32       
33            }
34
35
36
37//-----------------------------FFT Function----------------------------------------------//
38/*
39Code  to perform High speed (5-7 times faster) and low accuracy FFT on arduino,
40This  code compromises accuracy for speed,
41setup:
42
431. in[]     : Data array,  
442. N        : Number of sample (recommended sample size 2,4,8,16,32,64,128...)
453.  Frequency: sampling frequency required as input (Hz)
46
47It will by default return  frequency with max aplitude,
48
49If sample size is not in power of 2 it will  be clipped to lower side of number. 
50i.e, for 150 number of samples, code will  consider first 128 sample, remaining sample  will be omitted.
51For Arduino nano,  FFT of more than 256 sample not possible due to mamory limitation 
52Code by ABHILASH
53Contact:  abhilashpatel121@gmail.com 
54Documentation & deatails: https://www.instructables.com/member/abhilash_patel/instructables/
55*/
56
57
58float  Q_FFT(int in[],int N,float Frequency)
59{ 
60
61unsigned int Pow2[13]={1,2,4,8,16,32,64,128,256,512,1024,2048};  // declaring this as global array will save 1-2 ms of time
62
63
64int a,c1,f,o,x;         
65byte check=0;
66a=N;  
67                                 
68      for(int  i=0;i<12;i++)                 //calculating the levels
69         { if(Pow2[i]<=a){o=i;}  }
70      
71int out_r[Pow2[o]]={};   //real part of transform
72int out_im[Pow2[o]]={};  //imaginory part of transform
73           
74x=0;  
75      for(int b=0;b<o;b++)                     // bit reversal
76         {
77          c1=Pow2[b];
78          f=Pow2[o]/(c1+c1);
79                for(int j=0;j<c1;j++)
80                    {  
81                     x=x+1;
82                     out_im[x]=out_im[j]+f;
83                    }
84         }
85
86 
87      for(int i=0;i<Pow2[o];i++)            // update input array as per bit reverse order
88         {
89          out_r[i]=in[out_im[i]];  
90          out_im[i]=0;
91         }
92
93
94int i10,i11,n1,tr,ti;
95float  e;
96int c,s;
97    for(int i=0;i<o;i++)                                    //fft
98    {
99     i10=Pow2[i];              // overall values of sine/cosine  
100     i11=Pow2[o]/Pow2[i+1];    // loop with similar sine cosine
101     e=360/Pow2[i+1];
102     e=0-e;
103     n1=0;
104
105          for(int j=0;j<i10;j++)
106          {
107            c=e*j;
108  while(c<0){c=c+360;}
109  while(c>360){c=c-360;}
110
111          n1=j;
112          
113          for(int k=0;k<i11;k++)
114                 {
115
116       if(c==0) {  tr=out_r[i10+n1];
117                  ti=out_im[i10+n1];}
118  else if(c==90){  tr= -out_im[i10+n1];
119                  ti=out_r[i10+n1];}
120  else if(c==180){tr=-out_r[i10+n1];
121                  ti=-out_im[i10+n1];}
122  else if(c==270){tr=out_im[i10+n1];
123                  ti=-out_r[i10+n1];}
124  else if(c==360){tr=out_r[i10+n1];
125                  ti=out_im[i10+n1];}
126  else if(c>0  && c<90)   {tr=out_r[i10+n1]-out_im[i10+n1];
127                           ti=out_im[i10+n1]+out_r[i10+n1];}
128  else if(c>90  &&  c<180) {tr=-out_r[i10+n1]-out_im[i10+n1];
129                           ti=-out_im[i10+n1]+out_r[i10+n1];}
130  else if(c>180 && c<270) {tr=-out_r[i10+n1]+out_im[i10+n1];
131                           ti=-out_im[i10+n1]-out_r[i10+n1];}
132  else if(c>270 && c<360) {tr=out_r[i10+n1]+out_im[i10+n1];
133                           ti=out_im[i10+n1]-out_r[i10+n1];}
134          
135                 out_r[n1+i10]=out_r[n1]-tr;
136                 out_r[n1]=out_r[n1]+tr;
137                 if(out_r[n1]>15000 || out_r[n1]<-15000){check=1;}
138          
139                 out_im[n1+i10]=out_im[n1]-ti;
140                 out_im[n1]=out_im[n1]+ti;
141                 if(out_im[n1]>15000 || out_im[n1]<-15000){check=1;}          
142          
143                 n1=n1+i10+i10;
144                  }       
145             }
146
147    if(check==1){                                             //  scale the matrics if value higher than 15000 to prevent varible from overloading
148                for(int i=0;i<Pow2[o];i++)
149                    {
150                     out_r[i]=out_r[i]/100;
151                     out_im[i]=out_im[i]/100;    
152                    }
153                     check=0;  
154                }           
155
156     }
157
158/*
159for(int i=0;i<Pow2[o];i++)
160{
161Serial.print(out_r[i]);
162Serial.print("\	");                                     // un comment to print RAW o/p    
163Serial.print(out_im[i]);  Serial.println("i");      
164}
165*/
166
167//---> here onward out_r contains  amplitude and our_in conntains frequency (Hz)
168int fout,fm,fstp;
169float fstep;
170fstep=Frequency/N;
171fstp=fstep;
172fout=0;fm=0;
173
174    for(int i=1;i<Pow2[o-1];i++)               // getting amplitude from compex  number
175        {
176        if((out_r[i]>=0) && (out_im[i]>=0)){out_r[i]=out_r[i]+out_im[i];}
177   else if((out_r[i]<=0) && (out_im[i]<=0)){out_r[i]=-out_r[i]-out_im[i];}
178   else  if((out_r[i]>=0) && (out_im[i]<=0)){out_r[i]=out_r[i]-out_im[i];}
179   else if((out_r[i]<=0)  && (out_im[i]>=0)){out_r[i]=-out_r[i]+out_im[i];}
180   // to find peak sum of mod  of real and imaginery part are considered to increase speed
181        
182out_im[i]=out_im[i-1]+fstp;
183if  (fout<out_r[i]){fm=i; fout=out_r[i];}
184         /*
185         Serial.print(out_im[i]);Serial.print("Hz");
186         Serial.print("\	");                            // un comment to print  freuency bin    
187         Serial.println(out_r[i]); 
188          */
189        }
190
191
192float  fa,fb,fc;
193fa=out_r[fm-1];
194fb=out_r[fm]; 
195fc=out_r[fm+1];
196fstep=(fa*(fm-1)+fb*fm+fc*(fm+1))/(fa+fb+fc);
197
198return(fstep*Frequency/N);
199}
200    
201
202//------------------------------------------------------------------------------------//