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//--------------------------------------------------------------------------------//
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//--------------------------------------------------------------------------------//
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//--------------------------------------------------------------------------------//
Comments
Only logged in users can leave comments
gusto2
2 months ago
Thank you very much for the inspiration.
Anonymous user
2 years ago
interesting, i wonder how much could this or EasyFFT be expanded with Teensy 4.1
Anonymous user
2 years ago
This is great! One note though, I think there is an error on line 70 of ApproxFFT.ino and ApproxFFT above (but fixed in ApproxFFT V2). I think: `out_r[i]=0; out_r[i]=0;` should be: `out_r[i]=0; out_im[i]=0;`. In the current code, out_im[] never gets initialized, and it hangs without throwing any warnings. Otherwise, awesome code! Thanks for sharing!
Anonymous user
2 years ago
I like this efficient code! It would be good to not use native data types, so the type size is the same on all boards, so don't use int, but int16_t and so on..
Anonymous user
3 years ago
interesting, i wonder how much could this or EasyFFT be expanded with Teensy 4.1
Anonymous user
4 years ago
Nice. Does not use machine code and seems easily portable.
abhilashpatel121
6 Followers
•9 Projects
1
7
ApproxFFT: Fastest FFT Function for Arduino | Arduino Project Hub
carloneb
a month ago
Thanks for your work, I'm building a mini oscilloscope with stm32 (bluepill) and it seems to work well, now I'm implementing the FFT function with your last ApproxFFT_V2_060521. Everything seems to work with 1024 samples. With a sinusoid (at various test frequencies) I can correctly plot the fundamental frequency and the relative maximum peak. What I don't understand is that with a square wave input, in addition to the fundamental, I don't see all the "infinite" harmonics, what you see is similar to the plot of a sinusoidal wave. Am I doing things wrong? Do you have any suggestions? Thanks