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:
trdiff_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
/= 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 = 0.5 * ( even + *omega_iter * odd );
// the next line works because omega^(length/2) = -1 !!!
z[i + length/2] = 0.5 * ( even - *omega_iter * odd );
}
}
}
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 = even + *omega_iter * odd;
// the next line works because omega^(length/2) = -1 !!!
z[i + length/2] = even - *omega_iter * odd;
}
}
}
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:FT( true, z );
std::copy( z.begin(), z.end(),
std:stream_iterator< fft::complex_number >( std::cout, "\n" ) );
fft:FT( false, z );
std::copy( z.begin(), z.end(),
std:stream_iterator< fft::complex_number >( std::cout, "\n" ) );
fft::FFT( true, z.begin(), z.end() );
std::copy( z.begin(), z.end(),
std:stream_iterator< fft::complex_number >( std::cout, "\n" ) );
fft::FFT( false, z.begin(), z.end() );
std::copy( z.begin(), z.end(),
std:stream_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