/**************************************************************************************
 FILE
 	realdft2.c - C source for an example inplementation of the DFT/IDFT
	of two N-point real sequences using one N-point complex DFT/IDFT.
 
 **************************************************************************************

 DESCRIPTION
 	This program is an example implementation of an efficent way of computing
	the DFT/IDFT of two real-valued sequences.

 	Assume we have two real-valued sequneces of length N - x1[n] abd x2[n].  The 
	DFT of x1[n] and x2[n] can be computed with one complex-valued DFT of length
	N, as shown above, by following this algorithm.

	1.  Form the complex-valued sequence x[n] from x1[n] and x2[n]
	      
		xr[n] = x1[n]  and   xi[n] = x2[n],   0,1, ..., N-1

	Note, if the sequences x1[n] and x2[n] are coming from another algorithm 
	or a data aquisition driver, this step may be eliminated if these put the
	data in the complex-valued format correctly. 

	2.  Compute X[k] = DFT{x[n]}

	This can be the direct form DFT algorithm, or an FFT algorithm.  If using an
	FFT lgorithm, make sure the output is in normal order - bitereversal is 
	performed.

	3.  Compute the following equations to get the DFTs of x1[n] and x2[n].

          	X1r[0] = Xr[0]
		X1i[0] = 0

		X2r[0] = Xi[0]
		X2i[0] = 0

		X1r[N/2] = Xr[N/2]
		X1i[N/2] = 0
		  
		X2r[N/2] = Xi[N/2]
		X2i[N/2] = 0
		  
		for k = 1,2,3, ...., N/2-1
		       	X1r[k] = (Xr[k] + Xr[N-k])/2        
		       	X1i[k] = (Xi[k] - Xi[N-k])/2
			X1r[N-k] = X1r[k]
			X1i[N-k] = X1i[k]

			X2r[k] = (Xi[k] + Xi[N-k])/2 
			X2i[k] = (Xr[N-k] - Xr[k])/2
			X2r[N-k] = X2r[k]
			X2i[N-k] = X2i[k]


	4. Form X[k] from X1[k] and X2[k]
	   
	   	for k = 0,1, ..., N-1		
			Xr[k] = X1r[k] - X2i[k]
			Xi[k] = X1i[k] + X2r[k]

	5.  Compute x[n] = IDFT{X[k]}

	This can be the direct form IDFT algorithm, or an IFFT algorithm.  If using 
	an IFFT	algorithm, make sure the output is in normal order - bitereversal is
	performed

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



#include <math.h>	/* include the C RTS math library		*/
#include "params2.h"	/* include file with parameters  		*/
#include "params.h"	/* include file with parameters  		*/


extern short x1[];
extern short x2[];

void dft(int, COMPLEX *);
extern void split2(int, COMPLEX *, COMPLEX *, COMPLEX *);


main()
{
	int	n, k;
	
	COMPLEX X1[NUMDATA];		/* array of real-valued DFT output sequence, X1(k) */	
	COMPLEX X2[NUMDATA];		/* array of real-valued DFT output sequence, X2(k) */
	COMPLEX x[NUMPOINTS+1];		/* array of complex DFT data, X(k)                 */
	
	

/* Forward DFT */

/* From the two N-point real sequnces, x1(n) and x2(n), form the N-point complex 
   sequence, x(n) = x1(n) + jx2(n) */

	for (n=0; n<NUMDATA; n++)
	{
		x[n].real = x1[n];
		x[n].imag = x2[n];
	}
	

/* Compute the DFT of x(n), X(k) = DFT{x(n)}.  Note, the DFT can be any
   DFT implementation such as FFTs. */

	dft(NUMPOINTS, x);

/* Because of the periodicitty property of the DFT, we know that X(N+k)=X(k). */

	x[NUMPOINTS].real = x[0].real;
 	x[NUMPOINTS].imag = x[0].imag;


/* The split function performs the additional computations required to get
   X1(k) and X2(k) from X(k). */

	split2(NUMPOINTS, x, X1, X2);


/* Inverse DFT - We now want to get back x1(n) and x2(n) from X1(k) and X2(k) using one
   complex DFT */


/* Recall that x(n) = x1(n) + jx2(n).  Since the DFT operator is linear,
   X(k) = X1(k) + jX2(k).  Thus we can express X(k) in terms of X1(k) and X2(k). */

	for (k=0; k<NUMPOINTS; k++)
	{
		x[k].real = X1[k].real - X2[k].imag;
		x[k].imag = X1[k].imag + X2[k].real;	
 	}

/* Take the inverse DFT of X(k) to get x(n).  Note the inverse DFT could be any
   IDFT implementation, such as an IFFT. */

/* The inverse DFT can be calculated by using the foward DFT algorithm directly
   by complex conjugation - x(n) = (1/N)(DFT{X*(k)})*, where * is the complex conjugate
   operator.  */

/* Compute the complex conjugate of X(k). */

	for (k=0; k<NUMPOINTS; k++)
	{
		x[k].imag = -x[k].imag;	
		
 	}

/* Compute the DFT of X*(k). */

	dft(NUMPOINTS, x);

/* Complex conjugate the output of the DFT and divide by N to get x(n). */

	for (n=0; n<NUMPOINTS; n++)
	{
		x[n].real = x[n].real/16;
		x[n].imag = (-x[n].imag)/16;	
	}

/* x1(n) is the real part of x(n), and x2(n) is the imaginary part of x(n). */

	for (n=0; n<NUMDATA; n++)
	{
		x1[n] = x[n].real;
		x2[n] = x[n].imag;
	}

	return(0);

 }