fft library for mbed
Dependents: 4180_Tuner mbed_capstone 4180_EditThis_copy 4180_EditThis_copy_Demo_Test
Fork of FFT by
FFT.cpp@0:e3af07c00c13, 2011-02-11 (annotated)
- Committer:
- Suky
- Date:
- Fri Feb 11 00:35:40 2011 +0000
- Revision:
- 0:e3af07c00c13
Who changed what in which revision?
User | Revision | Line number | New contents of line |
---|---|---|---|
Suky | 0:e3af07c00c13 | 1 | /* |
Suky | 0:e3af07c00c13 | 2 | @file FFT.cpp |
Suky | 0:e3af07c00c13 | 3 | @version: 1.0 |
Suky | 0:e3af07c00c13 | 4 | @author: Suky |
Suky | 0:e3af07c00c13 | 5 | @web www.micros-designs.com.ar |
Suky | 0:e3af07c00c13 | 6 | @date 10/02/11 |
Suky | 0:e3af07c00c13 | 7 | */ |
Suky | 0:e3af07c00c13 | 8 | #include "FFT.h" |
Suky | 0:e3af07c00c13 | 9 | |
Suky | 0:e3af07c00c13 | 10 | // Extracted from Numerical Recipes in C |
Suky | 0:e3af07c00c13 | 11 | void vFFT(float data[], unsigned int nn){ |
Suky | 0:e3af07c00c13 | 12 | /*Replaces data[1..2*nn] by its discrete Fourier transform, if isign is input as 1; or replaces |
Suky | 0:e3af07c00c13 | 13 | data[1..2*nn] by nn times its inverse discrete Fourier transform, if isign is input as -1. |
Suky | 0:e3af07c00c13 | 14 | data is a complex array of length nn or, equivalently, a real array of length 2*nn. nn MUST |
Suky | 0:e3af07c00c13 | 15 | be an integer power of 2 (this is not checked for!).*/ |
Suky | 0:e3af07c00c13 | 16 | unsigned int n,mmax,m,j,istep,i; |
Suky | 0:e3af07c00c13 | 17 | double wtemp,wr,wpr,wpi,wi,theta; |
Suky | 0:e3af07c00c13 | 18 | float tempr,tempi; |
Suky | 0:e3af07c00c13 | 19 | |
Suky | 0:e3af07c00c13 | 20 | #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr |
Suky | 0:e3af07c00c13 | 21 | |
Suky | 0:e3af07c00c13 | 22 | n=nn << 1; |
Suky | 0:e3af07c00c13 | 23 | j=1; |
Suky | 0:e3af07c00c13 | 24 | for (i=1;i<n;i+=2) { |
Suky | 0:e3af07c00c13 | 25 | if(j>i){ |
Suky | 0:e3af07c00c13 | 26 | SWAP(data[j],data[i]); |
Suky | 0:e3af07c00c13 | 27 | SWAP(data[j+1],data[i+1]); |
Suky | 0:e3af07c00c13 | 28 | } |
Suky | 0:e3af07c00c13 | 29 | m=n >> 1; |
Suky | 0:e3af07c00c13 | 30 | while (m >= 2 &&j>m){ |
Suky | 0:e3af07c00c13 | 31 | j-=m; |
Suky | 0:e3af07c00c13 | 32 | m >>= 1; |
Suky | 0:e3af07c00c13 | 33 | } |
Suky | 0:e3af07c00c13 | 34 | j+=m; |
Suky | 0:e3af07c00c13 | 35 | } |
Suky | 0:e3af07c00c13 | 36 | |
Suky | 0:e3af07c00c13 | 37 | mmax=2; |
Suky | 0:e3af07c00c13 | 38 | while (n > mmax) { |
Suky | 0:e3af07c00c13 | 39 | istep=mmax << 1; |
Suky | 0:e3af07c00c13 | 40 | theta=(6.28318530717959/mmax); |
Suky | 0:e3af07c00c13 | 41 | wtemp=sin(0.5*theta); |
Suky | 0:e3af07c00c13 | 42 | wpr = -2.0*wtemp*wtemp; |
Suky | 0:e3af07c00c13 | 43 | wpi=sin(theta); |
Suky | 0:e3af07c00c13 | 44 | wr=1.0; |
Suky | 0:e3af07c00c13 | 45 | wi=0.0; |
Suky | 0:e3af07c00c13 | 46 | for (m=1;m<mmax;m+=2) { |
Suky | 0:e3af07c00c13 | 47 | for (i=m;i<=n;i+=istep) { |
Suky | 0:e3af07c00c13 | 48 | j=i+mmax; |
Suky | 0:e3af07c00c13 | 49 | tempr=wr*data[j]-wi*data[j+1]; |
Suky | 0:e3af07c00c13 | 50 | tempi=wr*data[j+1]+wi*data[j]; |
Suky | 0:e3af07c00c13 | 51 | data[j]=data[i]-tempr; |
Suky | 0:e3af07c00c13 | 52 | data[j+1]=data[i+1]-tempi; |
Suky | 0:e3af07c00c13 | 53 | data[i] += tempr; |
Suky | 0:e3af07c00c13 | 54 | data[i+1] += tempi; |
Suky | 0:e3af07c00c13 | 55 | } |
Suky | 0:e3af07c00c13 | 56 | wr=(wtemp=wr)*wpr-wi*wpi+wr; |
Suky | 0:e3af07c00c13 | 57 | wi=wi*wpr+wtemp*wpi+wi; |
Suky | 0:e3af07c00c13 | 58 | } |
Suky | 0:e3af07c00c13 | 59 | mmax=istep; |
Suky | 0:e3af07c00c13 | 60 | } |
Suky | 0:e3af07c00c13 | 61 | } |
Suky | 0:e3af07c00c13 | 62 | |
Suky | 0:e3af07c00c13 | 63 | // Extracted from Numerical Recipes in C |
Suky | 0:e3af07c00c13 | 64 | void vRealFFT(float data[], unsigned int n){ |
Suky | 0:e3af07c00c13 | 65 | /*Calculates the Fourier transform of a set of n real-valued data points. Replaces this data (which |
Suky | 0:e3af07c00c13 | 66 | is stored in array data[1..n]) by the positive frequency half of its complex Fourier transform. |
Suky | 0:e3af07c00c13 | 67 | The real-valued rst and last components of the complex transform are returned as elements |
Suky | 0:e3af07c00c13 | 68 | data[1] and data[2], respectively. n must be a power of 2. This routine also calculates the |
Suky | 0:e3af07c00c13 | 69 | inverse transform of a complex data array if it is the transform of real data. (Result in this case |
Suky | 0:e3af07c00c13 | 70 | must be multiplied by 2/n.)*/ |
Suky | 0:e3af07c00c13 | 71 | unsigned long i,i1,i2,i3,i4,np3; |
Suky | 0:e3af07c00c13 | 72 | float c1=0.5,c2,h1r,h1i,h2r,h2i; |
Suky | 0:e3af07c00c13 | 73 | double wr,wi,wpr,wpi,wtemp,theta; |
Suky | 0:e3af07c00c13 | 74 | theta=3.141592653589793/(double) (n>>1); |
Suky | 0:e3af07c00c13 | 75 | |
Suky | 0:e3af07c00c13 | 76 | c2 = -0.5; |
Suky | 0:e3af07c00c13 | 77 | vFFT(data,n>>1); |
Suky | 0:e3af07c00c13 | 78 | wtemp=sin(0.5*theta); |
Suky | 0:e3af07c00c13 | 79 | wpr = -2.0*wtemp*wtemp; |
Suky | 0:e3af07c00c13 | 80 | wpi=sin(theta); |
Suky | 0:e3af07c00c13 | 81 | wr=1.0+wpr; |
Suky | 0:e3af07c00c13 | 82 | wi=wpi; |
Suky | 0:e3af07c00c13 | 83 | np3=n+3; |
Suky | 0:e3af07c00c13 | 84 | for (i=2;i<=(n>>2);i++) { |
Suky | 0:e3af07c00c13 | 85 | i4=1+(i3=np3-(i2=1+(i1=i+i-1))); |
Suky | 0:e3af07c00c13 | 86 | h1r=c1*(data[i1]+data[i3]); |
Suky | 0:e3af07c00c13 | 87 | h1i=c1*(data[i2]-data[i4]); |
Suky | 0:e3af07c00c13 | 88 | h2r = -c2*(data[i2]+data[i4]); |
Suky | 0:e3af07c00c13 | 89 | h2i=c2*(data[i1]-data[i3]); |
Suky | 0:e3af07c00c13 | 90 | data[i1]=h1r+wr*h2r-wi*h2i; |
Suky | 0:e3af07c00c13 | 91 | data[i2]=h1i+wr*h2i+wi*h2r; |
Suky | 0:e3af07c00c13 | 92 | data[i3]=h1r-wr*h2r+wi*h2i; |
Suky | 0:e3af07c00c13 | 93 | data[i4] = -h1i+wr*h2i+wi*h2r; |
Suky | 0:e3af07c00c13 | 94 | wr=(wtemp=wr)*wpr-wi*wpi+wr; |
Suky | 0:e3af07c00c13 | 95 | wi=wi*wpr+wtemp*wpi+wi; |
Suky | 0:e3af07c00c13 | 96 | } |
Suky | 0:e3af07c00c13 | 97 | data[1] = (h1r=data[1])+data[2]; |
Suky | 0:e3af07c00c13 | 98 | data[2] = h1r-data[2]; |
Suky | 0:e3af07c00c13 | 99 | |
Suky | 0:e3af07c00c13 | 100 | } |
Suky | 0:e3af07c00c13 | 101 | |
Suky | 0:e3af07c00c13 | 102 | |
Suky | 0:e3af07c00c13 | 103 | void vCalPowerf(float Input[],float Power[], unsigned int n){ |
Suky | 0:e3af07c00c13 | 104 | unsigned char k,j; |
Suky | 0:e3af07c00c13 | 105 | |
Suky | 0:e3af07c00c13 | 106 | for(k=0,j=0;k<n;k++,j+=2){ |
Suky | 0:e3af07c00c13 | 107 | Power[k]=sqrt(Input[j]*Input[j]+Input[j+1]*Input[j+1]); |
Suky | 0:e3af07c00c13 | 108 | } |
Suky | 0:e3af07c00c13 | 109 | } |
Suky | 0:e3af07c00c13 | 110 | |
Suky | 0:e3af07c00c13 | 111 | void vCalPowerInt(float Input[],unsigned char Power[], unsigned int n){ |
Suky | 0:e3af07c00c13 | 112 | unsigned char k,j; |
Suky | 0:e3af07c00c13 | 113 | |
Suky | 0:e3af07c00c13 | 114 | for(k=0,j=0;k<n;k++,j+=2){ |
Suky | 0:e3af07c00c13 | 115 | Power[k]=sqrt(Input[j]*Input[j]+Input[j+1]*Input[j+1]); |
Suky | 0:e3af07c00c13 | 116 | } |
Suky | 0:e3af07c00c13 | 117 | } |
Suky | 0:e3af07c00c13 | 118 | |
Suky | 0:e3af07c00c13 | 119 | void vCalPowerLog(float Input[],unsigned char Power[], unsigned int n){ |
Suky | 0:e3af07c00c13 | 120 | unsigned char k,j; |
Suky | 0:e3af07c00c13 | 121 | float Temp; |
Suky | 0:e3af07c00c13 | 122 | |
Suky | 0:e3af07c00c13 | 123 | for(k=0,j=0;k<n;k++,j+=2){ |
Suky | 0:e3af07c00c13 | 124 | if((Input[j]!=0) && (Input[j+1]!=0)){ |
Suky | 0:e3af07c00c13 | 125 | Temp=sqrt(Input[j]*Input[j]+Input[j+1]*Input[j+1]); |
Suky | 0:e3af07c00c13 | 126 | Power[k]=10.0*log10(Temp); |
Suky | 0:e3af07c00c13 | 127 | }else{ |
Suky | 0:e3af07c00c13 | 128 | Power[k]=0; |
Suky | 0:e3af07c00c13 | 129 | } |
Suky | 0:e3af07c00c13 | 130 | } |
Suky | 0:e3af07c00c13 | 131 | |
Suky | 0:e3af07c00c13 | 132 | } |