Takehisa Oneta / UIT_FFT_Real_mod-ohneta

自分用に書き換えてしまったので forkします

Fork of UIT_FFT_Real by 不韋 呂

fftReal.cpp@3:9649d0e2bb4a, 2015-12-18 (annotated)

Committer:
MikamiUitOpen
Date:
Fri Dec 18 10:01:28 2015 +0000
Revision:
3:9649d0e2bb4a
Parent:
2:559a63853f3f 
3

Who changed what in which revision?

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