C++

Forum

 
How to combine fft and ifft function using class?

nurulhudaismail (43)
Hi, I am doing fast fourier transform(FFT) and inverse FFT (ifft) and I am using class in my code. Since I am dealing with Vivado HLS compiler, so I need to split my code into 3 files: header file, cpp file, and testbench file. How am I going to combine my FFT function and IFFT function into one main function? Thank you.

FFT file:
header file: 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 
//header.h
#include <iostream>
#include <complex>
#include <cmath>
#include <iterator>

using namespace std;

const double PI = 3.1415926536;

//template<typename R>
class Iter_T
{
	public:
		void fft(Iter_T a, Iter_T b, int log2n);

};


unsigned int bitReverse(unsigned int x, int log2n);



CPP file:
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
 
//fft.cpp
#include "lib.h" 


unsigned int bitReverse(unsigned int x, int log2n) {
	int n = 0;
	int mask = 0x1;
	for (int i=0; i < log2n; i++) {
		n <<= 1;
		n |= (x & 1);
		x >>= 1;
	}
	return n;
}


template<class Iter_T> 

void fft(Iter_T a, Iter_T b, int log2n)
{
	typedef typename iterator_traits<Iter_T>::value_type complex;
	const complex J(0, 1);
	int n = 1 << log2n;
	for (unsigned int i=0; i < n; ++i) {
		b[bitReverse(i, log2n)] = a[i];
	}
	for (int s = 1; s <= log2n; ++s) {
		int m = 1 << s;
		int m2 = m >> 1;
		complex w(1, 0);
		complex wm = exp(-J * (PI / m2));
		for (int j=0; j < m2; ++j) {
			for (int k=j; k < n; k += m) {
				complex t = w * b[k + m2];
				complex u = b[k];
				b[k] = u + t;
				b[k + m2] = u - t;
			}
			w *= wm;
		}
	}
}



testbench file:
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
 
//testbench.cpp
#include "lib.h"
                                                                                                              
template<class Iter_T> 
void fft(Iter_T a, Iter_T b, int log2n)
{
        typedef typename iterator_traits<Iter_T>::value_type complex;
        const complex J(0, 1);
        int n = 1 << log2n;
        for (unsigned int i=0; i < n; ++i) {
                b[bitReverse(i, log2n)] = a[i];
        }
        for (int s = 1; s <= log2n; ++s) {
                int m = 1 << s;
                int m2 = m >> 1;
                complex w(1, 0);
                complex wm = exp(-J * (PI / m2));
                for (int j=0; j < m2; ++j) {
                        for (int k=j; k < n; k += m) {
                                complex t = w * b[k + m2];
                                complex u = b[k];
                                b[k] = u + t;
                                b[k + m2] = u - t;
                        }
                        w *= wm;
                }
        }
}



int main() {

	typedef complex<double> cx;

	//cx a[] = { cx(1,0), cx(2,0), cx(3,0), cx(4,0), cx(5,0), cx(6,0), cx(7,0), cx(8,0) };
	cx a[] = { {64.0261,0.0},{5.96217,2.0254}, {0.0,0.0}, {5.96217,-2.0254} };
	cx b[8];

	fft(a, b, 3);

	for (int i=0; i<8; ++i) 
		cout << b[i] << "\n";
}


What I would like to add is the inverse fft (ifft) such that:

1
2
3
4
5
6
7
8
9
10
11
12
13
 
// Inverse Fast Fourier Transform
void iFFT( Complex ftilde[], Complex f[], int N )                    
{
   Complex* ftildeConjugate = new Complex[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] = conj( f[m] ) * factor;

   delete [] ftildeConjugate;
}
jonnin (11437)
sorry, carefully explain the requirement again. (Not the FFT part). You only ever have 'one main function' if you mean 'int main. That part is confusing.

Lets maybe ask a more basic question:
you run your program. what do you want it to DO, EXACTLY (gloss over the math. FFT gives me a headache, but I know what it means and does and don't really care. I want to know something like "print the FFT and iFFT for some problem" or words like that to explain what you want to accomplish.
And, to get to THAT point, are there any odd restrictions with your system? You can have multiple functions, right?

Last edited on
nurulhudaismail (43)
Sorry if I make you confuse. To be clear, I got the code to solve FFT from this forum by lastchance program as written below: 

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
100
101
102
103
 
#include <iostream>
#include <iomanip>
#include <complex>
using namespace std;

const double pi = 3.14159265358979323846;
using Complex = complex<double>;

void print(const char * prompt, Complex A[], int N);
void FFT(Complex f[], Complex ftilde[], int N, int step = 1);
void DFT(Complex f[], Complex ftilde[], int N, int step = 1);
void iFFT(Complex ftilde[], Complex f[], int N);

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

int main()
{
	Complex test[] = { 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 };
	int N = sizeof test / sizeof test[0];
	Complex * ft = new Complex[N];
	print("Original:", test, N);

	// Forward transform
	FFT(test, ft, N);   print("Transform:", ft, N);

	// Invert transform
	iFFT(ft, test, N);   print("Invert:", test, N);

	delete[] ft;
}


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

void print(const char * prompt, Complex A[], int N)
{
	cout << prompt << '\n' << fixed;
	for (int i = 0; i < N; i++) cout << A[i] << '\n';
}

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

void FFT(Complex f[], Complex ftilde[], int N, int step)           // Fast Fourier Transform
{
	if (N > 2 && N % 2 == 0)                                        // Recursion if a multiple of 2
	{
		int N2 = N >> 1;
		Complex* g = new Complex[N];

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

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

		delete[] g;
	}
	else                                                              // Otherwise just do a DFT
	{
		DFT(f, ftilde, N, step);
	}
}

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

void DFT(Complex f[], Complex ftilde[], int N, int step)           // Basic Discrete Fourier Transform
{
	Complex w = polar(1.0, -2.0 * pi / N);
	Complex wn = 1.0;
	for (int n = 0; n < N; n++)
	{
		if (n > 0) wn *= w;
		Complex wm = 1.0;
		ftilde[n] = 0.0;
		for (int m = 0, mpos = 0; m < N; m++, mpos += step)
		{
			if (m > 0) wm *= wn;
			ftilde[n] += f[mpos] * wm;
		}
	}
}

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

void iFFT(Complex ftilde[], Complex f[], int N)                    // Inverse Fast Fourier Transform
{
	Complex* ftildeConjugate = new Complex[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] = conj(f[m]) * factor;

	delete[] ftildeConjugate;
}


However, my system does not support recursion method. Finally, I found the FFT code which works with my hardware system. Unfortunately, it is not complete since it consists of FFT code only. What I would like to do is to apply the IFFT part as written above into the FFT code that I have in the last post. So that my code is complete (FFT and IFFT). My problem is I don't know how to integrate the IFFT code into FFT code. I hope it is more clear. Thank you.

'You can have multiple functions, right?'..Yes, I can have multiple functions.
Last edited on
jonnin (11437)
web search says this..


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 

typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;

void ifft(CArray& x)
{
    // conjugate the complex numbers
    x = x.apply(std::conj);
 
    // forward fft
    fft( x );
 
    // conjugate the complex numbers again
    x = x.apply(std::conj);
 
    // scale the numbers
    x /= x.size();
}


which assumes that FFT modifies X itself, so you may need to adjust that line.
nurulhudaismail (43)
Thank you for your quick response. I am so sorry, the system does not support valarray library. I have tried this code before, but it didnot work. That is why, I insist to use 

1
2
3
4
5
6
7
8
9
10
11
12
13
 
// Inverse Fast Fourier Transform
void iFFT( Complex ftilde[], Complex f[], int N )                    
{
   Complex* ftildeConjugate = new Complex[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] = conj( f[m] ) * factor;

   delete [] ftildeConjugate;
}


I just want to get help in term how to write the correct way to integrate with the exist FFT code.
Last edited on
lastchance (6980)
@Nurulhuda,
I don't understand your code, but if you want an FFT without any recursion then you can try the code below. For the forward transform it is a hybrid of what I understand of the FFT with reordering and the code that you have somehow got hold of. For the backward transform I have only made minor changes to my original code.

If that doesn't work on your picky system then please tell us exactly which fragments of the C++ language it can actually handle. You can use a typedef for an alias instead of the "using" statement. I don't intend to do any more refactoring.

Note that, in contrast to my original code here:
http://www.cplusplus.com/forum/general/263983/#msg1139656
this FFT with reordering will only work if N is a power of 2 and there is no simple "non-fast" discrete Fourier transform to compare it with.

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
 
#include <iostream>
#include <iomanip>
#include <complex>
using namespace std;

const double pi = 3.14159265358979323846;
using Complex = complex<double>;

void print( const char * prompt, Complex A[], int log2N );
void FFT ( Complex f[]     , Complex ftilde[], int log2N );
void iFFT( Complex ftilde[], Complex f[]     , int log2N );

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

int main()
{
   const int N = 8, log2N = 3;                                             // Hardwired for testing
   //Complex fN] = { 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 };
   Complex f[N] = { 1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 1.0, -1.0 };
   Complex ft[N];

   print( "Original:", f, log2N );

   // Forward transform
   FFT ( f, ft, log2N );   print( "\nTransform:", ft, log2N );

   // Invert transform
   iFFT( ft, f, log2N );   print( "\nInvert:",     f, log2N );
}


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

void print( const char * prompt, Complex A[], int log2N )
{
   int N = 1 << log2N;
   cout << prompt << '\n' << fixed;
   for ( int i = 0; i < N; i++ ) cout << A[i] << '\n';
}

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

void FFT( Complex f[], Complex ftilde[], int log2N )                 // Fast Fourier Transform
{
   int N = 1 << log2N;

   // Reorder
   for ( int i = 0; i < N; i++ )
   {
      int ii = 0, x = i;
      for (int j = 0; j < log2N; j++)
      {
         ii <<= 1;
         ii |= ( x & 1 );
         x >>= 1;
      }
      ftilde[ii] = f[i];
   }

   for ( int s = 1; s <= log2N; s++ )
   {
      int m = 1 << s;
      int m2 = m / 2;
      Complex w = 1.0;
      Complex wm = polar( 1.0, -pi / m2 );
      for ( int j = 0; j < m2; j++ )
      {
         for ( int k = j; k < N; k += m )
         {
            Complex t = w * ftilde[k+m2];
            Complex u =     ftilde[k   ];
            ftilde[k   ] = u + t;
            ftilde[k+m2] = u - t;
         }
         w *= wm;
      }
   }
}

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

void iFFT( Complex ftilde[], Complex f[], int log2N )                // Inverse Fast Fourier Transform
{
   int N = 1 << log2N;

   for ( int m = 0; m < N; m++ ) ftilde[m] = conj( ftilde[m] );      // Apply conjugate (reversed below)

   FFT( ftilde, f, log2N );

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

   for ( int m = 0; m < N; m++ ) ftilde[m] = conj( ftilde[m] );      // Only necessary to reinstate ftilde
}

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

Original:
(1.000000,0.000000)
(2.000000,0.000000)
(3.000000,0.000000)
(4.000000,0.000000)
(5.000000,0.000000)
(4.000000,0.000000)
(1.000000,0.000000)
(-1.000000,0.000000)

Transform:
(19.000000,0.000000)
(-8.949747,-4.121320)
(2.000000,-3.000000)
(0.949747,-0.121320)
(1.000000,0.000000)
(0.949747,0.121320)
(2.000000,3.000000)
(-8.949747,4.121320)

Invert:
(1.000000,-0.000000)
(2.000000,0.000000)
(3.000000,-0.000000)
(4.000000,-0.000000)
(5.000000,-0.000000)
(4.000000,0.000000)
(1.000000,0.000000)
(-1.000000,-0.000000)


Last edited on
Topic archived. No new replies allowed.