rtfft.c

FFT FUNCTION "C" SOURCE CODE

				{
					c2 = c1;
					if (nt == 1)
					{
						c1 = c1 * ck + s1 * sk;
						s1 = s1 * ck - c2 * sk;
					}
					else
					{
						c1 = (c1 - s1) * sq;
						s1 = (s1 + c2) * sq;
					}
					kb = kb + k4;
					if (kb <= kn)
						goto l1;
					else
						break;
				}
				if (nt == -1)
				{
					k = 2;
					continue;
				}
				if (cc[j] <= indx)
				{
					indx = indx - cc[j];
					j--;
					if (cc[j] <= indx)
					{
						indx = indx - cc[j];
						j--;
						k = k + 2;
					}
					else
					{
						indx = cc[j] + indx;
						j = mk;
					}
				}
				else
				{
					indx = cc[j] + indx;
					j = mk;
				}
				if (j < mk)
					continue;
				k = 0;
				nt = 3;
				kb = kb + k4;
				if (kb > kn)
					break;
			}
			while (1);
		}
		k = (m / 2) * 2;
		if (k != m)
		{
			k2 = kn;
			k0 = j = kn - cc[k];

			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 (k2 > j);
		}
	}
	while (kn < n);
	free (cc);
}

/**************************************************************************/
void rtreorder ( realtype *x, realtype *y, int n, int m, int ks, int reel)
{
	int	i, j, indx, k, kk, kb, k2, ku, lim, p;
	realtype  temp;

	int	*cc, *lst;

	cc =  (int far *)calloc (m + 1, sizeof (int));
	lst =  (int far *)calloc (m + 1, sizeof (int));
	
	cc [m] = ks;
	for (k = m; k > 0; k--)
		cc [k-1] = cc [k] / 2;
	p = j = m - 1;
	i = kb = 0;
	if (reel)
	{
		ku = n - 2;
		for (k = 0; k <= ku; k += 2)
		{
			temp = x [k + 1];
			x [k + 1] = y [k];
			y [k] = temp;
		}
	}
	else
		m--;
	lim = (m + 2) / 2;
	if (p > 0)
	do
	{
		ku = k2 = cc [j] + kb;
		indx = cc [m - j];
		kk = kb + indx;

		do
		{
			 k = kk + indx;
   	   		 do
			 {
		         temp = x [kk];
				  x [kk] = x [k2];
				  x [k2] = temp;
        		  temp = y [kk];
				  y [kk] = y [k2];
				  y [k2] = temp;
				  kk++;
				  k2++;
			 }
			 while (kk < k);

			 kk = kk + indx;
			 k2 = k2 + indx;
		}
		while (kk < ku);

		if (j > lim)
		{
			j--, i++;
			lst [i] = j;
			continue;
		}
		kb = k2;
		if (i > 0)
		{
			j = lst [i];
			i--;
			continue;
		}
		if (kb < n)
			j = p;
		else
			break;
    }
	while (1);

	free (cc);
	free (lst);
}

/**************************************************************************/

void rttrans (realtype *x, realtype *y, int n, int eval)
{
	int		  k, nk, nn;
	realtype  x1, ab, ba, y1, re, im, ck, sk, dc, ds, r;

	nn = n / 2;
	r = (realtype) 3.14159265359 / n;
	ds = (realtype) sin (r);
	r = (realtype) 2.0 * sin (0.5 * r);
	r = - r * r;
	dc = (realtype) -0.5 * r;
	ck = (realtype) 1.0;
	sk = (realtype) 0.0;

	if (!eval)
	{
		x[n] = x[0]; y [n] = y [0];
	}
	for (k = 0; k <= nn; k++)
	{
		nk = n - k;
		x1 = x[k] + x[nk];
		ab = x[k] - x[nk];
		ba = y[k] + y[nk];
		y1 = y[k] - y[nk];

		re = ck * ba + sk * ab;
		im = sk * ba - ck * ab;
		y [nk] = im - y1;
		y [k] = im + y1;
		x [nk] = x1 - re;
		x [k] = x1 + re;
		dc = r * ck + dc;
		ck = ck + dc;
		ds = r * sk + ds;
		sk = sk + ds;
	}
}

/**************************************************************************/

void rtfftcalc(realtype *xreal, realtype *yimag, int numdat)
{

   rtfft(xreal, yimag, numdat, 0);
}

void rtfftinvcalc(realtype *xreal, realtype *yimag, int numdat)
{

   rtfft(xreal, yimag, numdat, 1);
}

realtype rtsquareandsum (realtype a, realtype b)
{
   realtype fret;

   fret = a * a + b * b;
   return(fret);
}

void rtpowerspectrumcalc(realtype *xreal, realtype *yimag,
	int numdat, realtype delta)
{
	realtype timespan, normal, n;
	int i;

	n = numdat;
	normal =  (realtype) 1.0 / (n * n);
	rtfft (xreal, yimag, numdat, 0);
	xreal[0] = rtsquareandsum (xreal[0], yimag[0]) * normal;

	for (i = 1; i <= (numdat / 2) - 1; i++)
	{
		xreal[i] = normal * (rtsquareandsum(xreal[i], yimag[i]) +
		rtsquareandsum (xreal[numdat - i], yimag[numdat - i]));
	}
	i = numdat / 2;
	xreal[i] = rtsquareandsum(xreal[i], yimag[i]) * normal;
	timespan = (realtype) ((1.0*numdat) * delta);
	for (i = 0; i <= (numdat/2); i++)
	{
		yimag[i] = (realtype) (1.0*i) / timespan;
	}
}

/*--------------------------------------------------------------*/
/* Calculates the FFT magnitude associated with a given harmonic index. */
realtype rtfftmagnitude(realtype *xr, realtype *yi, int n, int i)
{
	realtype result, nr, rr, ii;

	nr = (realtype) n;
	if (i > n/2)
		result = (realtype) 0;
	else
	if (i == 0)
		result = sqrt(xr[0] * xr[0] + yi[0] * yi[0]) / nr;
	else
	{
		rr = fabs(xr[i]) + fabs(xr[n-i]);
		ii = fabs(yi[i]) + fabs(yi[n-i]);
		result =  sqrt(rr * rr + ii * ii) / nr;
	}
	return( result);
}

/*--------------------------------------------------------------*/
/* Calculates the FFT phase associated with a given harmonic index. */
realtype rtfftphase (realtype *xr, realtype *yi, int i)
{
	realtype result;

	result = atan2 (yi[i], xr[i]);
	return (result);
}

/*--------------------------------------------------------------*/
/* Calculates the FFT frequence associated with a given harmonic index
and sampling frequency. */
realtype rtfftfrequency (int n, realtype samplefreq, int i)
{
	realtype result;

	if (i <= n/2)
		result = samplefreq * (realtype)i/ (realtype)n;
	else
	{
		i = n - i;
		result = -samplefreq * (realtype)i/ (realtype)n;
	}
	return(result);
}



/* single pass of 2d fft done in row column order */
void rtfft2drccalc(matstruct *xreal,   /* input and output of real part of 2D array of data */
		 matstruct *yimag,   /* input and output of imaginary part of 2d data */
		 int c, int r,       /* number of columns, number of rows */
		 int flag)           /* 0 = forward fft, 1 = inverse fft */
{ int n,a,j;
  realtype *yvector, *xvector;
  if (r>=c) n=r; else n = c;
  xvector = rtgetmem(n); if (rtbadmem(xvector)) exit(0);
  yvector = rtgetmem(n); if (rtbadmem(yvector)) exit(0);
  for (j = 0; j <  r; j++)
  {
    for (a = 0; a < c; a++)
    {
      xvector[a] =rdmat(xreal,j,a);
      yvector[a] =rdmat(yimag,j,a);
    }
    if (flag == 0)
      rtfftcalc(xvector, yvector, c);
    else
      rtfftinvcalc(xvector, yvector, c);
    for (a = 0; a < c; a++)
    {
      wrmat(xreal,j,a,xvector[a]);
      wrmat(yimag,j,a,yvector[a]);
    }
  }
  free(xvector);
  free(yvector);
}

/* single pass of 2d fft done in column row order */
void rtfft2dcrcalc(matstruct *xreal,   /* input and output of real part of 2d array of data */
		   matstruct *yimag,   /* input and output of imaginary part of 2d data */
		 int c, int r,       /* number of columns, number of rows */
		 int flag)           /* 0 = forward fft, 1 = inverse fft */
{ int n,a,i;
  realtype *yvector, *xvector;

  if (r>=c) n=r; else n = c;
  xvector = rtgetmem(n); if (rtbadmem(xvector)) exit(0);
  yvector = rtgetmem(n); if (rtbadmem(yvector)) exit(0);
  for (i = 0; i < c; i++)
  {
    for (a = 0; a < r; a++)
    {
      xvector[a]=rdmat(xreal,a,i);
      yvector[a]=rdmat(yimag,a,i);
    }
    if (flag == 0)
      rtfftcalc(xvector, yvector, r);
    else
      rtfftinvcalc(xvector, yvector, r);
    for (a = 0; a < r; a++)
    {
      wrmat(xreal,a,i,xvector[a]);
      wrmat(yimag,a,i,yvector[a]);
    }
  }
  free(xvector);
  free(yvector);
}


void rtfft2dcalc (matstruct *xreal, matstruct *yimag,
			   int c, int r, int flag)
{
  if (flag)
  {
    rtfft2dcrcalc (xreal,yimag,c,r,flag);
    rtfft2drccalc (xreal,yimag,c,r,flag);
  }
  else
  {
    rtfft2drccalc (xreal,yimag,c,r,flag);
    rtfft2dcrcalc (xreal,yimag,c,r,flag);
  }
}

void rtfreqresponse (realtype *filtcoef, realtype *a, int n,int k )
{
	int i, j, m, n2;
	realtype att, am, q;

	q = (realtype) M_PI/k;
	am = (realtype) ((n+1)/2.0);
	m = (n+1) / 2;
	n2 = n / 2;
	for (j = 1; j <= k; j++)
	{
		att =  (realtype) 0.0;
		if (am == m)  att = filtcoef[m-1]/2.0;
		for (i = 1; i <= n2; i++)
			att = (realtype) (att + filtcoef[i-1] * cos(q*(am-i)*(j-1)));
		a[j-1] = (realtype) 2.0 * att;
	}
}

void rtfreqsamplefir (realtype *filtcoef, realtype fp, int n,
		int dc, int filttype, int win)
{
	int np, i, j, m, m1, n2;
	realtype lowbin, highbin, xt, am, q, *a;

	a = rtgetmem(n);  if (rtbadmem(a)) exit(0);

	m = (n+1) / 2;
	am = (realtype) ((n+1)/2.0);
	m1 = n / 2 + 1;
	q = (realtype) (2.0*M_PI/n);
	n2 = n / 2;
	if (dc == 1)
		np = (int) floor(n * fp + 0.5);
	else
		np = (int) floor(n * fp + 1.0);

	if (filttype == 1)
	{
		lowbin = (realtype) 0.0;
		highbin = (realtype) 1.0;
	}
	else
	{
		lowbin = (realtype) 1.0;
		highbin = (realtype) 0.0;
	}
	for (j = 0; j <= np; j++)
		a[j] = lowbin;

	for (j = np+1; j <= m1; j++)
		a[j] = highbin;

	if (dc == 0)
	{
		for (j = 1; j <= m; j++)
		{
			xt = a[1] / 2.0;
			for (i = 2; i <= m; i++)
				xt = xt + a[i]* cos(q * (am-j) * (i-1));
			filtcoef[j-1] = (realtype) (2.0 * xt / n);
			filtcoef[n-j] = filtcoef[j-1];
		}
	}
	else
	{
		for (j = 1; j <= m; j++)
		{
			xt = (realtype) 0.0;
			for (i = 1; i <= n2; i++)
				xt = xt + a[i] * cos(q * (am-j) * (i-0.5));
			if (am != m)
				xt = xt + a[m] * cos(M_PI * (am-j))/2.0;
			filtcoef[j-1] = (realtype) (2.0 * xt / n);
			filtcoef[n-j] = filtcoef[j-1];
		}
	}

	if (win > 0 ) rtwindowdata (filtcoef, n, win);
	free(a);
}


void rtconvolve (realtype *filtcoef, realtype *x, realtype *y,
	int filtlen, int alen)
{
	int n, m;
	realtype summ;

	for (n = 0; n <= (alen+filtlen-2); n++)
	{
		summ = (realtype) 0.0;
		for (m = 0; m <= filtlen-1; m++)
		{
			if ((n-m >= 0) && (n-m) <= (alen-1))
			{
				summ = summ + filtcoef[m] * x[n-m];
			}
		}
		y[n] = summ;
	}
}