ftp.nice.ch/pub/next/audio/apps/Tuner.NI.b.tar.gz#/Tuner/libfft.c

This is libfft.c in view mode; [Download] [Up]

/* libfft.c - fast Fourier transform library
**
** Copyright (C) 1989 by Jef Poskanzer.
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
*/

#include <stdio.h>
#include <math.h>

#include "libfft.h"

#define TWOPI ( M_PI + M_PI )

static int bitreverse[MAXFFTSIZE], bits;

/* initfft - initialize for fast Fourier transform
**
** b = power of two such that 2**nu = number of samples
*/
void initfft ( int b ) {
	register int i, j, k;

	if ( ( bits = b ) > LOG2_MAXFFTSIZE ) {
		fprintf ( stderr, "%d is too many bits, max is %d\n", bits, LOG2_MAXFFTSIZE );
	 	exit( 1 );
		}

	for ( i = ( 1 << bits ) - 1; i >= 0; --i ) {
		for ( j = 0, k = 0; j < bits; ++j ) {
			k *= 2;
			if ( i & ( 1 << j ) ) k++;
			}
		bitreverse[i] = k;
		}
	}

/* fft - a fast Fourier transform routine
**
** xr   real part of data to be transformed
** xi   imaginary part (normally zero, unless inverse transform in effect)
** inv  flag for inverse
*/

void fft ( float xr[], float xi[], int inv ) {
	int n, n2, i, k, kn2, l, p;
	float ang, s, c, tr, ti;

	n2 = ( n = ( 1 << bits ) ) / 2;

	for ( l = 0; l < bits; ++l ) {
		for ( k = 0; k < n; k += n2 ) {
			for ( i = 0; i < n2; ++i, ++k ) {
				p = bitreverse[k / n2];
				ang = TWOPI * p / n;
				c = cos( ang );
				s = sin( ang );
				kn2 = k + n2;
				if ( inv ) s = -s;
				tr = xr[kn2] * c + xi[kn2] * s;
				ti = xi[kn2] * c - xr[kn2] * s;
				xr[kn2] = xr[k] - tr;
				xi[kn2] = xi[k] - ti;
				xr[k] += tr;
				xi[k] += ti;
				}
			}
		n2 /= 2;
		}

	for ( k = 0; k < n; k++ ) {
		if ( ( i = bitreverse[k] ) <= k ) continue;
		tr = xr[k];
		ti = xi[k];
		xr[k] = xr[i];
		xi[k] = xi[i];
		xr[i] = tr;
		xi[i] = ti;
		}

	/* Finally, multiply each value by 1/n, if this is the forward transform. */
	if ( ! inv ) {
		register float f;

		for ( i = 0, f = 1.0 / n; i < n ; i++ ) {
			xr[i] *= f;
			xi[i] *= f;
			}
		}
	}

These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.