Zeheng Chen / FFT

fft library for mbed

Dependencies:   CMSIS_DSP_401

Dependents:   4180_Tuner mbed_capstone 4180_EditThis_copy 4180_EditThis_copy_Demo_Test

Fork of FFT by Ale C.-

FFT.cpp@1:e7e724e172dc, 2015-04-17 (annotated)

Committer:
zchen78
Date:
Fri Apr 17 03:39:52 2015 +0000
Revision:
1:e7e724e172dc
Parent:
0:e3af07c00c13 
initial commit of fft

Who changed what in which revision?

UserRevisionLine numberNew 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 }