By using this site, you agree to our updated Privacy Policy and our Terms of Use. Manage your Cookies Settings.
435,395 Members | 2,537 Online
Bytes IT Community
+ Ask a Question
Need help? Post your question and get tips & solutions from a community of 435,395 IT Pros & Developers. It's quick & easy.

C++ - STL usage 'Advice'

P: n/a

I'm perusing a slide with roughly 12 bullets spread across 3 pages.
Each bullet reflects 'advice'. I'm ok with all but 1 bullet, more
specifically the bullet that states:
" Avoid the STL unless absolutely necessary "

Now, I'm not acclimated with much C++ history, but something tells me
this is akin to trepidation that surrounded C++ during it's inception?
IOW, is this 'old school' rhetoric or ..... How do you refute this
argument?

I persued the web and there was discussion about code bloat (pre 2000
posts) but I'm not quite fure if I'm following that. Oh well, I
thought I'd inquire within. I'm a fan of lang.c++/learn and moderated
c++ and rest assured the vast majority of folks on here know what the
hell they're talking about so somebody will set me straight.

Granted, my background is signal proceesing, I've enjoyed the language
to include the STL and all it's glory (the Big O notation still
confuses me) but I'm tempted to go to the upcoming Vegas conference on
the language and ask Bjarne, Scott, Herb etc. this question :)
// formerly ma******@pegasus.cc.ucf.edu

Oct 11 '05 #1
Share this Question
Share on Google+
11 Replies


P: n/a
On 11 Oct 2005 05:12:37 -0700, ma******@gmail.com wrote:
I'm perusing a slide with roughly 12 bullets spread across 3 pages.
Each bullet reflects 'advice'. I'm ok with all but 1 bullet, more
specifically the bullet that states:
" Avoid the STL unless absolutely necessary "
That prescription sounds biased against the use of STL. Perhaps a more
appropriate warning should be: "Be aware of the cost and consequences of using
STL before using it".
Now, I'm not acclimated with much C++ history, but something tells me
this is akin to trepidation that surrounded C++ during it's inception?
IOW, is this 'old school' rhetoric or ..... How do you refute this
argument?


You're right about the "old school" bias in many ways. However, there are
valid reasons not to use STL, especially if you live in the embedded world
where dynamic memory usage is frowned upon, due to the significant risk of
heap fragmentation. I can't use STL containers in my embedded project because
of that very reason.

However, STL (if used properly) gives you a wonderful toolkit, so that you
don't have to reinvent the wheel each time. It's also very efficient--many
implementations are as efficient as hand-written code once it gets through the
optimizer--so there are many compelling reasons to use it. Plus, it's part of
the language specification, so you know it's going to be portable.

There are other practical reasons not to use STL, e.g. when your program
already uses a similar class library (like MFC's CString and xxxArray classes)
extensively.

Using any library or feature is a cost-benefit tradeoff. Take a look at your
project and evaluate each toolkit against your needs. Don't dismiss anything
out of hand.

-dr
Oct 11 '05 #2

P: n/a
Dave Rahardja wrote:
On 11 Oct 2005 05:12:37 -0700, ma******@gmail.com wrote:
Now, I'm not acclimated with much C++ history, but something tells me
this is akin to trepidation that surrounded C++ during it's inception?
IOW, is this 'old school' rhetoric or ..... How do you refute this
argument?


You're right about the "old school" bias in many ways. However, there are
valid reasons not to use STL, especially if you live in the embedded world
where dynamic memory usage is frowned upon, due to the significant risk of
heap fragmentation. I can't use STL containers in my embedded project because
of that very reason.


That's not a case for avoiding the Standard Library - it's a case for
avoiding dynamically sized data structures altogether. But what's
stopping you from calling reserve() during the initialization phase?
This may vary by implementation, but I think that effectively gives
your container a pool of memory to play with - saving you from
fragmenting the rest of the heap.

The whole point of C++ and the Standard Library is to be able to use
them in a variety of ways, including for embedded systems. It's just a
matter of figuring out the right way to use them in your situation.
Maybe you have to reserve memory in advance. Maybe you have to write
your own allocator (I've seen that done before). Don't give up on the
Standard Library just because someone told you to.

Kristo

Oct 11 '05 #3

P: n/a
Dave Rahardja wrote:
On 11 Oct 2005 05:12:37 -0700, ma******@gmail.com wrote:

I'm perusing a slide with roughly 12 bullets spread across 3 pages.
Each bullet reflects 'advice'. I'm ok with all but 1 bullet, more
specifically the bullet that states:
" Avoid the STL unless absolutely necessary "

That prescription sounds biased against the use of STL. Perhaps a more
appropriate warning should be: "Be aware of the cost and consequences of using
STL before using it".
[...]


Still sounds like STL (I presume the Standard library is meant) is somehow
a bad thing. I'd probably put it this way: "Weigh pros and cons of using
the Standard library mechanisms". Of course, being afraid of it should be
listed among "cons".

V
Oct 11 '05 #4

P: n/a
ma******@gmail.com wrote:
I'm perusing a slide with roughly 12 bullets spread across 3 pages.
Each bullet reflects 'advice'. I'm ok with all but 1 bullet, more
specifically the bullet that states:
" Avoid the STL unless absolutely necessary "

[snip]

It is a strangely worded directive. I'd say it's rarely
"absolutely necessary" to use the standard template
containers. It is often highly desirable, easier, faster
to design and code, easier to maintain, easier to document,
easier to test, etc. But I'd doubt that it is "absolutely
necessary" to use the STL. You can nearly always do the
same job "by hand" (by self-created code for example.)
It may be much more work to do it that way.

There may be important context missing here. Later (in the
part I snipped) you mention signal processing. It is just
possible that the intent of this directive is to achieve
something else that is not directly and trivially obvious
from the raw text. For example: If you were doing a system
that was going to be embedded in some device, it might make
sense. You might want to do things that made the smallest
possible stack size, just as an example. Or that used the
least possible amount of read/write memory. In many such
embedded systems there is only a very small amount of RAM,
and often that has to be shared among a lot of activity.
It may be difficult to have this level of minute detail
control if you pass off a lot of overhead to the standard
lib things.

So, by all means take into account the advice others have given.
But check with the powers-that-be in your project and try to
get them to tell you why this advice is there. It may be that
somebody got burnt by an old version of the standard lib that
did things poorly. Or there may in fact be a good reason for it.
Socks

Oct 11 '05 #5

P: n/a
ma******@gmail.com wrote:

How do you refute this
argument?


It's not an argument. It's a command. The way to refute it is to ignore
it. If whoever wrote that slide tells you why they wrote that particular
piece of [bad] advice then you can respond to their argument. Without
reasons, though, it's just noise.

--

Pete Becker
Dinkumware, Ltd. (http://www.dinkumware.com)
Oct 11 '05 #6

P: n/a
ma******@gmail.com wrote:
I'm perusing a slide with roughly 12 bullets spread across 3 pages.
Each bullet reflects 'advice'. I'm ok with all but 1 bullet, more
specifically the bullet that states:
" Avoid the STL unless absolutely necessary "

Now, I'm not acclimated with much C++ history, but something tells me
this is akin to trepidation that surrounded C++ during it's inception?
IOW, is this 'old school' rhetoric or ..... How do you refute this
argument?


What is the argument? In fact, what is the advice? As I understand it,
the advice is not to use the STL because it is outdated. Programs
should use the container and iterator types as defined in the Standard
Library in order to be current, to be portable and to be assured of
timely bug fixes if needed.

If the advice really is to avoid the STL and STL-derived libraries,
then what is the proposed alternative? C arrays? Everyone to write
their own container classes? In what ways are either of those two
alternatives more attractive than using an existing library?

Greg

Oct 11 '05 #7

P: n/a
[Dave Rahardja]
|| However, there are
|| valid reasons not to use STL, especially if you live in the embedded
world
|| where dynamic memory usage is frowned upon, due to the significant
risk of
|| heap fragmentation. I can't use STL containers in my embedded
project because
|| of that very reason.

[Kristo]
|| But what's stopping you from calling reserve() during the
initialization phase?
In the academic world. Advisors depend on grants, grants makes the
university look 'good'; the graduate students get to deal with the
'science project' received from the outside world; this make the
advsior look good, the student (us) get a pat on the back and our
Phd/Masters becomes more appealing. The folks in the embedded wold
are the ones who generally claim limited 'STL' on a continuum. This
brings about slides with the bullet in question. So now consider an
approach to a Discrete Fourier Transform (DFT):

// the 'C' way.
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <numeric>
#include <functional>
#include <vector>
#include <algorithm>
using namespace std;

enum {
TRUE=1,
FALSE=0
};

/*
Direct fourier transform
dir = 1 gives forward transform
dir = -1 gives reverse transform
*/
int DFT(int dir,int m,double *x1,double *y1)
{
long i,k;
double arg;
double cosarg,sinarg;
double *x2=NULL,*y2=NULL;

x2 = (double *)malloc(m*sizeof(double));
y2 = (double *)malloc(m*sizeof(double));
if (x2 == NULL || y2 == NULL)
return(FALSE);

for (i=0;i<m;i++) {
x2[i] = 0;
y2[i] = 0;
arg = - dir * 2.0 * 3.141592654 * (double)i / (double)m;
for (k=0;k<m;k++) {
cosarg = cos(k * arg);
sinarg = sin(k * arg);
x2[i] += (x1[k] * cosarg - y1[k] * sinarg);
y2[i] += (x1[k] * sinarg + y1[k] * cosarg);
}
}

/* Copy the data back */
if (dir == 1) {
for (i=0;i<m;i++) {
x1[i] = x2[i] / (double)m;
y1[i] = y2[i] / (double)m;
}
} else {
for (i=0;i<m;i++) {
x1[i] = x2[i];
y1[i] = y2[i];
}
}

free(x2);
free(y2);
return(TRUE);
}

I'm now a fan of the STL, I'm now programmed to think STL. STL
implementations of FFT/IFFTs etc are - nonexistent (I have yet to find
one). Having said, lets experiment with an STL version of the DFT
above:

class MinusWithCoeff {
public:
MinusWithCoeff(double arg_) : i(0), arg(arg_) {}
pair<double, double> operator()(double x, double y)
{
double cosarg = cos(i * arg);
double sinarg = sin(i * arg);
++i;
pair<double, double> res;
res.first = x*cosarg-y*sinarg;
res.second = x*sinarg + y*cosarg;
return res;
}
private:
long i;
double arg;
};

template <typename T>
class pair_plus {
public:
pair<T, T> operator()(const pair<T, T>& lh, const pair<T, T>& rh)
{
pair<T, T> res;
res.first=lh.first+rh.first;
res.second=lh.second+rh.second;
return res;
}
};

/*
Direct fourier transform
dir = 1 gives forward transform
dir = -1 gives reverse transform
*/
int DFT_CPP(int dir, int m, double *x1, double *y1)
{
long i, k;
double arg;
vector<double> x2(m), y2(m);

for (i=0; i<m; i++) {
arg = - dir * 2.0 * 3.141592654 * double(i) / double(m);
pair<double, double> res = inner_product(x1, x1+m, y1, pair<double,
double>(0., 0.),
pair_plus<double>(), MinusWithCoeff(arg));
x2[i]=res.first;
y2[i]=res.second;
}

/* Copy the data back */
if (dir == 1) {
for (i=0;i<m;i++) {
x1[i] = x2[i] / double(m);
y1[i] = y2[i] / double(m);
}
} else {
copy(x2.begin(), x2.end(), x1);
copy(y2.begin(), y2.end(), y1);
}
return(TRUE);
}

class TestUnit {
public:
void test_new_dft(void)
{
int m=2;
for (int dir=-1; dir<=1; dir+=2) {
vector<double> x1(m, 17.), y1(m, 18.);
vector<double> x2(m, 17.), y2(m, 18.);
if (!DFT_CPP(dir, m, &x1[0], &y1[0]) ||
!DFT(dir, m, &x2[0], &y2[0]) ||
!equal(x1, x2) || !equal(y1, y2))
cout<<"Test failed with dir="<<dir<<endl;
else
cout<<"Test passed with dir="<<dir<<endl;
}
}
private:
bool equal(const vector<double>& rh, const vector<double>& lh) {
static const double eps=1e-14;
vector<double>::const_iterator p, q;
for (p=rh.begin(), q=lh.begin(); p!=rh.end(), q!=lh.end(); ++p,
++q)
if (fabs(*p-*q)>eps)
return false;
return p==rh.end() && q==lh.end();
}
};

int main()
{
TestUnit().test_new_dft();
return EXIT_SUCCESS;
}
Granted, I haven't benchmarked this, I suspect it's slow compared to
the C style approach but from the looks of it, one's head is liable to
spin. Of coruse BiG O factors in here and might might make things slow
but nowithstanding executing speed, I don't see the code bloat.
Complexity, yes!!

Oct 11 '05 #8

P: n/a
ma******@gmail.com wrote:
[Dave Rahardja]
|| However, there are
|| valid reasons not to use STL, especially if you live in the embedded
world
|| where dynamic memory usage is frowned upon, due to the significant
risk of
|| heap fragmentation. I can't use STL containers in my embedded
project because
|| of that very reason.

[Kristo]
|| But what's stopping you from calling reserve() during the
initialization phase?

[snip] So now consider an
approach to a Discrete Fourier Transform (DFT):
[the 'C' way snipped]
I'm now a fan of the STL, I'm now programmed to think STL. STL
implementations of FFT/IFFTs etc are - nonexistent (I have yet to find
one). Having said, lets experiment with an STL version of the DFT
above:

class MinusWithCoeff {
public:
MinusWithCoeff(double arg_) : i(0), arg(arg_) {}
pair<double, double> operator()(double x, double y)
{
double cosarg = cos(i * arg);
double sinarg = sin(i * arg);
++i;
pair<double, double> res;
res.first = x*cosarg-y*sinarg;
res.second = x*sinarg + y*cosarg;
return res;
}
private:
long i;
double arg;
};

template <typename T>
class pair_plus {
public:
pair<T, T> operator()(const pair<T, T>& lh, const pair<T, T>& rh)
{
pair<T, T> res;
res.first=lh.first+rh.first;
res.second=lh.second+rh.second;
return res;
}
};

/*
Direct fourier transform
dir = 1 gives forward transform
dir = -1 gives reverse transform
*/
int DFT_CPP(int dir, int m, double *x1, double *y1)
{
long i, k;
double arg;
vector<double> x2(m), y2(m);

for (i=0; i<m; i++) {
arg = - dir * 2.0 * 3.141592654 * double(i) / double(m);
pair<double, double> res = inner_product(x1, x1+m, y1, pair<double,
double>(0., 0.),
pair_plus<double>(), MinusWithCoeff(arg));
x2[i]=res.first;
y2[i]=res.second;
}

/* Copy the data back */
if (dir == 1) {
for (i=0;i<m;i++) {
x1[i] = x2[i] / double(m);
y1[i] = y2[i] / double(m);
}
} else {
copy(x2.begin(), x2.end(), x1);
copy(y2.begin(), y2.end(), y1);
}
return(TRUE);
}
[class TestUnit snipped]
Granted, I haven't benchmarked this, I suspect it's slow compared to
the C style approach but from the looks of it, one's head is liable to
spin. Of coruse BiG O factors in here and might might make things slow
but nowithstanding executing speed, I don't see the code bloat.
Complexity, yes!!


You make it look so unnatural. Consider the following two versions:

A simple straight forward adaptation of your code just movin to complex
numbers:

#include <cstdlib>
#include <cmath>
#include <iostream>
#include <numeric>
#include <functional>
#include <vector>
#include <algorithm>
#include <complex>

typedef double real_number;
typedef std::complex< real_number > complex_number;
typedef std::vector< complex_number > complex_vector;

void DFT( bool do_forward, complex_vector & z ) {
typedef complex_vector::size_type size_type;

const real_number sign ( do_forward ? 1.0 : -1.0 );
const size_type length = z.size();

// allocate and initialize result
complex_vector result ( length, 0.0 );

// do matrix multiplication
for ( size_type i = 0; i < length; ++ i ) {
real_number arg =
- sign * 2.0 * 3.141592654 *
static_cast<real_number>(i)/ static_cast<real_number>( length );
for ( size_type k = 0; k < length; ++ k ) {
complex_number coefficient_ik
( std::cos( k * arg ), std::sin( k * arg ) );
result[i] += z[k] * coefficient_ik;
}
}

// rescale if forward
if ( do_forward ) {
for ( size_type i = 0; i < length; ++ i ) {
result[i] /= static_cast<real_number>( length );
}
}

// copy data back (swap is fast!)
z.swap( result );
}
A version where the coefficients of the transform are realized by fake
vectors:
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <numeric>
#include <functional>
#include <vector>
#include <algorithm>
#include <complex>
#include <iterator>

typedef double real_number;
typedef std::complex< real_number > complex_number;
typedef std::vector< complex_number > complex_vector;

class FFT_Pseudoiterator
: public std::iterator< std::input_iterator_tag,
complex_number,
std::ptrdiff_t,
complex_number const *,
complex_number const &>
{

typedef complex_vector::size_type size_type;

real_number arg;
size_type k;

public:

FFT_Pseudoiterator ( real_number sign, size_type i, size_type length )
: arg ( - sign * 2.0 * 3.141592654 * static_cast<real_number>(i)
/ static_cast<real_number>( length ) )
, k ( 0 )
{}

complex_number operator* ( void ) const {
return ( complex_number ( std::cos( k * arg ), std::sin( k * arg ) ) );
}

FFT_Pseudoiterator & operator++ ( void ) {
++ k;
return( *this );
}

FFT_Pseudoiterator operator++ ( int ) {
FFT_Pseudoiterator dummy ( *this );
++ k;
return( dummy );
}

};
void DFT( bool do_forward, complex_vector & z ) {
typedef complex_vector::size_type size_type;

const real_number sign ( do_forward ? 1.0 : -1.0 );
const size_type length = z.size();

// allocate result
complex_vector result;
result.reserve( length );

// do matrix multiplication
for ( size_type i = 0; i < length; ++ i ) {
FFT_Pseudoiterator coefficient_iter ( sign, i, length );
result.push_back( std::inner_product( z.begin(), z.end(),
coefficient_iter,
complex_number(0.0) ) );
}

// rescale if forward
if ( do_forward ) {
for ( size_type i = 0; i < length; ++ i ) {
result[i] /= static_cast<real_number>( length );
}
}

// copy data back (swap is fast!)
z.swap( result );
}
int main ( void ) {
complex_vector z;
z.push_back( complex_number( 0.5, 2 ) );
z.push_back( complex_number( 0.7, 1.4 ) );
z.push_back( complex_number( 0.9, -2.1 ) );
z.push_back( complex_number( 0.9, -2.1 ) );
z.push_back( complex_number( -1.2, 1.1 ) );

DFT( true, z );
std::copy( z.begin(), z.end(),
std::ostream_iterator< complex_number >( std::cout, "\n" ) );
DFT( false, z );
std::copy( z.begin(), z.end(),
std::ostream_iterator< complex_number >( std::cout, "\n" ) );
}
I think the C++ versions are much more easy to understand -- after all you
can see the complex arithmetic. Also, I would be surprised if they were
actually slower than the C version. After all this is a pretty straight
forward translation.

Best

Kai-Uwe Bux
Oct 12 '05 #9

P: n/a
|| I think the C++ versions are much more easy to understand -- after
all you
|| can see the complex arithmetic. Also, I would be surprised if they
were
|| actually slower than the C version. After all this is a pretty
straight
|| forward translation.

Impressive. I suspect it's time to reivew std::complex and iterators
in depth. I ended up perusing those chapters (oso lightly) in
Josuttis. Do you have anything like this for FFT and IFFT ;) Makes me
think my custom versions (at the lab right now) are also unnatural.

Thanks a lot.

Oct 12 '05 #10

P: n/a
ma******@gmail.com wrote:
|| I think the C++ versions are much more easy to understand -- after
all you
|| can see the complex arithmetic. Also, I would be surprised if they
were
|| actually slower than the C version. After all this is a pretty
straight
|| forward translation.

Impressive. I suspect it's time to reivew std::complex and iterators
in depth. I ended up perusing those chapters (oso lightly) in
Josuttis. Do you have anything like this for FFT and IFFT ;) Makes me
think my custom versions (at the lab right now) are also unnatural.


I don't have the faintest idea about FFT or IFFT -- I do not even know what
IFFT means. However, if you post some code, I will have a look.
Best

Kai-Uwe Bux

Oct 12 '05 #11

P: n/a
ma******@gmail.com wrote:
Do you have anything like this for FFT and IFFT ;)

So, I read up a little bit on this FFT thing. Seems that it is an easy
divide and conquer idea -- at least for powers of 2 (I did not bother
thinking hard about other primes). So, here is a first draft of how one
could go about this.

// FFT
// ===
/*
| Discrete Fourier Transform can be regarded as evaluating a
| polynomial of degree N-1 on the powers
| omega^0, omega, omega^2, ..., omega^(N-1) where omega is a the
| Nth root of unity.
|
| Given a polynomial of even degree
|
| p(t) = a_0 + a_1 t + a_2 t^2 + ...
|
| we find:
|
| p(t) = a_0 + a_2 t^2 + a_4 t^4 + ... + t( a_1 + a_3 t^2 + ... )
| = q_e( t^2 ) + t q_o ( t^2 )
|
| where q_e and q_o are polynomials formed with the even and odd
| coefficients of p(t). Thus, we get
|
| p(1) = q_e(1) + omega^0 q_o(1)
| p(omega) = q_e(omega^2) + omega^1 q_o(omega^2)
| p(omega^2) = q_e(omega^4) + omega^2 q_o(omega^4)
| ...
|
| Note how on the rhs the Fouriertransforms of q_e and q_o appear. Thus:
*/

#include <cstdlib>
#include <cmath>
#include <iostream>
#include <numeric>
#include <functional>
#include <vector>
#include <algorithm>
#include <complex>
#include <iterator>

namespace fft {

typedef double real_number;
typedef std::complex< real_number > complex_number;
typedef std::vector< complex_number > complex_vector;
typedef complex_vector::size_type size_type;

complex_number unit_complex ( real_number arg ) {
return complex_number ( std::cos( arg ), std::sin( arg ) );
}

real_number rational_arg ( size_type i, size_type m ) {
return( 2.0 * 3.141592654 * static_cast<real_number>(i)
/ static_cast<real_number>(m) );
}

bool is_odd ( size_type m ) {
return( m % 2 != 0 );
}

class DFT_Pseudoiterator
: public std::iterator< std::input_iterator_tag,
complex_number,
std::ptrdiff_t,
complex_number const *,
complex_number const &> {

complex_number omega_to_the_i;
complex_number omega_to_the_ik;
size_type k;

public:

DFT_Pseudoiterator ( real_number sign,
size_type i,
size_type length,
bool past_end = false )
: omega_to_the_i ( unit_complex ( - sign * rational_arg( i,
length ) ) )
, omega_to_the_ik ( 1.0 )
, k ( past_end ? length : 0 )
{}

complex_number operator* ( void ) const {
return ( omega_to_the_ik );
}

DFT_Pseudoiterator & operator++ ( void ) {
omega_to_the_ik *= omega_to_the_i;
++ k;
return( *this );
}

DFT_Pseudoiterator operator++ ( int ) {
DFT_Pseudoiterator dummy ( *this );
++ k;
return( dummy );
}

bool operator== ( DFT_Pseudoiterator const & other ) const {
return( ( this->omega_to_the_i == other.omega_to_the_i )
&&
( this->k == other.k ) );
}

bool operator!= ( DFT_Pseudoiterator const & other ) const {
return( ! ( *this == other ) );
}

}; // DFT_Pseudoiterator

void DFT( bool do_forward, complex_vector & z ) {
const real_number sign ( do_forward ? 1.0 : -1.0 );
const size_type length = z.size();

// allocate result
complex_vector result;
result.reserve( length );

// do matrix multiplication
for ( size_type i = 0; i < length; ++ i ) {
DFT_Pseudoiterator coefficient_iter ( sign, i, length );
result.push_back( std::inner_product( z.begin(), z.end(),
coefficient_iter,
complex_number(0.0) ) );
}

// rescale if forward
if ( do_forward ) {
for ( size_type i = 0; i < length; ++ i ) {
result[i] /= static_cast<real_number>( length );
}
}

// copy data back (swap is fast!)
z.swap( result );
}

void FFT ( complex_vector & z ) {
size_type const length = z.size();
if ( is_odd( length ) ) {
// too bad, were odd
DFT( true, z );
} else {
// good, we are even.
// Let's divide and conquer!
complex_vector odd;
complex_vector even;
odd.reserve( length / 2 );
even.reserve( length / 2 );
for ( complex_vector::const_iterator run = z.begin();
run != z.end(); ++ run ) {
even.push_back( *run );
++ run;
odd.push_back( *run );
}
FFT( even );
FFT( odd );
DFT_Pseudoiterator omega_iter ( 1.0, 1, length );
for ( size_type i = 0; i < length / 2; ++ i, ++ omega_iter ) {
z[i] = 0.5 * ( even[i] + *omega_iter * odd[i] );
// the next line works because omega^(length/2) = -1 !!!
z[i + length/2] = 0.5 * ( even[i] - *omega_iter * odd[i] );
}
}
}

void IFFT ( complex_vector & z ) {
size_type const length = z.size();
if ( is_odd( length ) ) {
// too bad, were odd
DFT( false, z );
} else {
// good, we are even.
// Let's divide and conquer!
complex_vector odd;
complex_vector even;
odd.reserve( length / 2 );
even.reserve( length / 2 );
for ( complex_vector::const_iterator run = z.begin();
run != z.end(); ++ run ) {
even.push_back( *run );
++ run;
odd.push_back( *run );
}
IFFT( even );
IFFT( odd );
DFT_Pseudoiterator omega_iter ( -1.0, 1, length );
for ( size_type i = 0; i < length / 2; ++ i, ++ omega_iter ) {
z[i] = even[i] + *omega_iter * odd[i];
// the next line works because omega^(length/2) = -1 !!!
z[i + length/2] = even[i] - *omega_iter * odd[i];
}
}
}

template < typename ComplexIterator >
void FFT ( bool do_forward,
ComplexIterator const & from, ComplexIterator const & to ) {
complex_vector z ( from, to );
if ( do_forward ) {
FFT( z );
} else {
IFFT( z );
}
std::copy( z.begin(), z.end(), from );
}

} // namespace fft

int main ( void ) {
fft::complex_vector z;
z.push_back( fft::complex_number( 0.5, 2 ) );
z.push_back( fft::complex_number( 0.7, 1.4 ) );
z.push_back( fft::complex_number( 0.9, -2.1 ) );
z.push_back( fft::complex_number( 0.9, -2.1 ) );
z.push_back( fft::complex_number( -1.2, 1.1 ) );
z.push_back( fft::complex_number( -1.2, 1.1 ) );

fft::DFT( true, z );
std::copy( z.begin(), z.end(),
std::ostream_iterator< fft::complex_number >( std::cout, "\n" ) );
fft::DFT( false, z );
std::copy( z.begin(), z.end(),
std::ostream_iterator< fft::complex_number >( std::cout, "\n" ) );
fft::FFT( true, z.begin(), z.end() );
std::copy( z.begin(), z.end(),
std::ostream_iterator< fft::complex_number >( std::cout, "\n" ) );
fft::FFT( false, z.begin(), z.end() );
std::copy( z.begin(), z.end(),
std::ostream_iterator< fft::complex_number >( std::cout, "\n" ) );
}

I did not bother putting in many comments. So, I cannot claim that I have
proof that this code is correct -- but I am reasonably confident, except
that I may have switched FFT and IFFT.

Also, please note that this code is in no way optimized for performance
(except for a slight change in DFT_Pseudoiterator that now avoids
unnecessary calls to trig-functions). I am pretty sure that one can
eliminate many of the copies that are made just to separate the even and
odd indices. I think, that one should be able to do most of this in place
(maybe some fancy slice tricks from the <valarray> header would be cool).
But any of that would probably make the code more difficult to understand.

Also, a final version should probably be completely templated on the
real_number type.
Best

Kai-Uwe Bux
Oct 12 '05 #12

This discussion thread is closed

Replies have been disabled for this discussion.