Takehisa Oneta / UIT_FFT_Real_mod

for LPC1114FN28 (Cortex-M0)

Fork of UIT_FFT_Real by 不韋 呂

fftReal.cpp@1:ba9ce95ec9a4, 2015-10-15 (annotated)

Committer:
ohneta
Date:
Thu Oct 15 15:37:24 2015 +0000
Revision:
1:ba9ce95ec9a4
Parent:
0:982a9acf3a07 
for Cortex-M0

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 0:982a9acf3a07 4 // Copyright (c) 2014 MIKAMI, Naoki, 2014/12/19
MikamiUitOpen 0:982a9acf3a07 5 //------------------------------------------------------------------------------
ohneta 1:ba9ce95ec9a4 6 // for LPPC1114FN28(Cortex-M0) modified by Takehisa Oneta(ohneta) 4th Aug.2015
MikamiUitOpen 0:982a9acf3a07 7
MikamiUitOpen 0:982a9acf3a07 8 #include "fftReal.hpp"
ohneta 1:ba9ce95ec9a4 9 #include "mbed.h"
MikamiUitOpen 0:982a9acf3a07 10
ohneta 1:ba9ce95ec9a4 11 extern Serial pc;
ohneta 1:ba9ce95ec9a4 12
ohneta 1:ba9ce95ec9a4 13 // Constructor
ohneta 1:ba9ce95ec9a4 14 FftReal::FftReal(int16_t n) : N_FFT_(n), N_INV_(1.0f / n)
MikamiUitOpen 0:982a9acf3a07 15 {
ohneta 1:ba9ce95ec9a4 16 // __clz(): Count leading zeros
ohneta 1:ba9ce95ec9a4 17 uint32_t shifted = n << (__clz(n) + 1);
ohneta 1:ba9ce95ec9a4 18 if (shifted != 0) {
ohneta 1:ba9ce95ec9a4 19 pc.printf("\r\nNot power of 2, in FftReal class.");
ohneta 1:ba9ce95ec9a4 20 pc.printf("\r\nForce to exit the program.");
ohneta 1:ba9ce95ec9a4 21 exit(EXIT_FAILURE); // Terminate program
ohneta 1:ba9ce95ec9a4 22 }
MikamiUitOpen 0:982a9acf3a07 23
ohneta 1:ba9ce95ec9a4 24 wTable_ = new Complex[n / 2];
ohneta 1:ba9ce95ec9a4 25 if (wTable_ == NULL) {
ohneta 1:ba9ce95ec9a4 26 pc.printf("wTable_ is NULL.\n");
ohneta 1:ba9ce95ec9a4 27 exit(EXIT_FAILURE);
ohneta 1:ba9ce95ec9a4 28 }
ohneta 1:ba9ce95ec9a4 29 bTable_ = new uint16_t[n];
ohneta 1:ba9ce95ec9a4 30 if (bTable_ == NULL) {
ohneta 1:ba9ce95ec9a4 31 pc.printf("bTable_ is NULL.\n");
ohneta 1:ba9ce95ec9a4 32 exit(EXIT_FAILURE);
ohneta 1:ba9ce95ec9a4 33 }
ohneta 1:ba9ce95ec9a4 34 u_ = new Complex[n];
ohneta 1:ba9ce95ec9a4 35 if (u_ == NULL) {
ohneta 1:ba9ce95ec9a4 36 pc.printf("u_ is NULL.\n");
ohneta 1:ba9ce95ec9a4 37 exit(EXIT_FAILURE);
ohneta 1:ba9ce95ec9a4 38 }
MikamiUitOpen 0:982a9acf3a07 39
ohneta 1:ba9ce95ec9a4 40 // calculation of twiddle factor
ohneta 1:ba9ce95ec9a4 41 Complex arg = Complex(0, -6.283185f / N_FFT_);
ohneta 1:ba9ce95ec9a4 42 for (int k = 0; k < N_FFT_ / 2; k++) {
ohneta 1:ba9ce95ec9a4 43 wTable_[k] = exp(arg * (float)k);
MikamiUitOpen 0:982a9acf3a07 44 }
MikamiUitOpen 0:982a9acf3a07 45
ohneta 1:ba9ce95ec9a4 46 // for bit reversal table
ohneta 1:ba9ce95ec9a4 47 uint16_t nShift = __clz(n) + 1;
ohneta 1:ba9ce95ec9a4 48 for (int k = 0; k < n; k++) {
ohneta 1:ba9ce95ec9a4 49 //bTable_[k] = __rbit(k) >> nShift; // __rbit(k): Reverse the bit order in a 32-bit word
ohneta 1:ba9ce95ec9a4 50 unsigned int x = (unsigned int)k;
ohneta 1:ba9ce95ec9a4 51 x = ((x & 0x55555555) << 1) | ((x & 0xAAAAAAAA) >> 1);
ohneta 1:ba9ce95ec9a4 52 x = ((x & 0x33333333) << 2) | ((x & 0xCCCCCCCC) >> 2);
ohneta 1:ba9ce95ec9a4 53 x = ((x & 0x0F0F0F0F) << 4) | ((x & 0xF0F0F0F0) >> 4);
ohneta 1:ba9ce95ec9a4 54 x = ((x & 0x00FF00FF) << 8) | ((x & 0xFF00FF00) >> 8);
ohneta 1:ba9ce95ec9a4 55 x = ((x & 0x0000FFFF) << 16) | ((x & 0xFFFF0000) >> 16);
ohneta 1:ba9ce95ec9a4 56 bTable_[k] = (x >> nShift);
ohneta 1:ba9ce95ec9a4 57 }
ohneta 1:ba9ce95ec9a4 58 }
ohneta 1:ba9ce95ec9a4 59
ohneta 1:ba9ce95ec9a4 60 // Destructor
ohneta 1:ba9ce95ec9a4 61 FftReal::~FftReal()
ohneta 1:ba9ce95ec9a4 62 {
ohneta 1:ba9ce95ec9a4 63 delete[] wTable_;
ohneta 1:ba9ce95ec9a4 64 delete[] bTable_;
ohneta 1:ba9ce95ec9a4 65 delete[] u_;
ohneta 1:ba9ce95ec9a4 66 }
ohneta 1:ba9ce95ec9a4 67
ohneta 1:ba9ce95ec9a4 68 // Execute FFT
ohneta 1:ba9ce95ec9a4 69 void FftReal::Execute(const float x[], Complex y[])
ohneta 1:ba9ce95ec9a4 70 {
ohneta 1:ba9ce95ec9a4 71 for (int n = 0; n < N_FFT_; n++) {
ohneta 1:ba9ce95ec9a4 72 u_[n] = x[n];
MikamiUitOpen 0:982a9acf3a07 73 }
MikamiUitOpen 0:982a9acf3a07 74
ohneta 1:ba9ce95ec9a4 75 // except for last stage
ohneta 1:ba9ce95ec9a4 76 ExcludeLastTtage();
MikamiUitOpen 0:982a9acf3a07 77
ohneta 1:ba9ce95ec9a4 78 // Last stage
ohneta 1:ba9ce95ec9a4 79 y[0] = u_[0] + u_[1];
ohneta 1:ba9ce95ec9a4 80 y[N_FFT_ / 2] = u_[0] - u_[1];
ohneta 1:ba9ce95ec9a4 81 for (int k = 2; k < N_FFT_; k += 2) {
ohneta 1:ba9ce95ec9a4 82 u_[k] = u_[k] + u_[k + 1];
MikamiUitOpen 0:982a9acf3a07 83 }
MikamiUitOpen 0:982a9acf3a07 84
ohneta 1:ba9ce95ec9a4 85 // Reorder to bit reversal
ohneta 1:ba9ce95ec9a4 86 for (int k = 1; k < N_FFT_ / 2; k++) {
ohneta 1:ba9ce95ec9a4 87 y[k] = u_[bTable_[k]];
ohneta 1:ba9ce95ec9a4 88 }
ohneta 1:ba9ce95ec9a4 89 }
ohneta 1:ba9ce95ec9a4 90
ohneta 1:ba9ce95ec9a4 91 #if 0
ohneta 1:ba9ce95ec9a4 92 // Execute IFFT
ohneta 1:ba9ce95ec9a4 93 void FftReal::ExecuteIfft(const Complex y[], float x[])
MikamiUitOpen 0:982a9acf3a07 94 {
MikamiUitOpen 0:982a9acf3a07 95 int half = N_FFT_/2;
MikamiUitOpen 0:982a9acf3a07 96
ohneta 1:ba9ce95ec9a4 97 for (int n=0; n<=half; n++) {
ohneta 1:ba9ce95ec9a4 98 u_[n] = y[n];
ohneta 1:ba9ce95ec9a4 99 }
ohneta 1:ba9ce95ec9a4 100 for (int n=half+1; n<N_FFT_; n++) {
MikamiUitOpen 0:982a9acf3a07 101 u_[n] = conj(y[N_FFT_-n]);
ohneta 1:ba9ce95ec9a4 102 }
MikamiUitOpen 0:982a9acf3a07 103
MikamiUitOpen 0:982a9acf3a07 104 // except for last stage
MikamiUitOpen 0:982a9acf3a07 105 ExcludeLastTtage();
MikamiUitOpen 0:982a9acf3a07 106
MikamiUitOpen 0:982a9acf3a07 107 // Last stage including bit reversal
MikamiUitOpen 0:982a9acf3a07 108 x[0] = N_INV_*(u_[0].real() + u_[1].real());
MikamiUitOpen 0:982a9acf3a07 109 x[half] = N_INV_*(u_[0].real() - u_[1].real());
MikamiUitOpen 0:982a9acf3a07 110
MikamiUitOpen 0:982a9acf3a07 111 for (int n=2; n<N_FFT_; n+=2)
MikamiUitOpen 0:982a9acf3a07 112 {
MikamiUitOpen 0:982a9acf3a07 113 float un = u_[n].real();
MikamiUitOpen 0:982a9acf3a07 114 float un1 = u_[n+1].real();
MikamiUitOpen 0:982a9acf3a07 115 x[Index(n)] = N_INV_*(un + un1);
MikamiUitOpen 0:982a9acf3a07 116 x[Index(n+1)] = N_INV_*(un - un1);
MikamiUitOpen 0:982a9acf3a07 117 }
ohneta 1:ba9ce95ec9a4 118 }
ohneta 1:ba9ce95ec9a4 119 #endif
MikamiUitOpen 0:982a9acf3a07 120
ohneta 1:ba9ce95ec9a4 121 // Processing except for last stage
ohneta 1:ba9ce95ec9a4 122 void FftReal::ExcludeLastTtage()
ohneta 1:ba9ce95ec9a4 123 {
ohneta 1:ba9ce95ec9a4 124 uint16_t nHalf = N_FFT_ / 2;
ohneta 1:ba9ce95ec9a4 125 for (int stg = 1; stg < N_FFT_ / 2; stg *= 2) {
ohneta 1:ba9ce95ec9a4 126 uint16_t nHalf2 = nHalf * 2;
ohneta 1:ba9ce95ec9a4 127 for (int kp = 0; kp < N_FFT_; kp += nHalf2) {
ohneta 1:ba9ce95ec9a4 128 uint16_t kx = 0;
ohneta 1:ba9ce95ec9a4 129 for (int k = kp; k <kp + nHalf; k++) {
ohneta 1:ba9ce95ec9a4 130 // Butterfly operation
ohneta 1:ba9ce95ec9a4 131 Complex uTmp = u_[k + nHalf];
ohneta 1:ba9ce95ec9a4 132 u_[k + nHalf] = (u_[k] - uTmp) * wTable_[kx];
ohneta 1:ba9ce95ec9a4 133 u_[k] = u_[k] + uTmp;
ohneta 1:ba9ce95ec9a4 134 kx = kx + stg;
MikamiUitOpen 0:982a9acf3a07 135 }
ohneta 1:ba9ce95ec9a4 136 }
ohneta 1:ba9ce95ec9a4 137 nHalf = nHalf / 2;
ohneta 1:ba9ce95ec9a4 138 }
MikamiUitOpen 0:982a9acf3a07 139 }
MikamiUitOpen 0:982a9acf3a07 140
ohneta 1:ba9ce95ec9a4 141