Mladen Adamovic / UIT_FFT_Real

Final version.

Dependents:   Spectrogram

Fork of UIT_FFT_Real by 不韋 呂

fftReal.cpp@3:7e761912ef10, 2017-08-16 (annotated)

Committer:
mladjo1993
Date:
Wed Aug 16 10:06:13 2017 +0000
Revision:
3:7e761912ef10
Parent:
2:9649d0e2bb4a
Child:
4:2e5de15ed0c2
Comments changed to english.

Who changed what in which revision?

UserRevisionLine numberNew contents of line
mladjo1993 3:7e761912ef10 1 /********************************************************************
mladjo1993 3:7e761912ef10 2 * FFT class for real data usind decimation-in-frequency algorithm
mladjo1993 3:7e761912ef10 3 * This class can execute FFT and IFFT
mladjo1993 3:7e761912ef10 4 *
mladjo1993 3:7e761912ef10 5 * Mladen Adamovic, 3326/2016
mladjo1993 3:7e761912ef10 6 ********************************************************************/
MikamiUitOpen 0:982a9acf3a07 7
MikamiUitOpen 0:982a9acf3a07 8 #include "fftReal.hpp"
MikamiUitOpen 0:982a9acf3a07 9
MikamiUitOpen 0:982a9acf3a07 10 namespace Mikami
MikamiUitOpen 0:982a9acf3a07 11 {
MikamiUitOpen 0:982a9acf3a07 12 // Constructor
MikamiUitOpen 0:982a9acf3a07 13 FftReal::FftReal(int16_t n)
MikamiUitOpen 0:982a9acf3a07 14 : N_FFT_(n), N_INV_(1.0f/n)
MikamiUitOpen 0:982a9acf3a07 15 {
MikamiUitOpen 0:982a9acf3a07 16 // __clz(): Count leading zeros
MikamiUitOpen 0:982a9acf3a07 17 uint32_t shifted = n << (__clz(n)+1);
MikamiUitOpen 0:982a9acf3a07 18 if (shifted != 0)
MikamiUitOpen 0:982a9acf3a07 19 {
MikamiUitOpen 0:982a9acf3a07 20 fprintf(stderr, "\r\nNot power of 2, in FftReal class.");
MikamiUitOpen 0:982a9acf3a07 21 fprintf(stderr, "\r\nForce to exit the program.");
MikamiUitOpen 0:982a9acf3a07 22 exit(EXIT_FAILURE); // Terminate program
MikamiUitOpen 0:982a9acf3a07 23 }
MikamiUitOpen 0:982a9acf3a07 24
MikamiUitOpen 0:982a9acf3a07 25 wTable_ = new Complex[n/2];
MikamiUitOpen 0:982a9acf3a07 26 bTable_ = new uint16_t[n];
MikamiUitOpen 0:982a9acf3a07 27 u_ = new Complex[n];
MikamiUitOpen 0:982a9acf3a07 28
MikamiUitOpen 0:982a9acf3a07 29 // calculation of twiddle factor
MikamiUitOpen 0:982a9acf3a07 30 Complex arg = Complex(0, -6.283185f/N_FFT_);
MikamiUitOpen 0:982a9acf3a07 31 for (int k=0; k<N_FFT_/2; k++)
MikamiUitOpen 0:982a9acf3a07 32 wTable_[k] = exp(arg*(float)k);
MikamiUitOpen 0:982a9acf3a07 33
MikamiUitOpen 0:982a9acf3a07 34 // for bit reversal table
MikamiUitOpen 0:982a9acf3a07 35 uint16_t nShift = __clz(n) + 1;
MikamiUitOpen 0:982a9acf3a07 36 for (int k=0; k<n; k++)
MikamiUitOpen 0:982a9acf3a07 37 // __rbit(k): Reverse the bit order in a 32-bit word
MikamiUitOpen 0:982a9acf3a07 38 bTable_[k] = __rbit(k) >> nShift;
MikamiUitOpen 0:982a9acf3a07 39 }
MikamiUitOpen 0:982a9acf3a07 40
MikamiUitOpen 0:982a9acf3a07 41 // Destructor
MikamiUitOpen 0:982a9acf3a07 42 FftReal::~FftReal()
MikamiUitOpen 0:982a9acf3a07 43 {
MikamiUitOpen 0:982a9acf3a07 44 delete[] wTable_;
MikamiUitOpen 0:982a9acf3a07 45 delete[] bTable_;
MikamiUitOpen 0:982a9acf3a07 46 delete[] u_;
MikamiUitOpen 0:982a9acf3a07 47 }
MikamiUitOpen 0:982a9acf3a07 48
MikamiUitOpen 0:982a9acf3a07 49 // Execute FFT
MikamiUitOpen 0:982a9acf3a07 50 void FftReal::Execute(const float x[], Complex y[])
MikamiUitOpen 0:982a9acf3a07 51 {
MikamiUitOpen 0:982a9acf3a07 52 for (int n=0; n<N_FFT_; n++) u_[n] = x[n];
MikamiUitOpen 0:982a9acf3a07 53
MikamiUitOpen 0:982a9acf3a07 54 // except for last stage
MikamiUitOpen 1:559a63853f3f 55 ExcludeLastStage();
MikamiUitOpen 0:982a9acf3a07 56
MikamiUitOpen 0:982a9acf3a07 57 // Last stage
MikamiUitOpen 0:982a9acf3a07 58 y[0] = u_[0] + u_[1];
MikamiUitOpen 0:982a9acf3a07 59 y[N_FFT_/2] = u_[0] - u_[1];
MikamiUitOpen 0:982a9acf3a07 60 for (int k=2; k<N_FFT_; k+=2)
MikamiUitOpen 0:982a9acf3a07 61 u_[k] = u_[k] + u_[k+1];
MikamiUitOpen 0:982a9acf3a07 62
MikamiUitOpen 0:982a9acf3a07 63 // Reorder to bit reversal
MikamiUitOpen 0:982a9acf3a07 64 for (int k=1; k<N_FFT_/2; k++)
MikamiUitOpen 0:982a9acf3a07 65 y[k] = u_[bTable_[k]];
MikamiUitOpen 0:982a9acf3a07 66 }
MikamiUitOpen 0:982a9acf3a07 67
MikamiUitOpen 0:982a9acf3a07 68 // Execute IFFT
MikamiUitOpen 0:982a9acf3a07 69 void FftReal::ExecuteIfft(const Complex y[], float x[])
MikamiUitOpen 0:982a9acf3a07 70 {
MikamiUitOpen 0:982a9acf3a07 71 int half = N_FFT_/2;
MikamiUitOpen 0:982a9acf3a07 72
MikamiUitOpen 0:982a9acf3a07 73 for (int n=0; n<=half; n++) u_[n] = y[n];
MikamiUitOpen 0:982a9acf3a07 74 for (int n=half+1; n<N_FFT_; n++)
MikamiUitOpen 0:982a9acf3a07 75 u_[n] = conj(y[N_FFT_-n]);
MikamiUitOpen 0:982a9acf3a07 76
MikamiUitOpen 0:982a9acf3a07 77 // except for last stage
MikamiUitOpen 1:559a63853f3f 78 ExcludeLastStage();
MikamiUitOpen 0:982a9acf3a07 79
MikamiUitOpen 0:982a9acf3a07 80 // Last stage including bit reversal
MikamiUitOpen 0:982a9acf3a07 81 x[0] = N_INV_*(u_[0].real() + u_[1].real());
MikamiUitOpen 0:982a9acf3a07 82 x[half] = N_INV_*(u_[0].real() - u_[1].real());
MikamiUitOpen 0:982a9acf3a07 83
MikamiUitOpen 0:982a9acf3a07 84 for (int n=2; n<N_FFT_; n+=2)
MikamiUitOpen 0:982a9acf3a07 85 {
MikamiUitOpen 0:982a9acf3a07 86 float un = u_[n].real();
MikamiUitOpen 0:982a9acf3a07 87 float un1 = u_[n+1].real();
MikamiUitOpen 0:982a9acf3a07 88 x[Index(n)] = N_INV_*(un + un1);
MikamiUitOpen 0:982a9acf3a07 89 x[Index(n+1)] = N_INV_*(un - un1);
MikamiUitOpen 0:982a9acf3a07 90 }
MikamiUitOpen 0:982a9acf3a07 91 }
MikamiUitOpen 0:982a9acf3a07 92
MikamiUitOpen 0:982a9acf3a07 93 // Processing except for last stage
MikamiUitOpen 1:559a63853f3f 94 void FftReal::ExcludeLastStage()
MikamiUitOpen 0:982a9acf3a07 95 {
MikamiUitOpen 0:982a9acf3a07 96 uint16_t nHalf = N_FFT_/2;
MikamiUitOpen 0:982a9acf3a07 97 for (int stg=1; stg<N_FFT_/2; stg*=2)
MikamiUitOpen 0:982a9acf3a07 98 {
MikamiUitOpen 0:982a9acf3a07 99 uint16_t nHalf2 = nHalf*2;
MikamiUitOpen 0:982a9acf3a07 100 for (int kp=0; kp<N_FFT_; kp+=nHalf2)
MikamiUitOpen 0:982a9acf3a07 101 {
MikamiUitOpen 0:982a9acf3a07 102 uint16_t kx = 0;
MikamiUitOpen 0:982a9acf3a07 103 for (int k=kp; k<kp+nHalf; k++)
MikamiUitOpen 0:982a9acf3a07 104 {
MikamiUitOpen 0:982a9acf3a07 105 // Butterfly operation
MikamiUitOpen 0:982a9acf3a07 106 Complex uTmp = u_[k+nHalf];
MikamiUitOpen 0:982a9acf3a07 107 u_[k+nHalf] = (u_[k] - uTmp)*wTable_[kx];
MikamiUitOpen 0:982a9acf3a07 108 u_[k] = u_[k] + uTmp;
MikamiUitOpen 0:982a9acf3a07 109 kx = kx + stg;
MikamiUitOpen 0:982a9acf3a07 110 }
MikamiUitOpen 0:982a9acf3a07 111 }
MikamiUitOpen 0:982a9acf3a07 112 nHalf = nHalf/2;
MikamiUitOpen 0:982a9acf3a07 113 }
MikamiUitOpen 0:982a9acf3a07 114 }
MikamiUitOpen 0:982a9acf3a07 115 }
MikamiUitOpen 0:982a9acf3a07 116