/*	Colmft.c

	Rocca's "row" Fourier transform, from Claerbout (1985) p. 73-74,
	modified for the indexing scheme used by C.

	Usage:		int n1, n2;
			float sign2, scale;
			struct complex *cx;

			colmft(n1, n2, cx, sign2, scale);

	Compile:	cc prog.c colmft.c -fswitch -lm -o prog.c

	Arguments:	n1	number of rows of 2-d matrix
				MUST BE A POWER OF 2
			n2	number of columns of matrix, fastest varying
			cx	matrix of complex elements cx[n1][n2]
			sign2	sign of the transform
			scale	scale factor to apply to data before transform
									*/
/* The usual includes */
#include <stdio.h>
#include <math.h>

/* Define complex multiplication and conjugate for declared struct complex */
#define  rmul(x,y)      (x.re * y.re - x.im * y.im)
#define  imul(x,y)      (x.im * y.re + x.re * y.im)
#define rcmul(x,y)      (x.re * y.re + x.im * y.im)
#define icmul(x,y)      (x.im * y.re - x.re * y.im)

/* Declare the structure to hold complex numbers */
struct complex 
	{
	float re; 
	float im;
	};

colmft(n1, n2, cx, sign2, scale)
int n1, n2;
float sign2, scale;
struct complex *cx;
	{
	int i, i1, i2, j, m, lstep, istep;
	float arg;
	struct complex cw, cdel, ctemp;

	/* Multiply each element in matrix by scale factor */
	for (i1=0; i1<n1; i1++)
		for (i2=0; i2<n2; i2++)
			{
			cx[i1*n2+i2].re *= scale;
			cx[i1*n2+i2].im *= scale;
			}

	j = 0;
	for (i=0; i<n1; i++)
		{
		if (i <= j)
			twid1(n2, &cx[i*n2+0], &cx[j*n2+0]);
		m = n1/2;
		while (j > m-1)
			{
			j = j - m;
			m = m/2;
			if (m < 1)
				break;
			}
		j = j + m;
		}
	lstep = 1;
	do
		{
		istep = 2*lstep;
		cw.re = 1;
		cw.im = 0;
		arg = sign2*3.14159265/lstep;
		cdel.re = (float)(cos((double)(arg)));
		cdel.im = (float)(sin((double)(arg)));
		for (m=0; m<lstep; m++)
			{
			for (i=m; i<n1; i+=istep)
				twid2(n2, &cw, &cx[i*n2+0],&cx[(i+lstep)*n2+0]);
			ctemp.re = rmul(cw,cdel);
			ctemp.im = imul(cw,cdel);
			cw.re = ctemp.re;
			cw.im = ctemp.im;
			}
		lstep = istep;
		}
	while (lstep < n1);

	return(0);
	}

twid1(n, cx, cy)
int n;
struct complex *cx, *cy;
	{
	int i;
	struct complex ct;

	for (i=0; i<n; i++)
		{
		ct.re = cx[i].re;
		ct.im = cx[i].im;
		cx[i].re = cy[i].re;
		cx[i].im = cy[i].im;
		cy[i].re = ct.re;
		cy[i].im = ct.im;
		}

	return(0);
	}

twid2(n, cwp, cx, cy)
int n;
struct complex *cwp, *cx, *cy;
	{
	int i;
	struct complex ctemp, cw;

	cw.re = cwp->re;
	cw.im = cwp->im;

	for (i=0; i<n; i++)
		{
		ctemp.re = rmul(cw,cy[i]);
		ctemp.im = imul(cw,cy[i]);
		cy[i].re = cx[i].re - ctemp.re;
		cy[i].im = cx[i].im - ctemp.im;
		cx[i].re += ctemp.re;
		cx[i].im += ctemp.im;
		}

	return(0);
	}