C++

Forum

 
Help with code analysis

alexBB (131)
Hi there,

I do numerical analysis in 2D. Part of it is Fast Fourier Transform (FFT). I've chosen so called Librow FFT.

http://www.librow.com/articles/article-10

It is really fast. However, recently I discovered that it reverses the image. It works fine with symmetrical figures when I test them but if I use asymmetrical figure like a Scalene triangle leaning to one side, when I restore the image with Inverse FFT, it will look the other way.

I want to explain why I need the Forward Transform to be correct. I actually don't need the Inverse transform, I use it for control only, I need coefficients. I wrote a Discrete Fourier Transform DFT without accelerating algorithms; it works fine, albeit very slow, and shows no inversion, HOWEVER, the Fourier coefficients after the Forward Transform are different from Librow's. I really need Fourier coefficients obtained fast but such that they will match my DFT.

I am trying to analyze the code and find out how it happens but so far had no success. I wonder if anybody will try to help.

Thanks, -
Last edited on
lastchance (6980)
Are you looking for 2D or 1D Fourier transforms?

Note that, if {Fn} is the Discrete Fourier transform of {fk} (irrespective of whether it is "fast" or not) then the inverse can be computed (in 1D) as
f = (1/N) conj( DFT( conj F ) )
where conj is complex conjugate.

In other words, you never have to code an inverse transform - you can simply use the routines for forward transform with some appropriate conjugation. So, try the library algorithm without using the inverse routine at all (it might be buggy) and see if your image is reversed.

To check coefficients just try a very small array by hand and see whether your (and the library coefficients) agree with them.

Fourier transforms sometimes have non-standard scaling factors and sign conventions.


alexBB (131)
Thank you but I don't think it is an answer to my question. Yes, I know that but I simply change the sign for sin (theta) from minus to plus. In a sense I take the complex conjugate. The problem is the whole set of Librow's routines reverses the image and this means, I think, that coefficients are wrong. Yes, I use 2-D but the second dimension is Associated Legendre Polynomials. It is all on 2-sphere.

There is a chance I am close to the answer. I emailed to one of the programmers who wrote the software and he responded. He tested as you said a very small array and there was no inversion, HOWEVER, he used different routines. He used inplace routines, I used routines with two separate arrays, for Input and Output. That could be the difference.

Thank you for posting. I will keep you posted :-).

Thanks.
lastchance (6980)
If you want simple 1-D DFT and/or FFT routines alone then the following header and cpp file may help. At least you can check with your own routine.

Call FFT to do the forward transform and iFFT to do the reverse. It will use the recursion to speed it up for every factor of 2 in N. If N is odd it will just do a DFT.

Original 1-d array is f[]; discrete Fourier transform is ftilde[]. Both arrays are complex<double>.

fft.h -
1
2
3
4
5
6
7
8
9
 
#include <complex>
using std::complex;

const double pi = 3.14159265358979323846;

void FFT ( int N, complex<double> *f     , complex<double> *ftilde );
void iFFT( int N, complex<double> *ftilde, complex<double> *f      );
void DFT ( int N, complex<double> *f     , complex<double> *ftilde );
void iDFT( int N, complex<double> *ftilde, complex<double> *f      );



fft.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
 
#include <complex>
#include "fft.h"
using namespace std;


//======================================================================


void FFT( int N, complex<double> *f, complex<double> *ftilde )       // fast fourier transform
{
   if ( N % 2 == 0 )                                                 // recursion if a multiple of 2
   { 
      int N2 = N / 2;
      complex<double> *even = new complex<double>[N2];
      complex<double> *odd  = new complex<double>[N2];
      complex<double> *g    = new complex<double>[N2];
      complex<double> *h    = new complex<double>[N2];

      for ( int m = 0; m < N2; m++ )
      {
         even[m] = f[2*m  ];
         odd [m] = f[2*m+1];
      }

      FFT( N2, even, g );
      FFT( N2, odd , h );

      complex<double> w = polar( 1.0, -2.0 * pi / N );
      complex<double> wn = complex<double>( 1.0 );
      for ( int n = 0; n < N2; n++ )
      {
         int nplus = n + N2;
         ftilde[n]     = g[n] + wn * h[n];
         ftilde[nplus] = g[n] - wn * h[n];
         wn *= w;
      }

      delete [] even;
      delete [] odd ;
      delete [] g   ;
      delete [] h   ;
   }

   else                                                              // otherwise just do a DFT
   {
      DFT( N, f, ftilde );
   }
}


//======================================================================


void iFFT( int N, complex<double> *ftilde, complex<double> *f )      // inverse fast fourier transform
{
   complex<double> * ftildeConjugate = new complex<double>[N];

   for ( int m = 0; m < N; m++ ) ftildeConjugate[m] = conj( ftilde[m] );

   FFT( N, ftildeConjugate, f );

   for ( int m = 0; m < N; m++ ) f[m] = conj( f[m] ) / (double) N;

   delete [] ftildeConjugate;
}


//======================================================================


void DFT( int N, complex<double> *f, complex<double> *ftilde )       // discrete Fourier transform
{
   for ( int n = 0; n < N; n++ )
   {
      ftilde[n] = complex<double>( 0.0 );
      complex<double> w = polar( 1.0, -2.0 * pi * n / N );
      complex<double> wm( 1.0 );
      for ( int m = 0; m < N; m++ )
      {
         ftilde[n] += f[m] * wm;
         wm *= w;
      }
   }
}


//======================================================================


void iDFT( int N, complex<double> *ftilde, complex<double> *f )      // inverse discrete Fourier transform
{
   complex<double> * ftildeConjugate = new complex<double>[N];

   for ( int m = 0; m < N; m++ ) ftildeConjugate[m] = conj( ftilde[m] );
   DFT( N, ftildeConjugate, f );
   for ( int m = 0; m < N; m++ ) f[m] = conj( f[m] ) / (double) N;

   delete [] ftildeConjugate;
}
lastchance (6980)
Here's a sample code to use with the routines above. It discrete-Fourier-transforms a function (in this case, an asymmetric, right-leaning triangle) and then inverts the transform to check that the original function is recovered (columns 2 and 6 in bold in the output - the real parts of the original function).

All the array elements have real and imaginary parts; the initial function starts pure real.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
 
#include <iostream>
#include <iomanip>
#include <complex>
#include "fft.h"
using namespace std;


double func( double x );     // returns function to be Fourier-transformed


int main()
{
   const int N = 16;
   complex<double> f[N];     // original function
   complex<double> ft[N];    // discrete Fourier transform
   complex<double> g[N];     // inverted Fourier transform - hopefully the same as f[]

   // Populate f with some values
   double L = 1.0;
   double dx = L / N;
   for ( int i = 0; i < N; i++ ) f[i] = func( i * dx );

   // Forward transform
   FFT( N, f, ft );

   // Invert transform
   iFFT( N, ft, g );

   // Tabulate output to check
   #define SH << setw( 15 ) <<
   #define SD << setw( 15 ) << fixed << setprecision( 6 ) <<
   #define NL << endl;
   cout SH " " SH "Original"  SH " " SH "Transform" SH " " SH "Inverted" NL
   cout SH "x";
   for ( int i = 0; i < 3; i++ ) cout SH "Real" SH "Imag";
   cout NL
   for ( int i = 0; i < N; i++ ) 
   {
      cout SD i * dx SD f[i].real() SD f[i].imag() SD ft[i].real() SD ft[i].imag() SD g[i].real() SD g[i].imag() NL
   }
}


double func( double x )      // right-leaning triangle
{
   if ( x > 0.1 && x <= 0.7 ) return       ( x - 0.1 ) / 0.6;
   if ( x > 0.7 && x <  0.9 ) return 1.0 - ( x - 0.7 ) / 0.2;
   return 0.0;
}


                      Original                     Transform                      Inverted
              x           Real           Imag           Real           Imag           Real           Imag
       0.000000       0.000000       0.000000       6.416667       0.000000       0.000000      -0.000000
       0.062500       0.000000       0.000000      -3.010477       1.780169       0.000000       0.000000
       0.125000       0.041667       0.000000      -0.745812      -0.716350       0.041667      -0.000000
       0.187500       0.145833       0.000000       0.340956      -0.304765       0.145833       0.000000
       0.250000       0.250000       0.000000       0.166667       0.125000       0.250000      -0.000000
       0.312500       0.354167       0.000000      -0.024374      -0.011872       0.354167      -0.000000
       0.375000       0.458333       0.000000       0.079146      -0.049683       0.458333       0.000000
       0.437500       0.562500       0.000000       0.027227       0.073062       0.562500       0.000000
       0.500000       0.666667       0.000000      -0.083333       0.000000       0.666667       0.000000
       0.562500       0.770833       0.000000       0.027227      -0.073062       0.770833      -0.000000
       0.625000       0.875000       0.000000       0.079146       0.049683       0.875000       0.000000
       0.687500       0.979167       0.000000      -0.024374       0.011872       0.979167      -0.000000
       0.750000       0.750000       0.000000       0.166667      -0.125000       0.750000       0.000000
       0.812500       0.437500       0.000000       0.340956       0.304765       0.437500       0.000000
       0.875000       0.125000       0.000000      -0.745812       0.716350       0.125000      -0.000000
       0.937500       0.000000       0.000000      -3.010477      -1.780169       0.000000       0.000000



EDIT: OK, I modified this code to run with the Librow FFT libraries instead (a bit of a nuisance: they have their own non-standard complex class).

It gave exactly the same answers, and inverted correctly - so I have confidence in their implementation (and mine).
Last edited on
alexBB (131)
@lastchance, thank you very much. I've seen the code you posted but now I will try it. The Librow engineer I was talking about helped me to find a bug in my code. I modified his software incorrectly and now everything works OK. Thank you very much. It is an exciting area.
lastchance (6980)
Hi @alexBB,
some slightly improved routines for you - judicious use of pointers enables less use of new/delete. I don't go in for rearranging indices as in the Librow code.

fft.h 
1
2
3
4
5
6
7
8
9
 
#include <complex>
using std::complex;

const double pi = 3.14159265358979323846;
typedef complex<double> cmplx;

void FFT ( cmplx* f     , cmplx* ftilde, int N, int step = 1 );
void DFT ( cmplx* f     , cmplx* ftilde, int N, int step = 1 );
void iFFT( cmplx* ftilde, cmplx* f     , int N );



fft.cpp 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
 
#include <complex>
#include "fft.h"
using namespace std;

typedef complex<double> cmplx;

//======================================================================


void FFT( cmplx* f, cmplx* ftilde, int N, int step )                 // fast fourier transform
{
   if ( N % 2 == 0 )                                                 // recursion if a multiple of 2
   { 
      int N2 = N / 2;
      cmplx* g = new cmplx[N2];
      cmplx* h = new cmplx[N2];

      int step2 = step + step;
      FFT( f       , g, N2, step2 );
      FFT( f + step, h, N2, step2 );

      cmplx w = polar( 1.0, -2.0 * pi / N );
      cmplx wn = 1.0;
      for ( int n = 0; n < N2; n++ )
      {
         ftilde[n]    = g[n] + wn * h[n];
         ftilde[n+N2] = g[n] - wn * h[n];
         wn *= w;
      }

      delete [] g;
      delete [] h;
   }

   else                                                              // otherwise just do a DFT
   {
      DFT( f, ftilde, N, step );
   }
}

//======================================================================

void DFT( cmplx* f, cmplx* ftilde, int N, int step )                 // basic discrete Fourier transform
{
   for ( int n = 0; n < N; n++ )
   {
      ftilde[n] = 0.0;
      cmplx w = polar( 1.0, -2.0 * pi * n / N );
      cmplx wm = 1.0;
      for ( int m = 0, mpos = 0; m < N; m++, mpos += step )
      {
         ftilde[n] += f[mpos] * wm;
         wm *= w;
      }
   }
}

//======================================================================

void iFFT( cmplx* ftilde, cmplx* f, int N )                          // inverse fast fourier transform
{
   cmplx* ftildeConjugate = new cmplx[N];

   for ( int m = 0; m < N; m++ ) ftildeConjugate[m] = conj( ftilde[m] );

   FFT( ftildeConjugate, f, N );

   double factor = 1.0 / N;
   for ( int m = 0; m < N; m++ ) f[m] = factor * conj( f[m] );

   delete [] ftildeConjugate;
}


testfft.cpp 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
 
#include <iostream>
#include <iomanip>
#include <complex>
#include <ctime>
#include "fft.h"
using namespace std;

typedef complex<double> cmplx;

double func( double x );     // returns function to be Fourier-transformed


int main()
{
   int N = 16;
   cout << "Input N: ";
   cin >> N;

   cmplx* f  = new cmplx[N];     // original function
   cmplx* ft = new cmplx[N];     // discrete Fourier transform
   cmplx* g  = new cmplx[N];     // inverted Fourier transform - hopefully the same as f[]

   // Populate f with some values
   double L = 1.0;
   double dx = L / N;
   for ( int i = 0; i < N; i++ ) f[i] = func( i * dx );

   // Start the clock
   clock_t t = clock();

   // Forward transform
   FFT( f, ft, N );

   // Invert transform
   iFFT( ft, g, N );

   // End the clock
   cout << "Processor time taken = " << (double)( clock() - t ) * 1000 / CLOCKS_PER_SEC << " ms\n";

   // Tabulate output to check
   #define SH << setw( 15 ) <<
   #define SD << setw( 15 ) << fixed << setprecision( 6 ) <<
   #define NL << endl;
   cout SH " " SH "Original"  SH " " SH "Transform" SH " " SH "Inverted" NL
   cout SH "x";
   for ( int i = 0; i < 3; i++ ) cout SH "Real" SH "Imag";
   cout NL
   int step = N / 16;    if ( step < 1 ) step = 1;
   for ( int i = 0; i < N; i += step )
   {
      cout SD i * dx SD f[i].real() SD f[i].imag() SD ft[i].real() SD ft[i].imag() SD g[i].real() SD g[i].imag() NL
   }
}


double func( double x )      // right-leaning triangle
{
   if ( x > 0.1 && x <= 0.7 ) return       ( x - 0.1 ) / 0.6;
   if ( x > 0.7 && x <  0.9 ) return 1.0 - ( x - 0.7 ) / 0.2;
   return 0.0;
}



Lest you doubt the advantages of the fast Fourier transform, note:
Input N: 16384
Processor time taken = 74 ms
Input N: 16383
Processor time taken = 27917 ms

The latter has no factor of 2, so no possibility of acceleration.
Last edited on
alexBB (131)
@lastchance, thank you for the code. Sorry for delayed response. Circumstances :-) BTW, a Librow engineer with whom I managed to communicate praised your code. I will try your code now. Marking the thread solved.

Thank you again, - Alex
alexBB (131)
@lastchance, I need your help again. Your code is fast and restores the function well. However, paradoxically, it is not what I need. The following simple code which I wrote myself months ago is a DFT. I know you have a similar routine as well.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
 
 class homeFourier 
{   // straightforward DFT, not FFT, works well June 1, 2017
  public : void ForwardFourier (cmplx *const Input, cmplx *const &Output, int NN)
  {
    cmplx XX,prom;   // WORKS VERY WELL June 1, 2017
    for (int nn = 0; nn < NN; nn++) {
      XX = cmplx(0.0,0.0);
      for (int kk = 0; kk < NN; kk++) {
        prom = cmplx (cos (kk*2*M_PI*nn/NN), -sin (kk*2*M_PI*nn/NN));
        XX += *(Input + kk) * prom;
      }
      *(Output + nn) = 2.0*XX/cmplx (NN,NN);   // Factor 2.0 is empirical, without it the inverse function is half of magnitude.
    }
  }
 
  public : void InverseFourier (cmplx *const Input, cmplx *const &Output, int NN)
  {
    cmplx xx,prom;      // WORKS VERY WELL June 1, 2017
    for (int nn = 0; nn < NN; nn++) {
      xx = cmplx(0.0,0.0);
      for (int kk = 0; kk < NN; kk++) {
        prom = cmplx (cos (kk*2*M_PI*nn/NN), sin (kk*2*M_PI*nn/NN));
        xx += *(Input + kk) * prom;
      }
      *(Output + nn) = xx;
    }
  }
};   // homeFourier


When I run my input function through this class I get perfect restoration. When I run the same input function though your FFT I get perfect restoration as well HOWEVER, your coefficients (the result of the forward FFT) do not match my coefficients after my forward DFT. Not only the absolute values of the corresponding coefficients are different, the signs are different as well. I need MY coefficients but your speed, if it makes sense.

How can I achieve this?

The same happens when I do Librow routines. I haven't tried your DFT yet.

Thanks, - Alex
Last edited on
lastchance (6980)
On your line 12 your complex number in the normalisation appears to be N+Ni; it should be just the real number N. I don't recognise your "empirical" factor of 2.0 in that line either.

Both I, and the Librow routines put the division by N in the inverse transform, rather than the forward one, although you are obviously entitled to split it between either, particularly if you are trying to approximate the integral transform.
alexBB (131)
@lastchance thank you. I've made corrections in the forward Fourier routine.

1
2
3
4
5
6
7
8
9
10
11
12
 
public : void ForwardFourier (cmplx *const Input, cmplx *const &Output, int NN)
  {
    cmplx XX,prom;   // WORKS VERY WELL June 1, 2017
    for (int nn = 0; nn < NN; nn++) {
      XX = cmplx(0.0,0.0);
      for (int kk = 0; kk < NN; kk++) {
        prom = cmplx (cos (kk*2*M_PI*nn/NN), -sin (kk*2*M_PI*nn/NN));
        XX += *(Input + kk) * prom;
      }
      *(Output + nn) = XX/cmplx (NN,0);
    }
  }


I ran three packages for comparison your, mine and Librow's. They all work well but the Fourier coefficients are all different. I wonder how it could be.

Thanks, - AlexBB
Last edited on
lastchance (6980)
When I ran my original test case I checked against Librow's and got exactly the same coefficients. I've also checked against various samples on the web and got the answers cited in those.

Your DFT has the factor 1/N in the forward transform; both mine and Librow's put it in the inverse transform (see equations (3) and (13) in the Librow link that you gave at the start of the thread) - so there is bound to be a difference between yours and the other two there.

Be aware that Librow use their own non-standard complex class. I've just used the one from the C++ standard.

I suspect that you are simply calling the various routines incorrectly. Possibly you are not sure how you are splitting the 1/N factor between forward and inverse transforms. Also, in your code I am not a good enough programmer to know if function parameters written as, e.g., cmplx *const &Output are correct: this is a very unfamiliar coding style to me; you shouldn't need to put a & for an array parameter as you are already passing a pointer to the start of the array and you aren't changing the memory address associated with that pointer ... just the content of that memory.


A comment on your own code. Writing
prom = cmplx (cos (kk*2*M_PI*nn/NN), -sin (kk*2*M_PI*nn/NN))
is an extremely expensive way of calculating
e-2 PI i k n / N
This is the same as
wk
where
w = e-2 PI i n / N
w only needs to be calculated once, with the kth power from successive complex multiplies.


Finally, note that an FFT should give exactly the same end result as a standard DFT ... just a lot faster. There are some irritating commercial packages (including, I believe, versions of matlab) that pad the original array with extra zeroes in order to allow it to do a supremely fast FFT. In my opinion - based on what it does to noisy lab data - this is wrong: it puts spurious extra frequency components in the signal and can lead to wrong conclusions about the frequencies you need to damp.
Last edited on
alexBB (131)
@lastchance, thank you. I will check everything out. Your replies are very helpful. - A.

P.S. I am tempted to mark it "solved," but I may have more questions in a couple of days, so I will wait.
Topic archived. No new replies allowed.