26 #include "FastFourierTransform.h"
31 FastFourierTransform::FastFourierTransform(){
33 computeMagnitude =
true;
36 windowFunction = RECTANGULAR_WINDOW;
39 infoLog.setProceedingText(
"[FastFourierTransform]");
40 warningLog.setProceedingText(
"[WARNING FastFourierTransform]");
41 errorLog.setProceedingText(
"[ERROR FastFourierTransform]");
46 FastFourierTransform::FastFourierTransform(
const FastFourierTransform &rhs){
48 this->initialized = rhs.initialized;
49 this->computeMagnitude = rhs.computeMagnitude;
50 this->computePhase = rhs.computePhase;
51 this->windowSize = rhs.windowSize;
52 this->windowFunction = rhs.windowFunction;
53 this->averagePower = 0;
55 this->infoLog = rhs.infoLog;
56 this->warningLog = rhs.warningLog;
57 this->errorLog = rhs.errorLog;
59 if( rhs.initialized ){
60 this->init(rhs.windowSize,rhs.windowFunction,rhs.computeMagnitude,rhs.computePhase);
63 for(UINT i=0; i<this->windowSize; i++){
64 this->fftReal[i] = rhs.fftReal[i];
65 this->fftImag[i] = rhs.fftImag[i];
66 this->magnitude[i] = rhs.magnitude[i];
67 this->phase[i] = rhs.phase[i];
68 this->power[i] = rhs.power[i];
73 FastFourierTransform::~FastFourierTransform() {
77 FastFourierTransform& FastFourierTransform::operator=(
const FastFourierTransform &rhs){
80 this->initialized = rhs.initialized;
81 this->computeMagnitude = rhs.computeMagnitude;
82 this->computePhase = rhs.computePhase;
83 this->windowSize = rhs.windowSize;
84 this->windowFunction = rhs.windowFunction;
85 this->averagePower = 0;
88 if( rhs.initialized ){
89 this->init(rhs.windowSize,rhs.windowFunction,rhs.computeMagnitude,rhs.computePhase);
92 for(UINT i=0; i<this->windowSize; i++){
93 this->fftReal[i] = rhs.fftReal[i];
94 this->fftImag[i] = rhs.fftImag[i];
95 this->magnitude[i] = rhs.magnitude[i];
96 this->phase[i] = rhs.phase[i];
97 this->power[i] = rhs.power[i];
104 bool FastFourierTransform::init(
const unsigned int windowSize,
const unsigned int windowFunction,
const bool computeMagnitude,
const bool computePhase){
110 if( !isPowerOfTwo( windowSize ) ){
114 if( windowFunction != RECTANGULAR_WINDOW && windowFunction != BARTLETT_WINDOW &&
115 windowFunction != HAMMING_WINDOW && windowFunction != HANNING_WINDOW ){
121 this->windowSize = windowSize;
122 this->windowFunction = windowFunction;
123 this->computeMagnitude = computeMagnitude;
124 this->computePhase = computePhase;
127 fftReal.resize( windowSize );
128 fftImag.resize( windowSize );
129 tmpReal.resize( windowSize/2 );
130 tmpImag.resize( windowSize/2 );
131 magnitude.resize( windowSize );
132 phase.resize( windowSize );
133 power.resize( windowSize );
137 for(UINT i=0; i<windowSize/2; i++){
141 for(UINT i=0; i<windowSize; i++){
155 bool FastFourierTransform::computeFFT( VectorDouble &data ){
162 if( !windowData( data ) ){
167 realFFT(data, &fftReal[0], &fftImag[0]);
171 for(
unsigned int i = 0; i<windowSize/2; i++){
173 if( computeMagnitude ){
174 power[i] = fftReal[i]*fftReal[i] + fftImag[i]*fftImag[i];
175 averagePower += power[i];
176 magnitude[i] = 2.0*sqrt( power[i] );
180 phase[i] = atan2(fftImag[i],fftReal[i]);
185 averagePower = averagePower / (double)(windowSize/2);
190 bool FastFourierTransform::windowData( VectorDouble &data ){
192 if( data.size() != windowSize ){
193 errorLog <<
"The size of the data vector (" << data.size() <<
") does not match the windowSize: " << windowSize << endl;
197 switch( windowFunction ){
198 case RECTANGULAR_WINDOW:
201 case BARTLETT_WINDOW:
202 for(
unsigned int i=0; i<windowSize/2; i++) {
203 data[i] *= (i / (double) (windowSize / 2));
204 data[i + (windowSize/2)] *= (1.0 - (i / (
double) (windowSize/2)));
209 for(
unsigned int i=0; i<windowSize; i++)
210 data[i] *= 0.54 - 0.46 * cos(2 * PI * i / (windowSize - 1));
214 for(
unsigned int i=0; i <windowSize; i++)
215 data[i] *= 0.50 - 0.50 * cos(2 * PI * i / (windowSize - 1));
226 VectorDouble FastFourierTransform::getMagnitudeData(){
228 if( !initialized )
return VectorDouble();
230 VectorDouble magnitudeData(windowSize/2);
232 for(
unsigned int i=0; i<windowSize/2; i++){
233 magnitudeData[i] = magnitude[i];
236 return magnitudeData;
239 VectorDouble FastFourierTransform::getPhaseData(){
240 if( !initialized )
return VectorDouble();
242 VectorDouble phaseData(windowSize/2);
244 for(
unsigned int i=0; i<windowSize/2; i++){
245 phaseData[i] = phase[i];
251 VectorDouble FastFourierTransform::getPowerData(){
252 if( !initialized )
return VectorDouble();
254 VectorDouble powerData(windowSize/2);
256 for(
unsigned int i=0; i<windowSize/2; i++){
257 powerData[i] = power[i];
263 double FastFourierTransform::getAveragePower(){
267 double* FastFourierTransform::getMagnitudeDataPtr(){
268 return &magnitude[0];
271 double* FastFourierTransform::getPhaseDataPtr(){
275 double* FastFourierTransform::getPowerDataPtr(){
294 bool FastFourierTransform::realFFT(
const VectorDouble &realIn,
double *realOut,
double *imagOut ){
295 int NumSamples = (int)windowSize;
296 int Half = NumSamples / 2;
299 double theta = PI / Half;
301 for (i = 0; i < Half; i++) {
302 tmpReal[i] = realIn[2 * i];
303 tmpImag[i] = realIn[2 * i + 1];
306 if( !FFT(Half, 0, &tmpReal[0], &tmpImag[0], realOut, imagOut) ){
310 double wtemp = double(sin(0.5 * theta));
312 double wpr = -2.0 * wtemp * wtemp;
313 double wpi = double (sin(theta));
314 double wr = 1.0 + wpr;
319 double h1r, h1i, h2r, h2i;
321 for (i = 1; i < Half / 2; i++) {
325 h1r = 0.5 * (realOut[i] + realOut[i3]);
326 h1i = 0.5 * (imagOut[i] - imagOut[i3]);
327 h2r = 0.5 * (imagOut[i] + imagOut[i3]);
328 h2i = -0.5 * (realOut[i] - realOut[i3]);
330 realOut[i] = h1r + wr * h2r - wi * h2i;
331 imagOut[i] = h1i + wr * h2i + wi * h2r;
332 realOut[i3] = h1r - wr * h2r + wi * h2i;
333 imagOut[i3] = -h1i + wr * h2i + wi * h2r;
335 wr = (wtemp = wr) * wpr - wi * wpi + wr;
336 wi = wi * wpr + wtemp * wpi + wi;
339 realOut[0] = (h1r = realOut[0]) + imagOut[0];
340 imagOut[0] = h1r - imagOut[0];
345 bool FastFourierTransform::FFT(
int numSamples,
bool inverseTransform,
double *realIn,
double *imagIn,
double *realOut,
double *imagOut){
348 int BlockSize, BlockEnd;
350 double angle_numerator = 2.0 * PI;
353 if( !isPowerOfTwo(numSamples) ) {
354 fprintf(stderr,
"%d is not a power of two\n", numSamples);
358 if( bitTable.size() == 0 ) initFFT();
360 if( inverseTransform ) angle_numerator = -angle_numerator;
362 NumBits = numberOfBitsNeeded(numSamples);
365 for(i = 0; i < numSamples; i++) {
366 j = fastReverseBits(i, NumBits);
367 realOut[j] = realIn[i];
368 imagOut[j] = (imagIn == NULL) ? 0.0 : imagIn[i];
373 for (BlockSize = 2; BlockSize <= numSamples; BlockSize <<= 1) {
375 double delta_angle = angle_numerator / (double) BlockSize;
377 double sm2 = sin(-2 * delta_angle);
378 double sm1 = sin(-delta_angle);
379 double cm2 = cos(-2 * delta_angle);
380 double cm1 = cos(-delta_angle);
382 double ar0, ar1, ar2, ai0, ai1, ai2;
384 for (i = 0; i < numSamples; i += BlockSize) {
391 for (j = i, n = 0; n < BlockEnd; j++, n++) {
401 tr = ar0 * realOut[k] - ai0 * imagOut[k];
402 ti = ar0 * imagOut[k] + ai0 * realOut[k];
404 realOut[k] = realOut[j] - tr;
405 imagOut[k] = imagOut[j] - ti;
412 BlockEnd = BlockSize;
416 if( inverseTransform ){
417 double denom = (double) numSamples;
419 for(i = 0; i < numSamples; i++) {
428 int FastFourierTransform::numberOfBitsNeeded(
int powerOfTwo)
430 for (
int i = 0;; i++){
431 if (powerOfTwo & (1 << i)){
438 int FastFourierTransform::reverseBits(
int index,
int numBits)
442 for (i = rev = 0; i < numBits; i++) {
443 rev = (rev << 1) | (index & 1);
450 void FastFourierTransform::initFFT()
452 bitTable.resize( MAX_FAST_BITS );
455 for (
int b = 1; b <= MAX_FAST_BITS; b++) {
457 bitTable[b - 1].resize(len);
459 for (
int i = 0; i < len; i++)
460 bitTable[b - 1][i] = reverseBits(i, b);
466 inline int FastFourierTransform::fastReverseBits(
const int i,
const int numBits)
468 if (numBits <= MAX_FAST_BITS)
469 return bitTable[numBits - 1][i];
471 return reverseBits(i, numBits);
474 inline bool FastFourierTransform::isPowerOfTwo(
const unsigned int x){
475 if (x < 2)
return false;
476 if (x & (x - 1))
return false;