26 #define GRT_DLL_EXPORTS
27 #include "FastFourierTransform.h"
31 FastFourierTransform::FastFourierTransform(){
33 computeMagnitude =
true;
35 enableZeroPadding =
true;
37 windowFunction = RECTANGULAR_WINDOW;
40 infoLog.setProceedingText(
"[FastFourierTransform]");
41 warningLog.setProceedingText(
"[WARNING FastFourierTransform]");
42 errorLog.setProceedingText(
"[ERROR FastFourierTransform]");
49 this->initialized = rhs.initialized;
50 this->computeMagnitude = rhs.computeMagnitude;
51 this->computePhase = rhs.computePhase;
52 this->enableZeroPadding = rhs.enableZeroPadding;
53 this->windowSize = rhs.windowSize;
54 this->windowFunction = rhs.windowFunction;
55 this->averagePower = 0;
57 this->infoLog = rhs.infoLog;
58 this->warningLog = rhs.warningLog;
59 this->errorLog = rhs.errorLog;
61 if( rhs.initialized ){
62 this->init(rhs.windowSize,rhs.windowFunction,rhs.computeMagnitude,rhs.computePhase);
65 for(UINT i=0; i<this->windowSize; i++){
66 this->fftReal[i] = rhs.fftReal[i];
67 this->fftImag[i] = rhs.fftImag[i];
68 this->magnitude[i] = rhs.magnitude[i];
69 this->phase[i] = rhs.phase[i];
70 this->power[i] = rhs.power[i];
75 FastFourierTransform::~FastFourierTransform() {
82 this->initialized = rhs.initialized;
83 this->computeMagnitude = rhs.computeMagnitude;
84 this->computePhase = rhs.computePhase;
85 this->enableZeroPadding = rhs.enableZeroPadding;
86 this->windowSize = rhs.windowSize;
87 this->windowFunction = rhs.windowFunction;
88 this->averagePower = 0;
91 if( rhs.initialized ){
92 this->init(rhs.windowSize,rhs.windowFunction,rhs.computeMagnitude,rhs.computePhase);
95 for(UINT i=0; i<this->windowSize; i++){
96 this->fftReal[i] = rhs.fftReal[i];
97 this->fftImag[i] = rhs.fftImag[i];
98 this->magnitude[i] = rhs.magnitude[i];
99 this->phase[i] = rhs.phase[i];
100 this->power[i] = rhs.power[i];
107 bool FastFourierTransform::init(
const unsigned int windowSize,
const unsigned int windowFunction,
const bool computeMagnitude,
const bool computePhase,
const bool enableZeroPadding){
113 if( !isPowerOfTwo( windowSize ) ){
117 if( windowFunction != RECTANGULAR_WINDOW && windowFunction != BARTLETT_WINDOW &&
118 windowFunction != HAMMING_WINDOW && windowFunction != HANNING_WINDOW ){
124 this->windowSize = windowSize;
125 this->windowFunction = windowFunction;
126 this->computeMagnitude = computeMagnitude;
127 this->computePhase = computePhase;
128 this->enableZeroPadding = enableZeroPadding;
131 fftReal.
resize( windowSize );
132 fftImag.
resize( windowSize );
133 tmpReal.
resize( windowSize/2 );
134 tmpImag.
resize( windowSize/2 );
135 magnitude.
resize( windowSize );
136 phase.
resize( windowSize );
137 power.
resize( windowSize );
141 for(UINT i=0; i<windowSize/2; i++){
145 for(UINT i=0; i<windowSize; i++){
159 bool FastFourierTransform::computeFFT(
VectorFloat &data ){
166 if( !enableZeroPadding ){
167 if( (
unsigned int)data.size() != windowSize ){
168 errorLog <<
"The size of the data Vector (" << data.size() <<
") does not match the windowSize: " << windowSize << std::endl;
172 if( (
unsigned int)data.size() > windowSize ){
173 errorLog <<
"The size of the data Vector (" << data.size() <<
") is greater than the windowSize: " << windowSize << std::endl;
179 if( !windowData( data ) ){
184 if( enableZeroPadding ){
185 if( ((
unsigned int)data.size()) != windowSize ){
186 const unsigned int oldSize = (
unsigned int)data.size();
187 data.
resize( windowSize );
188 for(
unsigned int i=oldSize; i<windowSize; i++){
195 realFFT(data, &fftReal[0], &fftImag[0]);
199 for(
unsigned int i = 0; i<windowSize/2; i++){
201 if( computeMagnitude ){
202 power[i] = fftReal[i]*fftReal[i] + fftImag[i]*fftImag[i];
203 averagePower += power[i];
204 magnitude[i] = 2.0*sqrt( power[i] );
208 phase[i] = atan2(fftImag[i],fftReal[i]);
213 averagePower = averagePower / (Float)(windowSize/2);
218 bool FastFourierTransform::windowData(
VectorFloat &data ){
220 const unsigned int N = (
unsigned int)data.size();
221 const unsigned int K = N/2;
223 switch( windowFunction ){
224 case RECTANGULAR_WINDOW:
227 case BARTLETT_WINDOW:
228 for(
unsigned int i=0; i<K; i++) {
229 data[i] *= (i / (Float) (K));
230 data[i + K] *= (1.0 - (i / (Float)K));
235 for(
unsigned int i=0; i<N; i++)
236 data[i] *= 0.54 - 0.46 * cos(2 * PI * i / (N - 1));
240 for(
unsigned int i=0; i <N; i++)
241 data[i] *= 0.50 - 0.50 * cos(2 * PI * i / (N - 1));
252 VectorFloat FastFourierTransform::getMagnitudeData(){
256 const unsigned int N = windowSize/2;
259 for(
unsigned int i=0; i<N; i++){
260 magnitudeData[i] = magnitude[i];
263 return magnitudeData;
269 const unsigned int N = windowSize/2;
272 for(
unsigned int i=0; i<N; i++){
273 phaseData[i] = phase[i];
284 for(
unsigned int i=0; i<windowSize/2; i++){
285 powerData[i] = power[i];
291 Float FastFourierTransform::getAveragePower(){
295 Float* FastFourierTransform::getMagnitudeDataPtr(){
296 return &magnitude[0];
299 Float* FastFourierTransform::getPhaseDataPtr(){
303 Float* FastFourierTransform::getPowerDataPtr(){
322 bool FastFourierTransform::realFFT(
const VectorFloat &realIn, Float *realOut, Float *imagOut ){
323 int NumSamples = (int)windowSize;
324 int Half = NumSamples / 2;
327 Float theta = PI / Half;
329 for (i = 0; i < Half; i++) {
330 tmpReal[i] = realIn[2 * i];
331 tmpImag[i] = realIn[2 * i + 1];
334 if( !
FFT(Half, 0, &tmpReal[0], &tmpImag[0], realOut, imagOut) ){
338 Float wtemp = Float(sin(0.5 * theta));
340 Float wpr = -2.0 * wtemp * wtemp;
341 Float wpi = Float (sin(theta));
342 Float wr = 1.0 + wpr;
347 Float h1r, h1i, h2r, h2i;
349 for (i = 1; i < Half / 2; i++) {
353 h1r = 0.5 * (realOut[i] + realOut[i3]);
354 h1i = 0.5 * (imagOut[i] - imagOut[i3]);
355 h2r = 0.5 * (imagOut[i] + imagOut[i3]);
356 h2i = -0.5 * (realOut[i] - realOut[i3]);
358 realOut[i] = h1r + wr * h2r - wi * h2i;
359 imagOut[i] = h1i + wr * h2i + wi * h2r;
360 realOut[i3] = h1r - wr * h2r + wi * h2i;
361 imagOut[i3] = -h1i + wr * h2i + wi * h2r;
363 wr = (wtemp = wr) * wpr - wi * wpi + wr;
364 wi = wi * wpr + wtemp * wpi + wi;
367 realOut[0] = (h1r = realOut[0]) + imagOut[0];
368 imagOut[0] = h1r - imagOut[0];
373 bool FastFourierTransform::FFT(
int numSamples,
bool inverseTransform,Float *realIn, Float *imagIn, Float *realOut, Float *imagOut){
376 int BlockSize, BlockEnd;
378 Float angle_numerator = 2.0 * PI;
381 if( !isPowerOfTwo(numSamples) ) {
382 fprintf(stderr,
"%d is not a power of two\n", numSamples);
386 if( bitTable.size() == 0 ) initFFT();
388 if( inverseTransform ) angle_numerator = -angle_numerator;
390 NumBits = numberOfBitsNeeded(numSamples);
393 for(i = 0; i < numSamples; i++) {
394 j = fastReverseBits(i, NumBits);
395 realOut[j] = realIn[i];
396 imagOut[j] = (imagIn == NULL) ? 0.0 : imagIn[i];
401 for (BlockSize = 2; BlockSize <= numSamples; BlockSize <<= 1) {
403 Float delta_angle = angle_numerator / (Float) BlockSize;
405 Float sm2 = sin(-2 * delta_angle);
406 Float sm1 = sin(-delta_angle);
407 Float cm2 = cos(-2 * delta_angle);
408 Float cm1 = cos(-delta_angle);
410 Float ar0, ar1, ar2, ai0, ai1, ai2;
412 for (i = 0; i < numSamples; i += BlockSize) {
419 for (j = i, n = 0; n < BlockEnd; j++, n++) {
429 tr = ar0 * realOut[k] - ai0 * imagOut[k];
430 ti = ar0 * imagOut[k] + ai0 * realOut[k];
432 realOut[k] = realOut[j] - tr;
433 imagOut[k] = imagOut[j] - ti;
440 BlockEnd = BlockSize;
444 if( inverseTransform ){
445 Float denom = (Float) numSamples;
447 for(i = 0; i < numSamples; i++) {
456 int FastFourierTransform::numberOfBitsNeeded(
int powerOfTwo)
458 for (
int i = 0;; i++){
459 if (powerOfTwo & (1 << i)){
466 int FastFourierTransform::reverseBits(
int index,
int numBits)
470 for (i = rev = 0; i < numBits; i++) {
471 rev = (rev << 1) | (index & 1);
478 void FastFourierTransform::initFFT()
480 bitTable.
resize( MAX_FAST_BITS );
483 for (
int b = 1; b <= MAX_FAST_BITS; b++) {
485 bitTable[b - 1].
resize(len);
487 for (
int i = 0; i < len; i++)
488 bitTable[b - 1][i] = reverseBits(i, b);
494 inline int FastFourierTransform::fastReverseBits(
const int i,
const int numBits)
496 if (numBits <= MAX_FAST_BITS)
497 return bitTable[numBits - 1][i];
499 return reverseBits(i, numBits);
502 inline bool FastFourierTransform::isPowerOfTwo(
const unsigned int x){
503 if (x < 2)
return false;
504 if (x & (x - 1))
return false;
virtual bool resize(const unsigned int size)