replace valarray library with orthodox c++ version

[nurulhudaismail]

Hi, I would like to do Fast Fourier Transform (FFT) on the hardware. However,I need to use classical c++ code which is more practical to the hardware. Therefore, my question is how to replace 'valarray' library with something which is more classic and practical. The code is as below. Thank you.

#include <iostream>
#include <iomanip>
#include <complex>
#include <valarray>

const double PI = 3.141592653589793238460;

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

void fft(CArray &x)
{
    // DFT
    unsigned int N = x.size(), k = N, n;
    double thetaT = 3.14159265358979323846264338328 / N;
    Complex phiT = Complex(cos(thetaT), -sin(thetaT)), T;
    while (k > 1)
    {
        n = k;
        k >>= 1;
        phiT = phiT * phiT;
        T = 1.0L;
        for (unsigned int l = 0; l < k; l++)
        {
            for (unsigned int a = l; a < N; a += n)
            {
                unsigned int b = a + k;
                Complex t = x[a] - x[b];
                x[a] += x[b];
                x[b] = t * T;
            }
            T *= phiT;
        }
    }

	// Decimate
	unsigned int m = (unsigned int)log2(N);
	for (unsigned int a = 0; a < N; a++)
	{
		unsigned int b = a;
		// Reverse bits
		b = (((b & 0xaaaaaaaa) >> 1) | ((b & 0x55555555) << 1));
		b = (((b & 0xcccccccc) >> 2) | ((b & 0x33333333) << 2));
		b = (((b & 0xf0f0f0f0) >> 4) | ((b & 0x0f0f0f0f) << 4));
		b = (((b & 0xff00ff00) >> 8) | ((b & 0x00ff00ff) << 8));
		b = ((b >> 16) | (b << 16)) >> (32 - m);
		if (b > a)
		{
			Complex t = x[a];
			x[a] = x[b];
			x[b] = t;
		}
	}
}

int main()
{
	const Complex test1[] = { 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 };

// forward fft
	CArray data1(test1, 8);
	fft(data1);

	std::cout << "fft" << std::endl;
	for (int i = 0; i < 8; ++i)
	{
		std::cout << data1[i] << std::endl;
	}

std::cout << "\n";
	system("pause");
	return 0;
}

[lastchance]

I have absolutely no idea why you would want to switch away from valarray here. There's nothing unorthodox about valarray: it has been in the c++ standard from the start. The valarray library was designed to be very fast and efficient for purely numerical arrays and a lot of scientists and engineers use it. I don't know what you mean by "classical c++", but if you want to switch your containers for vector instead of valarray then you can do just:

Change line 4 to
#include <vector>

Change line 9 to
typedef std::vector<Complex> CArray;

Change line 61 to
CArray data1(test1, test1+8);



Please explain what you mean by "more practical to the hardware". Exactly what limitations is it placing on you?

From "classical c++" (did the ancient greeks and romans have c++?) I suspect that you are after something written in C, not C++.

[nurulhudaismail]

Thank you for your comment. I would like to simulate my c++ code in Xilinx Vivado HLS which is used to estimate a latency, throughput, estimed resources (FlipFlop, LUT, BRAM) etc. Unfortunately, Vivado HLS can only understand the basic c++ instructions as a low-level language and the HLS libraries are limited. From my understanding, 'valarray' is not included in the HLS libraries. That is why I insist to replace 'valarray' with other basic libraries.

[lastchance]

Does your simulator have any support for a complex number class? Looking at its limited web documentation, does it actually have support for classes at all? If it doesn't, then C++ would be a misnomer.

I think you need to find an FFT library in C for your restricted environment.

[nurulhudaismail]

One more thing, may I know how to write the basic c++ to represent line 14 which relates to size?

Thank you.

[jonnin]

if you do not have valarray, do you even have vectors?
the C way is to use an array or a pointer and track size yourself in a variable, or lump up a struct with the size and data together.

the way to write line 14 is this:


void fft(CArray &x, unsigned int size)
{
// DFT
unsigned int N = size, k = N, n;

or, just make the parameter name N and bypass copying it an extra time:

void fft(CArray &x, unsigned int N)
{
// DFT
unsigned int k = N, n;

I did not see any valarray tricks here (eg slice). Watch out for those though; you will need to do some work to re-create that. Simple storage/array like use is easy to re-do.
its likely your machine has a CPU reverse bit instruction that could tear up that big pile of junk in the middle. Consider calling that directly.

[nurulhudaismail]

I have tried to change as suggested by it does not work. I wrote line 8 and line 9 as below:

[code]
typedef std::complex<double> Complex;
typedef std::vector<Complex> CArray;
[\code].

However, the system cannot recognise 'cArray'. Might be I wrote it wrong. With a very minimum knowledge in c++, I think another option is to take out the fft function. That means I put all the codes under 'main' function. Am I correct? If yes, how it should be written?

Thank you.

[lastchance]

nurulhudaismail wrote:

However, the system cannot recognise 'cArray'

There's a world of a difference between cArray and CArray - c++ is case-sensitive. Check your typing.

nurulhudaismail wrote:

That means I put all the codes under 'main' function. Am I correct?

If it didn't work before then it won't be any better doing that.

[nurulhudaismail]

Is that mean that pointer and addressing such as *x and &x must be written as a function?

[jonnin]

if you must you can use a pure C array:
complex<double> derp[1000]; // <--- a C style array. This is if vector isnt working or available.

if vector is working, you don't need to pass size around as I was showing. That is only needed if size isnt built in, like the above C style array you need that 1000 to be available or if you used 23 of the slots of the 1000, you need that 23 available, etc.

I highly recommend you don't use the same word, with only a caps letter apart, for anything ever. Its legal and you can do it, but youll be sorry when you typo one of them. Its not as catastrophic when its a type Type; declare, but the habit will get you in trouble someday.

[nurulhudaismail]

Hi,

With <valarray> library, I use 'apply' as my instruction such as x = x.apply(std::conj). May I know what should I replace with if I use <vector> library?

// inverse fft (in-place)
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();
}

Thank you.

[lastchance]

You could use std::for_each() or you write explicit loops.

Are you saying that this strange implementation has std::vector but not std::valarray?

[JLBorges]

#include <iostream>
#include <complex>
#include <algorithm>
#include <iterator>
#include <vector>
#include <list>
#include <deque>

template < typename SEQUENCE, typename UNARY_FN >
void apply_fn( SEQUENCE& seq, UNARY_FN fn )
{ std::transform( std::begin(seq), std::end(seq), std::begin(seq), fn ) ; }

int main()
{
    const auto print = []( const auto& seq )
    { for( const auto& v : seq ) std::cout << v << "  " ; std::cout << '\n' ; };

    std::complex<double> array[] { {1.1,2.2}, {3.3,4.4}, {5.5,6.6}, {7.7,8.8} } ;
    print(array) ;

    apply_fn( array, std::log<double> ) ;
    print(array) ;

    std::vector< std::complex<double> > vec( std::begin(array), std::end(array) ) ;
    apply_fn( vec, std::exp<double> ) ;
    print(vec) ;

    std::list< std::complex<double> > lst( std::begin(vec), std::end(vec) ) ;
    apply_fn( lst, []( const auto& v ) { return v*9.9 ; } ) ;
    print(lst) ;

    std::deque< std::complex<double> > deq( std::begin(lst), std::end(lst) ) ;
    apply_fn( deq, []( const auto& v ) { return v/9.9 ; } ) ;
    print(deq) ;
}

http://coliru.stacked-crooked.com/a/6682242782196e4d

[lastchance]

#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;
}

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

Original:
(1.000000,0.000000)
(1.000000,0.000000)
(1.000000,0.000000)
(1.000000,0.000000)
(0.000000,0.000000)
(0.000000,0.000000)
(0.000000,0.000000)
(0.000000,0.000000)
Transform:
(4.000000,0.000000)
(1.000000,-2.414214)
(0.000000,0.000000)
(1.000000,-0.414214)
(0.000000,0.000000)
(1.000000,0.414214)
(0.000000,0.000000)
(1.000000,2.414214)
Invert:
(1.000000,-0.000000)
(1.000000,0.000000)
(1.000000,-0.000000)
(1.000000,-0.000000)
(0.000000,-0.000000)
(0.000000,-0.000000)
(0.000000,0.000000)
(0.000000,0.000000)

[nurulhudaismail]

Thank you. Really appreciate it. Unfortunately, I just know that <vector> and <valarray> are not supported by Vivado. It can be simulated but cannot be synthesized.

May I know how to do FFT in a very very basic c++? Thank you.

[keskiverto]

It does sound like your "very very basic c++" means "C with (barely any) classes".
Which version of C++ Standard does "Vivado" understand? I bet it isn't C++17.

There are no vectors nor valarrays in lastchance's program.

[jonnin]

May I know how to do FFT in a very very basic c++? Thank you.

you do it the same way.
there are literally hundreds of old C++ or (age unimportant) C programs that do FFT out there on the web; maybe a pure C one will work on your compiler (?).
Its going to involve either arrays or pointer-arrays as your container of choice, and the math is the same as its always been.

I don't recall if C has complex; if not, you may have to roll your own (if you end up rewriting it) with a struct or class. Ideally you would have operator overloading so you can make it look like a double or int without any extras but however you handle it, its doable. See if your compiler has a complex?