Plan 9 from Bell Labs’s /usr/web/sources/contrib/cnielsen/bladeenc/subs.c

Copyright © 2021 Plan 9 Foundation.
Distributed under the MIT License.
Download the Plan 9 distribution.


/*
			(c) Copyright 1998-2001 - Tord Jansson
			======================================

		This file is part of the BladeEnc MP3 Encoder, based on
		ISO's reference code for MPEG Layer 3 compression, and might
		contain smaller or larger sections that are directly taken
		from ISO's reference code.

		All changes to the ISO reference code herein are either
		copyrighted by Tord Jansson ([email protected])
		or sublicensed to Tord Jansson by a third party.

	BladeEnc is free software; you can redistribute this file
	and/or modify it under the terms of the GNU Lesser General Public
	License as published by the Free Software Foundation; either
	version 2.1 of the License, or (at your option) any later version.



	------------    Changes    ------------

	2000-11-04  Andre Piotrowski

	-	speed up: don� calculate the unneeded phi values!

	2000-11-22	ap

	-	speed up: slightly improved way to handle the 'special case'

	2000-12-05  ap

	-	speed up: implemented prepacking

	2000-12-11  ap

	-	speed up: avoid data reordering

	2000-12-12  ap

	-	moved fft configuration switches to "l3psy.h"

	2001-01-12  ap

	-	use some explicit type casting to avoid compiler warnings

	2001-04-06  ap

	-	implemented radix-4 FFT
*/

/*****************************************************************************
 * FFT computes fast fourier transform of BLKSIZE samples of data            *
 *   based on the decimation-in-frequency algorithm described in "Digital    *
 *   Signal Processing" by Oppenheim and Schafer, refer to pages 304         *
 *   (flow graph) and 330-332 (Fortran program in problem 5).                *
 *                                                                           *
 *   required constants:                                                     *
 *         PI          3.14159265358979                                      *
 *         BLKSIZE     should be 2^^M for a positive integer M               *
 *                                                                           *
 *****************************************************************************/

#include	"common.h"
#include	"encoder.h"
#include	"l3psy.h"





/* The switches have been moved to "l3psy.h" */





int						fInit_fft;





#if NORMAL_FFT   /* handling sign exceptions */
#define P +
#define M -
#else
#define P -
#define M +
#endif


void fft (FLOAT x_real[], FLOAT x_imag[], FLOAT energy[], FLOAT phi[], int N)
{
 	int						i, j, k, off;
#if USED_VALUES_ONLY
	int						min_i, max_i;
#endif
	int						m, mr;
	double					xr0, xr1;
	double					xi0, xi1;
	double					ur0, ur1, ur2, ur3;
	double					ui0, ui1, ui2, ui3;
	double					tr0, tr1, tr2, tr3;
	double					ti0, ti1, ti2, ti3;
	double					er1, er2, er3;
	double					ei1, ei2, ei3;
	double					t_real, t_imag;
	double					sqrt05;
	static double			w_real[BLKSIZE*7/8], w_imag[BLKSIZE*7/8];
	int						N_ORG;

#if !REORDER_DATA
	static	int				swap_l[BLKSIZE/2+1];
	static	int				swap_s[BLKSIZE_s/2+1];
	int						*pSwap, a, b;
#endif

	double	t1, t2, t3, t4, t5, t6;

	if (fInit_fft == 0)
	{
		for (i = 0;  i < BLKSIZE*7/8;  i++)
		{
			w_real[i] =   cos(2*PI*i/BLKSIZE);
/* M */		w_imag[i] = M sin(2*PI*i/BLKSIZE);
		}

		sqrt05 = sqrt (0.5);

#if !REORDER_DATA
		j = 0;
		for (i = 0;  i < BLKSIZE/2-1;  i++)
		{
			swap_l[i] = j;  k = BLKSIZE/4;  while (k <= j) {j -= k;  k >>= 1;}  j += k;
		}
		swap_l[i] = i;  swap_l[i+1] = i+1;

		j = 0;
		for (i = 0;  i < BLKSIZE_s/2-1;  i++)
		{
			swap_s[i] = j;  k = BLKSIZE_s/4;  while (k <= j) {j -= k;  k >>= 1;}  j += k;
		}
		swap_s[i] = i;  swap_s[i+1] = i+1;
#endif

		fInit_fft++;
	}


#if REAL_SEQUENCE
	N_ORG = N;
	N >>= 1;
#if !PREPACKED
	/* packing the sequence to the half length */
	for (i = 0;  i < N;  i++)
	{
		x_real[i] = x_real[2*i];
		x_imag[i] = x_real[2*i+1];
	}
#endif
#endif



	off = BLKSIZE/N;

	for (m = N;  m >= 4;  m = mr)
	{
		mr = m >> 2;

			for (i = 0;  i < N;  i += m)
			{
				ur0 = x_real[i     ];  ui0 = x_imag[i     ];
				ur1 = x_real[i+mr  ];  ui1 = x_imag[i+mr  ];
				ur2 = x_real[i+mr*2];  ui2 = x_imag[i+mr*2];
				ur3 = x_real[i+mr*3];  ui3 = x_imag[i+mr*3];

				xr0 = ur0 + ur2;  xi0 = ui0 + ui2;
				xr1 = ur1 + ur3;  xi1 = ui1 + ui3;
				tr0 = xr0 + xr1;  ti0 = xi0 + xi1;
				tr2 = xr0 - xr1;  ti2 = xi0 - xi1;

				xr0 = ur0 - ur2;  xi0 = ui0 - ui2;
/* P */			xr1 = P (ui1 - ui3);  xi1 = P (ur3 - ur1);
				tr1 = xr0 + xr1;  ti1 = xi0 + xi1;
				tr3 = xr0 - xr1;  ti3 = xi0 - xi1;

				x_real[i     ] = tr0;  x_imag[i     ] = ti0;
				x_real[i+mr*2] = tr1;  x_imag[i+mr*2] = ti1;
				x_real[i+mr*1] = tr2;  x_imag[i+mr*1] = ti2;
				x_real[i+mr*3] = tr3;  x_imag[i+mr*3] = ti3;
			}

		k = off;
		for (j = 1;  j < mr;  j++)
		{
			er1 = w_real[  k];  ei1 = w_imag[  k];
			er2 = w_real[2*k];  ei2 = w_imag[2*k];
			er3 = w_real[3*k];  ei3 = w_imag[3*k];

			for (i = j;  i < N;  i += m)
			{
				ur0 = x_real[i     ];  ui0 = x_imag[i     ];
				ur1 = x_real[i+mr  ];  ui1 = x_imag[i+mr  ];
				ur2 = x_real[i+mr*2];  ui2 = x_imag[i+mr*2];
				ur3 = x_real[i+mr*3];  ui3 = x_imag[i+mr*3];

				xr0 = ur0 + ur2;  xi0 = ui0 + ui2;
				xr1 = ur1 + ur3;  xi1 = ui1 + ui3;
				tr0 = xr0 + xr1;  ti0 = xi0 + xi1;
				tr2 = xr0 - xr1;  ti2 = xi0 - xi1;

				xr0 = ur0 - ur2;  xi0 = ui0 - ui2;
/* P */			xr1 = P (ui1 - ui3);  xi1 = P (ur3 - ur1);
				tr1 = xr0 + xr1;  ti1 = xi0 + xi1;
				tr3 = xr0 - xr1;  ti3 = xi0 - xi1;

				x_real[i     ] = tr0;                    x_imag[i     ] = ti0;
				x_real[i+mr*2] = tr1 * er1 - ti1 * ei1;  x_imag[i+mr*2] = tr1 * ei1 + ti1 * er1;
				x_real[i+mr*1] = tr2 * er2 - ti2 * ei2;  x_imag[i+mr*1] = tr2 * ei2 + ti2 * er2;
				x_real[i+mr*3] = tr3 * er3 - ti3 * ei3;  x_imag[i+mr*3] = tr3 * ei3 + ti3 * er3;
			}
			k += off;
		}
		off <<= 2;
	}

	if (m == 2)
	{
		for (i = 0;  i < N;  i += 2)
		{
			xr0 = x_real[i+1];  x_real[i+1] = x_real[i] - xr0;  x_real[i] += xr0;
			xi0 = x_imag[i+1];  x_imag[i+1] = x_imag[i] - xi0;  x_imag[i] += xi0;
		}
	}


#if REORDER_DATA
	/* this section reorders the data to the correct ordering */
	j = 0;
	for (i = 0;  i < N-1;  i++)
	{
		if (i < j)
		{
			t_real = x_real[j];
			t_imag = x_imag[j];
			x_real[j] = x_real[i];
			x_imag[j] = x_imag[i];
			x_real[i] = (FLOAT) t_real;
			x_imag[i] = (FLOAT) t_imag;
		}
		k=N/2;
		while (k <= j)
		{
			j -= k;
			k >>= 1;
		}
		j += k;
	}
#else
	/*
		We don� reorder the data to the correct ordering,
		but access the data by the bit reverse order index array.
	*/
	pSwap = (N_ORG == BLKSIZE) ? swap_l : swap_s;
#endif


#if REAL_SEQUENCE
	/* unpacking the sequence */
	t_real = x_real[0];
	t_imag = x_imag[0];
	x_real[0] = (FLOAT) (t_real+t_imag);
	x_imag[0] = 0.0;
	x_real[N] = (FLOAT) (t_real-t_imag);
	x_imag[N] = 0.0;

	k = off = BLKSIZE/N_ORG;
	for (i = 1;  i < N/2;  i++)
	{
#if REORDER_DATA
#define		a	i
#define		b	(N-i)
#else
		a = pSwap[i];
		b = pSwap[N-i];
#endif
		t1 = x_real[a] + x_real[b];
		t2 = x_real[a] - x_real[b];
		t3 = x_imag[a] + x_imag[b];
		t4 = x_imag[a] - x_imag[b];
		t5 = t2*w_imag[k] + t3*w_real[k];
		t6 = t3*w_imag[k] - t2*w_real[k];

		x_real[a] = (FLOAT) (t1+t5) / 2.0;
		x_imag[a] = (FLOAT) (t6+t4) / 2.0;
		x_real[b] = (FLOAT) (t1-t5) / 2.0;
		x_imag[b] = (FLOAT) (t6-t4) / 2.0;

		k += off;
	}
	/* x_real[N/2] doesn� change */
	/* x_imag[N/2] changes the sign in case of a normal fft */
#if (NORMAL_FFT)
#if REORDER_DATA
	x_imag[i] = -x_imag[i];
#else
	x_imag[pSwap[i]] *= -1.0;
#endif   /* REORDER_DATA */
#endif   /* NORMAL_FFT */
	N = N_ORG;
#endif   /* REAL_SEQUENCE */


	/* calculating the energy and phase, phi */
#if USED_VALUES_ONLY
	if (N == BLKSIZE)
	{
		min_i = LONG_FFT_MIN_IDX;
		max_i = LONG_FFT_MAX_IDX;
	}
	else
	{
		min_i = SHORT_FFT_MIN_IDX;
		max_i = SHORT_FFT_MAX_IDX;
	}
#endif
#if REAL_SEQUENCE
	for (i = 0;  i <= N/2;  i++)
#else
	for (i = 0;  i < N;  i++)
#endif
	{
#if REORDER_DATA
#define		a	i
#else
		a = pSwap[i];
#endif
		energy[i] = x_real[a]*x_real[a] + x_imag[a]*x_imag[a];
		if(energy[i] <= 0.0005)
		{
			energy[i] = (FLOAT) 0.0005;        /* keep the identity */
			x_real[a] = (FLOAT) sqrt(0.0005);  /* energy[i] * cos(phi[i]) == x_real[i] */
			x_imag[a] = 0.0;                   /* energy[i] * sin(phi[i]) == x_imag[i] */
		}
#if USED_VALUES_ONLY
		if (i >= min_i  &&  i <= max_i)
#endif
		phi[i] = (FLOAT) atan2((double) x_imag[a], (double) x_real[a]);
	}
}


#undef P
#undef M




Bell Labs OSI certified Powered by Plan 9

(Return to Plan 9 Home Page)

Copyright © 2021 Plan 9 Foundation. All Rights Reserved.
Comments to [email protected].