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 N1 on the powers
 omega^0, omega, omega^2, ..., omega^(N1) 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 Fouriertransfor ms 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<rea l_number>(i)
/ static_cast<rea l_number>(m) );
}
bool is_odd ( size_type m ) {
return( m % 2 != 0 );
}
class DFT_Pseudoitera tor
: public std::iterator< std::input_iter ator_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_Pseudoitera tor ( 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_Pseudoitera tor & operator++ ( void ) {
omega_to_the_ik *= omega_to_the_i;
++ k;
return( *this );
}
DFT_Pseudoitera tor operator++ ( int ) {
DFT_Pseudoitera tor dummy ( *this );
++ k;
return( dummy );
}
bool operator== ( DFT_Pseudoitera tor const & other ) const {
return( ( this>omega_to_the _i == other.omega_to_ the_i )
&&
( this>k == other.k ) );
}
bool operator!= ( DFT_Pseudoitera tor const & other ) const {
return( ! ( *this == other ) );
}
}; // DFT_Pseudoitera tor
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_Pseudoitera tor coefficient_ite r ( sign, i, length );
result.push_bac k( std::inner_prod uct( z.begin(), z.end(),
coefficient_ite r,
complex_number( 0.0) ) );
}
// rescale if forward
if ( do_forward ) {
for ( size_type i = 0; i < length; ++ i ) {
result[i] /= static_cast<rea l_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_Pseudoitera tor 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_Pseudoitera tor 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_ve ctor z;
z.push_back( fft::complex_nu mber( 0.5, 2 ) );
z.push_back( fft::complex_nu mber( 0.7, 1.4 ) );
z.push_back( fft::complex_nu mber( 0.9, 2.1 ) );
z.push_back( fft::complex_nu mber( 0.9, 2.1 ) );
z.push_back( fft::complex_nu mber( 1.2, 1.1 ) );
z.push_back( fft::complex_nu mber( 1.2, 1.1 ) );
fft::DFT( true, z );
std::copy( z.begin(), z.end(),
std::ostream_it erator< fft::complex_nu mber >( std::cout, "\n" ) );
fft::DFT( false, z );
std::copy( z.begin(), z.end(),
std::ostream_it erator< fft::complex_nu mber >( std::cout, "\n" ) );
fft::FFT( true, z.begin(), z.end() );
std::copy( z.begin(), z.end(),
std::ostream_it erator< fft::complex_nu mber >( std::cout, "\n" ) );
fft::FFT( false, z.begin(), z.end() );
std::copy( z.begin(), z.end(),
std::ostream_it erator< fft::complex_nu mber >( 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_Pseudoitera tor that now avoids
unnecessary calls to trigfunctions). 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
KaiUwe Bux