Отправлено Илья Гаврилов 29 мая 2003 г. 19:44 В ответ на: по поводу FFT(+) отправлено pahmed 29 мая 2003 г. 13:11

Вот пример реализации, может не самой эффективной, зато читаемой (если форматирование не съедет).

#include 
#include "complex.h"
double const pi = 3.14159265358979323846;
double const pi2 = 2 * pi;
void MakeTestData( double data[], int size )
{
	double w1 = 0.3 * pi2;
	double a1 = 0.2;
	double w2 = 0.33 * pi2;
	double a2 = 0.5;
	for ( int j = 0; j < size; j++ )
		data[j] = a1 * cos( w1 * j ) + a2 * cos( w2 * j );
}
void MakeWindow( double w[], int size )			// Welch, not normalized
{
	int sz2 = size / 2;
	for ( int j = 0; j < sz2; j++ )
	{
		double t = sz2 - 1 - j;
		t /= sz2;
		int k = size - 1 - j;
		w[j] = w[k] = 1 - t * t;
	}
}
void MakeTwiddle( complex tw[], int fft_size )
{
	int size = fft_size / 2;
	double theta = pi2 / fft_size;
	for ( int j = 0; j < size; j++ )
		tw[j] = complex( cos( j*theta ), sin( j*theta ) );
}
int BitRev( int j, int N )		// j - index to bit-reverse
{								// N - size of array: if N == 1024
								// the last 10 bits of j are reversed
	int rev = 0;
	
	for ( int k = 1; k < N; k <<= 1, j >>= 1 )
		rev = ( rev << 1 ) | ( j & 1 );
	return rev;
}
void Scramble( complex in[], complex out[], int N )
{
// place elements of in[] to out[] in bit-reversed order
// can't be done in place !!!
// N should be a power of 2 
	for ( int j = 0; j < N; j++ )
		out[BitRev(j,N)] = in[j];
}
void InplaceScramble( complex data[], int N )
{
	int j, k;
	for ( j = 0; j < N; j++ )
		if ( ( k = BitRev( j, N ) ) > j )
		{
			complex t = data[j];
			data[j] = data[k];
			data[k] = t;
		}
}
void WinScramble( complex in[], complex out[], double w[], int N )
{
	for ( int j = 0; j < N; j++ )
		out[BitRev(j,N)] = w[j] * in[j];
}
void RealToComplexWinScramble( double in[], complex out[], double w[], int N )
{
	for ( int j = 0; j < N; j++ )
		out[BitRev(j,N)] = complex( w[j] * in[j], 0 );
}
void RealWinScramble( double in[], complex out[], double w[], int N )
{
	for ( int j = 0; j < N; j += 2 )
		out[BitRev(j/2,N/2)] = complex( w[j] * in[j], w[j+1] * in[j+1] );
}

void ComplexFFT( complex data[], complex twiddle[], int N, int do_bitrev = 0 )
{
// data[] - array of N complex values in scrambled or normal order
// twiddle[] - array of N/2 complex twiddle factors
// N should be a power of 2
// do_bitrev - 0 if data[] already bit-reversed, 1 if in normal order
	if ( do_bitrev )
		InplaceScramble( data, N );
// pass 1 & 2
	int j, k;
	for ( j = 0; j < N; j += 4 )
	{
		complex a = data[j] + data[j+1];		// pass 1
		complex b = data[j] - data[j+1];
		complex c = data[j+2] + data[j+3];
		complex d = data[j+2] - data[j+3];
		data[j] = a + c;						// pass 2
		data[j+1] =	complex( b.re - d.im, b.im + d.re );
		data[j+2] = a - c;
		data[j+3] = complex( b.re + d.im, b.im - d.re );
	}
// passes 3 to log2(N)
	int groups;		// number of partial Fourier transforms
	int bflys;		// butterflies in each group
// each butterfly has two legs so ( groups * bflys * 2 ) == N
	for ( bflys = 4, groups = N >> 3; groups != 0; bflys <<= 1, groups >>= 1 )
	{
		int bn;		// butterfly number in a group
		int f;		// phase angle
		for ( bn = f = 0; bn < bflys; bn++, f += groups )
		{
			complex w = twiddle[f];			// same-numbered bflys in all groups
											// use the same twiddle factor
			for ( j = bn; j < N; j += bflys*2 )	// first leg of butterfly
			{									// do bfly bn in all groups
				k = j + bflys;				// second leg of butterfly
				complex t = data[k]*w;
				data[k] = data[j] - t;
				data[j] += t;
			}
		}
	}
}