An efficient storage scheme for fourier data
The fourier transform (FT), of course, maps a real signal to the complex plane. Unfortunately, this also means that the result is twice as large as the original signal. As it turns out, though, fully half of that information is completely redundant! Let's see how this is possible. Suppose we are working with a signal containing 8 samples. After performing the FT, we see a pattern emerge:
Reals:
[W][A][b][C][X][C][b][A]
Imaginaries:
[Y][D][E][F][Z][-F][-E][-D]
Notice how all of the coefficients straddling X are "reflections" of eachother, and those of Z are "negative reflections". Also note that Y and Z are *always* zero (or extremely close to it). In other words, all we really need to reconstruct the two planes are:
Composite:
[W][A][b][C][X][-F][-E][-D]
Translated to C:
The only constraint, of course, it that the sample size must be a power of two.Code:int is_power_of_two( int n ) { int bits = 0; while( n ) { if( n & 1 ) ++bits; n >>= 1; } return bits == 1; } void compress_fft_result( const double reals[ ], const double imags[ ], double composite[ ], int size ) { if( !is_power_of_two( size ) ) return; for( int idx = 0, mid = size >> 1; idx < size; ++idx ) { if( idx <= mid ) composite[ idx ] = reals[ idx ]; else composite[ idx ] = imags[ idx ]; } } void decompress_fft_result( const double composite[ ], double reals[ ], double imags[ ], int size ) { if( !is_power_of_two( size ) ) return; int mid = size >> 1; for( int idx = 0; idx < size; ++idx ) { if( idx <= mid ) { reals[ idx ] = composite[ idx ]; imags[ idx ] = -composite[ mid + ( mid - idx ) ]; } else { reals[ idx ] = composite[ mid - ( idx - mid ) ]; imags[ idx ] = composite[ idx ]; } } imags[ 0 ] = imags[ mid ] = 0; }
Cheers!
Originally Posted by Salem
