1
votes

FFT works fine, but when I want to take IFFT I always see the same graph from its results. Results are complex and graph always the same regardless of the original signal.

in real part graph is a -sin with period = frame size

in imaginary part it is a -cos with the same period

Where can be a problem?

original signal: original signal

IFFT real value (on pics are only half of frame): IFFT real value

Algorithm FFT that I use.

double** FFT(double** f, int s, bool inverse) {
    if (s == 1) return f;
    int sH = s / 2;

    double** fOdd = new double*[sH];
    double** fEven = new double*[sH];
    for (int i = 0; i < sH; i++) {
       int j = 2 * i;
       fOdd[i] = f[j];
       fEven[i] = f[j + 1];
    }

    double** sOdd = FFT(fOdd, sH, inverse);
    double** sEven = FFT(fEven, sH, inverse);

    double**spectr = new double*[s];

    double arg = inverse ? DoublePI / s : -DoublePI / s;
    double*oBase = new double[2]{ cos(arg),sin(arg) };
    double*o = new double[2]{ 1,0 };

    for (int i = 0; i < sH; i++) {
        double* sO1 = Mul(o, sOdd[i]);

        spectr[i] = Sum(sEven[i], sO1);
        spectr[i + sH] = Dif(sEven[i], sO1);

        o = Mul(o, oBase);
    }

    return spectr;
}
2
Side-note: your code allocates a lot of objects on the heap using new but your code never calls delete without changing ownership of the objects before they leave scope - so your program will leak memory.Dai
@Dai, thank you for answer, I'll try to fix it when it will work. But now it's not the biggest my problem.Pavlo Holotiuk
imho it is the biggest problem. All that pointers and news make the code unnecesarily hard to read and to find bugs like the one you have463035818_is_not_a_number
Anyone would think that STL had never been invented...Paul R

2 Answers

2
votes

The "butterfly" portion is applying the coefficients incorrectly:

for (int i = 0; i < sH; i++) {
    double* sO1 = sOdd[i];
    double* sE1 = Mul(o, sEven[i]);

    spectr[i] = Sum(sO1, sE1);
    spectr[i + sH] = Dif(sO1, sE1);

    o = Mul(o, oBase);
}

Side Note:

I kept your notation but it makes things confusing:

fOdd has indexes 0, 2, 4, 6, ... so it should be fEven

fEven has indexes 1, 3, 5, 7, ... so it should be fOdd

really sOdd should be sLower and sEven should be sUpper since they correspond to the 0:s/2 and s/2:s-1 elements of the spectrum respectively:

sLower = FFT(fEven, sH, inverse); // fEven is 0, 2, 4, ...
sUpper = FFT(fOdd, sH, inverse); // fOdd is 1, 3, 5, ...

Then the butterfly becomes:

for (int i = 0; i < sH; i++) {
    double* sL1 = sLower[i];
    double* sU1 = Mul(o, sUpper[i]);

    spectr[i] = Sum(sL1, sU1);
    spectr[i + sH] = Dif(sL1, sU1);

    o = Mul(o, oBase);
}

When written like this it is easier to compare to this pseudocode example on wikipedia.

And @Dai is correct you are going to leak a lot of memory

1
votes

Regarding the memory, you can use the std::vector to encapsulate dynamically-allocated arrays and to ensure they're deallocated when execution leaves scope. You could use unique_ptr<double[]> but the performance gains are not worth it IMO and you lose the safety of the at() method.

(Based on @Robb's answer)

A few other tips:

  • Avoid cryptic identifiers - programs should be readable, and names like "f" and "s" make your program harder to read and maintain.
  • Type-based Hungarian notation is frowned upon as modern editors show type information automatically so it adds unnecessary complication to identifier names.
  • Use size_t for indexes, not int
  • The STL is your friend, use it!
  • Preemptively prevent bugs by using const to prevent accidental mutation of read-only data.

Like so:

#include <vector>

using namespace std;

vector<double> fastFourierTransform(const vector<double> signal, const bool inverse) {

    if( signal.size() < 2 ) return signal;
    const size_t half = signal.size() / 2;

    vector<double> lower; lower.reserve( half );
    vector<double> upper; upper.reserve( half );

    bool isEven = true;
    for( size_t i = 0; i < signal.size(); i++ ) {
        if( isEven ) lower.push_back( signal.at( i ) );
        else         upper.push_back( signal.at( i ) );

        isEven = !isEven;
    }

    vector<double> lowerFft = fastFourierTransform( lower, inverse );
    vector<double> upperFft = fastFourierTransform( upper, inverse );

    vector<double> result;
    result.reserve( signal.size() );

    double arg = ( inverse ? 1 : -1 ) * ( DoublePI / signal.size() );

    // Ideally these should be local `double` values passed directly into `Mul`.
    unique_ptr<double[]> oBase = make_unique<double[]>( 2 );
    oBase[0] = cos(arg);
    oBase[1] = sin(arg);

    unique_ptr<double[]> o = make_unique<double[]>( 2 );
    o[0] = 0;
    o[1] = 0;

    for( size_t i = 0; i < half; i++ ) {

        double* lower1 = lower.at( i );
        double* upper1 = Mul( o, upper.at( i ) );

        result.at( i        ) = Sum( lower1, upper1 );
        result.at( i + half ) = Dif( lower1, upper1 );

        o = Mul( o, oBase );
    }

    // My knowledge of move-semantics of STL containers is a bit rusty - so there's probably a better way to return the output 'result' vector.
    return result;
}