rtfft.c

FFT FUNCTION "C" SOURCE CODE

/*****************************************************************/
/*							         */
/*   copyright (c) quinn-curtis, 1992                            */
/*							         */
/*   filename:  fft.c                                         */
/*   revision:  3.0                                              */
/*   date:      2/04/92                                          */
/*							         */
/*   description: -fourier analysis 	*/
/*							         */
/*****************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "rtstdhdr.h"

#ifndef RECTANG
#define RECTANG 0
#define PARZEN  1
#define HANNING 2
#define WELCH   3
#define HAMMING 4
#define EXACTB  5
#endif


void rtfft2 (realtype *aa, realtype *bb, int n, int m, int ks);
void rtrfft2 (realtype *aa, realtype *bb, int n, int m, int ks);
void rtreorder (realtype *aa, realtype *bb, int n, int m, int ks, int reel);
void rttrans (realtype *aa, realtype *bb, int n, int eval);

realtype *rtgetmem(int n)
{      realtype *xt;
       xt = (realtype *) calloc(n,sizeof(realtype));
       return(xt);
}

int rtbadmem(void *p)
{
       if (p==NULL){
	    printf("out of memory\n");
	    return(1);
       }
       else
	    return(0);

}

realtype rtparzen (realtype n, realtype j)
{
	realtype fret;

	fret = (realtype) 1.0 - fabs((j - 0.5 * (n - 1.0)) / (0.5 * (n - 1.0)));
	return(fret);
}


realtype rthanning (realtype n, realtype j)
{
	realtype fret;

	fret =  (realtype) 0.5 * (1.0 - cos(2.0 * M_PI * j / (n-1.0)));
	return(fret);
}

realtype rthamming (realtype n, realtype j)
{
	realtype fret;

	fret =  (realtype) 0.54 - 0.46 * cos(2.0 * M_PI * j / (n-1.0));
	return(fret);
}

realtype rtexactblackman (realtype n, realtype j)
{
	realtype fret;

	fret =  (realtype) 0.42 - 0.50 * cos(2.0 * M_PI * j /(n-1.0)) + 0.08 * cos(4.0 * M_PI * j / (n-1.0));
	return(fret);
}

realtype rtwelch (realtype n, realtype j)
{
	realtype fret;
	fret =  (realtype) (j - 0.5 * (n - 1.0)) / (0.5 * (n + 1.0));
	fret = (realtype) 1.0 - fret*fret;
	return(fret);
}

void rtwindowdata (realtype *x, int numdat, int window)
{
   int i;
   realtype multiplier;

   for ( i = 0; i < numdat; i++ )
   {
	  switch ( window )
	  {
		 case RECTANG:
			multiplier = (realtype) 1.0;
			break;
		 case PARZEN:
			multiplier = rtparzen (numdat, i);
			break;
		 case HANNING:
			multiplier = rthanning (numdat, i);
			break;
		 case WELCH:
			multiplier = rtwelch (numdat, i);
			break;
		 case HAMMING:
			multiplier = rthamming (numdat, i);
			break;
		 case EXACTB:
			multiplier = rtexactblackman (numdat, i);
			break;
		 default:
			multiplier =  (realtype) 1.0;
			break;
	  }
	  x[i] = multiplier * x[i];
   }
}

void rtwindowfftdata (realtype *xdata, realtype *ydata, int numdat, int window)
{
	rtwindowdata (xdata, numdat, window);
	rtwindowdata (ydata, numdat, window );
}



/**************************************************************************/
/* complex fourier transform
/* number of points n must be a power of 2
/**************************************************************************/
void rtfft (realtype *xr, realtype *yi, int n, int inverse)
{
	int nn, j;
	int		 m;
	realtype pp, qq;

	m = -1;
	nn = n;

	/* determine power of 2 */
	while (nn > 0)
	{
		nn = nn >> 1;
		m++;
	}

	rtfft2 (xr, yi, n, m, n);
	rtreorder (xr, yi, n, m, n, 0);
	if (inverse)
	{
		pp = (realtype) 1.0 / n;
		qq = -pp;
		for (j = 0; j < n; j++)
		{
			xr [j] *= pp;
			yi [j] *= qq;
		}
	}
	else
	{
		for (j = 0; j < n; j++)
			yi [j] = -yi[j];
	}
}

/**************************************************************************/
/* 	real array is passed in x, cos and sin coefficients are returned in
/*  x and sinc
/**************************************************************************/
void rtrealfft ( realtype *x, realtype *sinc, int n, int inverse)
{
	int 	 nn, j;
	int		 m;
	realtype pp;

	m = -1;
	nn = n / 2;

	/* determine power of 2 */
	while (nn > 0)
	{
		nn = nn >> 1;
		m++;
	}
	nn = n / 2;

	if (inverse)
	{
		pp = (realtype) 0.5 / nn;
		rttrans (x, sinc, nn, 1);
		rtfft2 (x, sinc, nn, m, nn);
		for (j = 0; j < nn; j++)
		{
			x [j] *= pp;
			sinc [j] *= -pp;
		}
		rtreorder (x, sinc, nn, m, nn, 1);
		for (j = 0; j < nn; j++)		/* combine */
			x [j + nn] = sinc [j];
	}
	else
	{
		for (j = 0; j < nn; j++)		/* split array in 2 */
			sinc [j] = x [j + nn];

		rtreorder (x, sinc, nn, m, nn, 1);
		rtrfft2 (x, sinc, nn, m, 1);
		for (j = 0; j < nn; j++)
		{
			x [j] *= 0.5;
			sinc [j] *= 0.5;
		}
		rttrans (x, sinc, nn, 0);
		for (j = 0; j < nn; j++)
			sinc [j] = -sinc[j];

		for (j = 1; j < nn; j++)
		{
			x [j + nn] = x [nn - j];
			sinc [j + nn] = -sinc [nn - j];
		}
	}
}
/**************************************************************************/
/* fft for one variable of dimension n = 2^m
/**************************************************************************/
void rtfft2 (realtype *x, realtype *y, int n, int m, int ks)
{
	int  	   indx, k0, k1, k2, k3, span, j, k, kb, kn, mm, mk;
	realtype	rad, c1, c2, c3, s1, s2, s3, ck, sk, sq;
	realtype	x0, x1, x2, x3, y0, y1, y2, y3;
	int	*cc;

	cc = (int far *) calloc (m + 1, sizeof (int));

	sq = (realtype) 0.707106781187;
	sk = (realtype) 0.382683432366;
	ck = (realtype) 0.92387953251;

	cc[m] = ks;
	mm = (m / 2) * 2;
	kn = 0;
	for (k = m - 1; k >= 0; k--)
		cc[k] = cc[k+1] / 2;
	rad = (realtype) 6.28318530718 / (cc[0] * ks);
	mk = m - 5;

	do
	{
		kb = kn;
		kn = kn + ks;

		if (mm != m)
		{
			k2 = kn;
			k0 = cc[mm] + kb;
			do
			{
				k2--;
				k0--;
				x0 = x [k2];
				y0 = y [k2];
				x[k2] = x[k0] - x0;
				x[k0] = x[k0] + x0;
				y[k2] = y[k0] - y0;
				y[k0] = y[k0] + y0;
			}  while (k0 > kb);
		}

		c1 = (realtype) 1.0, s1 = (realtype) 0.0;
		indx = 0;
		k = mm - 2;
		j = 3;

		if (k >= 0)
			goto l2;
		else
			if (kn < n)
				continue;
			else
				break;

l1:	   	while (1)
		{
			if (cc[j] <= indx)
			{
				indx = indx - cc[j];
				j--;
				if (cc[j] <= indx)
				{
					indx = indx - cc[j];
					j--;
					k += 2;
				}
			}
			else
				break;
		}

		indx = cc[j] + indx;
		j = 3;

l2:		span = cc[k];
		if (indx != 0)
		{
			c2 = indx * span * rad;
			c1 = (realtype) cos (c2);
			s1 = (realtype) sin (c2);
l5:
			c2 = c1 * c1 - s1 * s1;
			s2 = (realtype) 2.0 * s1 * c1;
			c3 = c2 * c1 - s2 * s1;
			s3 = c2 * s1 + s2* c1;
		}
		for (k0 = kb + span - 1; k0 >= kb; k0--)
		{
			k1 = k0 + span;
			k2 = k1 + span;
			k3 = k2 + span;
			x0 = x[k0];
			y0 = y [k0];
			if (s1 == 0)
			{
				x1 = x[k1]; y1 = y [k1];
				x2 = x[k2]; y2 = y [k2];
				x3 = x[k3]; y3 = y [k3];
			}
			else
			{
				x1 = x[k1] * c1 - y [k1] * s1;
				y1 = x[k1] * s1 + y [k1] * c1;
				x2 = x[k2] * c2 - y [k2] * s2;
				y2 = x[k2] * s2 + y [k2] * c2;
				x3 = x[k3] * c3 - y [k3] * s3;
				y3 = x[k3] * s3 + y [k3] * c3;
			}
			x [k0] = x0 + x2 + x1 + x3; y[k0] = y0 + y2 + y1 + y3;
			x [k1] = x0 + x2 - x1 - x3; y[k1] = y0 + y2 - y1 - y3;
			x [k2] = x0 - x2 - y1 + y3; y[k2] = y0 - y2 + x1 - x3;
			x [k3] = x0 - x2 + y1 - y3; y[k3] = y0 - y2 - x1 + x3;
		}
		if (k > 0)
		{
			k = k - 2;
			goto l2;
		}
		kb = k3 + span;
		if (kb < kn)
		{
			if (j == 0)
			{
				k = 2;
				j = mk;
				goto l1;
			}
			j--;
			c2 = c1;
			if (j == 1)
			{
				c1 = c1 * ck + s1 * sk;
				s1 = s1 * ck - c2 * sk;
			}
			else
			{
				c1 = (c1 - s1) * sq;
				s1 = (s1 + c2) * sq;
			}
			goto	l5;
		}
	}
	while (kn < n);

	free (cc);
}
/**************************************************************************/

void rtrfft2 ( realtype *x, realtype *y, int n, int m, int ks)
{
	int  		k0, k1, k2, k3, k4, span, j, indx, k, kb, nt, kn, mk;
	realtype	rad, c1, c2, c3, s1, s2, s3, ck, sk, sq;
	realtype	x0, x1, x2, x3, y0, y1, y2, y3, re, im;
	int	*cc;

	cc =  (int far *)calloc (m + 1, sizeof (int));

	sq = (realtype) 0.707106781187;
	sk = (realtype) 0.382683432366;
	ck = (realtype) 0.92387953251;

	cc[0] = ks;
	kn = 0;
	k4 = 4 * ks;
	mk = m - 4;

	for (k = 1; k <= m; k++)
	{
		ks = ks + ks;
		cc[k] = ks;
	}
	rad = (realtype) 3.14159265359 / (cc[0] * ks);

	do
	{
		kb = kn + 4;
		kn = kn + ks;
		if (m != 1)
		{
			k = indx = 0;
			j = mk;
			nt = 3;
			c1 = (realtype) 1.0, s1 = (realtype) 0.0;

			do
			{
			   span = cc[k];
				if (indx != 0)
				{
					c2 = indx * span * rad;
					c1 = (realtype) cos (c2);
					s1 = (realtype) sin (c2);
l1:
					c2 = c1 * c1 - s1 * s1;
					s2 = (realtype) 2.0 * s1 * c1;
					c3 = c2 * c1 - s2 * s1;
					s3 = c2 * s1 + s2 * c1;
				}
				else
					s1 = 0;
				k3 = kb - span;
				do
				{
		           k2 = k3 - span;
					k1 = k2 - span;
					k0 = k1 - span;
		
					x0 = x[k0]; y0 = y [k0];
					x1 = x[k1]; y1 = y [k1];
					x2 = x[k2]; y2 = y [k2];
					x3 = x[k3]; y3 = y [k3];

					x [k0] = x0 + x2 + x1 + x3; y[k0] = y0 + y2 + y1 + y3;
					if (s1 == 0)
					{
						x [k1] = x0 - x1 - y2 + y3;
						y[k1] = y0 - y1 + x2 - x3;
						x [k2] = x0 + x1 - x2 - x3;
			y[k2] = y0 + y1 - y2 - y3;
						x [k3] = x0 - x1 + y2 - y3;
                        y[k3] = y0 - y1 - x2 + x3;
					}
					else
					{
						re = x0 - x1 - y2 + y3;
						im = y0 - y1 + x2 - x3;
						x [k1] = re * c1 - im * s1;
                        y[k1] = re * s1 + im * c1;
						re = x0 + x1 - x2 - x3;
						im = y0 + y1 - y2 - y3;
						x [k2] = re * c2 - im * s2;
						y[k2] = re * s2 + im * c2;
						re = x0 - x1 + y2 - y3;
						im = y0 - y1 - x2 + x3;
						x [k3] = re * c3 - im * s3;
                        y[k3] = re * s3 + im * c3;
					}
					k3++;
		      		}
				while (k3 < kb);

				nt--;
				if (nt >= 0)