Fft.h

/*
 * fft.h
 *
 *  Created on: Jan 21, 2016
 *  Author: mryan
 */
#ifndef FFT_H_
#define FFT_H_

#include <math.h>
#include <time.h>
#include <unistd.h>

// assuming basic number type is float (for now)
# ifndef fft_scalar
#   define fft_scalar float
# endif

// (in-phase, quadrature) data is treated as a complex number
typedef struct {
    fft_scalar re;					// real, i.e. in-phase
    fft_scalar im;					// imaginary, i.e. quadrature
}fft_cpx;

// each CPU has a struct of this type set up in SRAM
typedef struct {
    unsigned int status;
    unsigned int cinadrs;
    unsigned int lenP2;				// power of 2 that gives length
    unsigned int begin;
    unsigned int nextbegin;
    unsigned int oddadd;
    unsigned int startP2;			// power of 2 that gives butterfly to start with
    unsigned int twadrs;
    unsigned int toffset;
    unsigned int tinc;
    unsigned int normalise;
    unsigned int arrange;
    unsigned int realonly;
    unsigned int secondpass;
    unsigned int spare1;			// bring up to a 16 word boundary
    unsigned int spare2;
}fft_cpu;

#define FFT_WORK_MASK   1			// pick off last bit
#define FFT_NOWORK		0			// nothing to be done/finished
#define FFT_START		1			// work to be done/working

#define MAX_CPUS 		4
#define LAST_CPU		8			// 0x1000
#define START_CPU_INDEX	0			// N.B. fft_multi() assumes is 0, CPU that starts on power up

#define IRQMP_ADRS		0x80000210  // address of multiprocessor status register
#define CPUS_PM1_MASK	0xf0000000	// tells how many CPUs in system minus 1
#define CPUS_PM1_SHIFTS 28			// shift them to low order
#define	CPUS_MASK    	0xf			// bits used to power-up CPUs
#define CPU_START_PAUSE	1000		// pause between and after CPU startups

// set up space for the various arrays. S698PM SRAM is 32 Mbyte from 0x40000000.
#define FFT_SRAM 	 	0x60000000				// half way
#define FFT_FLAGS	 	FFT_SRAM+0x80000
#define FFT_CIN  	 	FFT_SRAM+0x100000		// space for 2 power 20 bytes,
#define FFT_RAIN	 	FFT_SRAM+0x200000		// i.e. 2 power 18 doubles
#define FFT_TWIDDLES 	FFT_SRAM+0x400000		// i.e. 2 power 17 complex nos
#define FFT_COUT	 	FFT_SRAM+0x600000
#define FFT_MAP			FFT_SRAM+0x800000		// maps input index to rearranged index

#define FFT_MIN_LENP2	5			// minimum input 32 (in-phase,quadrature)
#define FFT_MAX_LENP2	17			// for consistency with space allocation above

/*
 * FUNCTION DECLARATIONS
 */

/*
 *  The functions declared below are used to create an FFT
 *  using 1, 2 or 4 processors acting in parallel.
 *
 *	At present the input array length must be a power of 2.
 *
 *  At present the input data in each (in-phase, quadrature) pair must be C floats.
 *
 *	The processing is based on arrays of complex numbers corresponding to
 *	the (in-phase, quadrature) data.
 *
 *	If the input data is in-phase only, the input array will have 0 in each entry's
 *	quadrature component. An option is provided to allow this be exploited so as
 *	to get the result roughly twice as fast
 *
 *  Usage:
 *		0.	Create an array holding the in-phase and quadrature data
 *			as complex numbers (fft_cpx).
 *
 *			At present, the length of this array must be a power of 2.
 *			At present, data elements should be of the C float type.
 *
 *		1.	The complex array must be arranged in an order suitable
 *			for 'discrimination in time' (DIT) processing.
 *
 *			This can be done as the entries are made one at a time.
 *			An index array can first be created for use to map
 *			successive input to the correct position in the
 *			array. This can be stored in ROM, and reused to allocate data
 *			as it comes in to the appropriate place in the array.
 *
 *			The function fftArrangeIndex() below is provided for
 *			this purpose.
 *
 *			If the input data is not rearranged as it arrives,
 *			a parameter can be set to cause rearrangement to be done
 *			in parallel in the first phase of FFT processing on each CPU.
 *			In that case the FFT result does not overwrite the input data,
 *			but can be found in the FFT_RAIN area defined below.
 *
 *			Alternatively, if the array already exists (as an array of complex numbers),
 *			it can be rearranged in situ using fftArrangeIn(),
 *			or as a new complex array using fftArrangOut().
 *
 *			The data in the rearranged array is overwritten by the FFT
 *			results, so if it is desired to keep the original data
 *			either the parallel rearrangement parameter should be set
 *			or the fftArrangeOut() function should be used.
 *
 *			To avoid the overhead of rearrangement it is best to rearrange
 *			the data as it arrives, if possible.
 *
 *		3.	The 'twiddles' array, containing complex roots of unity,
 *			must be created. This is half the length of the input
 *			array of complex numbers, as due to symmetry only half
 *			of the roots need to be created and stored.
 *			The function fillTwiddles() is provided for this purpose.
 *
 *		N.B.	If the size of the data array does not change,
 *				the twiddles array can be stored in ROM.
 *				Two versions are needed in that case, one for the forward FFT,
 *				one for the inverse FFT. Both can be generated using fillTwiddles().
 *
 *		4.	The appropriately arranged complex input array is passed to fftMulti(),
 *			together with the twiddles array and which processors to use.
 *			The result overwrites the original data in the rearranged array.
 *
 *		5.	Other than input rearrangement, no copying is involved in passing data
 *			to the different processors, each is given a reference to a contiguous
 *			portion of the input array.
 *			A reference to the twiddles array is also passed, though the entries
 *			used in this by a processor are not contiguous.
 *
 *		6.	All stages of processing up to the last two (4 CPUs) or last one (2 CPUs)
 *			are done independently on the CPUs, with no interaction between them.
 *			The last stage or stages require results generated by a different CPU,
 *			so must wait for the previous stages to complete before progressing.
 *
 *		7.	If the 'real only', i.e in-phase only, option is selected any 'imaginary',
 *			i.e quadrature, component in the input is ignored and the input array compressed
 *			into half its length with the entries in the remaining half set to (0.0,0.0).
 *			The processing is done roughly twice as fast.
 *			Only the positive frequency components are produced,
 *			the negative frequency components are easily derived if needed by taking
 *			the complex conjugates of the positive frequency components,
 *			i.e. changing the sign of the imaginary component and leave the rest the same.
 */

/*
 * FUNCTION DECLARATIONS
 */
/*
 * This ASM-function returns the CPU index (0 to 3) of the current S698PM CPU.
 */
 static inline unsigned int get_asr17(void);


/*
 * Starting CPUs
 *
 * All CPUs have the same entry point.
 * The linker defaults to 0x40000000.
 * DMON's ep command is used to set each CPU's initial PC to this.
 *
 * This function is entered immediately on entering main().
 * It determines the index (0 to 3) of the CPU running the code.
 * If this is the START_CPU_INDEX it initializes the data areas used by all CPUs
 * and returns to main().
 * Otherwise, it checks to see if any work is waiting, does it,
 * and halts the CPU. When work becomes available, the CPU is restarted,
 * does the work, and halts again. This function never returns if running on
 * these CPUs.
 *
 */
int fftCpusStart(void);


/*
 * Preliminary start of the selected CPUs if other than the one running this.
 * Each such CPU starts, runs through its initialization code, and powers down
 * Each CPU shares the same initialization code, so there is a pause between
 * each startup just in case ...
 *
 * A CPU is restarted when there is work for it to do
 * (restarts are done directly, not by this).
 *
 *  Inputs:
 *  	const int		which cpus		from 0 to 0xf, the current CPU, if selected, is ignored
 *
 *  Returns:							None
 *
 * 						N.B. This routine should not be used to restart CPUs
 */
void startCpus(const int whichcpus);


/*
 * INPUT ARRAY REARRANGEMENT
 *
 * This is used with a discrimination in time (DIT) implementation of the FFT.
 * Separate parts of the rearranged array can be processed independently
 * on different CPUs without any interaction between the CPUs being required.
 * The last two stages (4 CPUS) or just the last stage (2 CPUS)can only be done when
 * the previous stages are finished, and are again distributed across all CPUs.
 *
 * The best approach to input rearrangement is to do it as the data arrives.
 *
 * Alternatively, by setting a parameter rearrangement can be done in parallel
 * one each CPU.
 *
 * Three very similar functions, not parallel, are given to allow a complete input
 * array be rearranged:
 * 	fftArrangeOut		an out of place rearrange, the input array is unchanged
 * 	fftArrangeIn		an in place rearrange, the input array is rearranged in situ
 * 	fftaArrangeIndex	creates an array of indexes. Entry [i] gives the index where
 * 						the (in-phase, quadrature) complex pair should be stored when
 * 						data is read one pair at a time
 *						- the input array will then not need further rearrangement
 *
 * All three functions assume the array length is a power of 2, and this power is what
 * is input to give the array length.
 *
 * (Thanks for bit-swap method used to graphics.stanford.edu/~seander/bithacks.html#ReverseParallel)
 *
 */

/*
 *	Out of place rearrangement of input array to prepare for FFT.
 *	Assumes array length is a power of 2.
 *
 *	Inputs:
 *		fft_cpx* 		cin		pointer to an array of complex numbers
 *								(the length of this array must be a power of 2)
 *		int				lenP2	log to base 2 of array length. 0 <= lenP2 <= 30.
 *
 *		fft_cpx* 		result	pointer to resulting array of complex numbers
 *								(same length as cin array)
 *
 *	Result:				A rearranged version of cin input array in the result array
 *
 *	N.B.				Space for the result array of complex numbers must be allocated
 *						before calling this function, same length as the cin input array.
 *
 */
void fftArrangeOut(fft_cpx* cin, fft_cpx* result, int lenP2);


 /*
  *	In place rearrangement of input array to prepare for FFT.
  *	Assumes array length is a power of 2.
  *
  *	Inputs:
  *		fft_cpx* 		cin		pointer to an array of complex numbers
  *								(the length of this array must be a power of 2)
  *		int				lenP2	log to base 2 of array length.
  *
  *	Result:				A rearranged version of cin input array in the same space
  *
  */
 void fftArrangeIn(fft_cpx* x, int lenP2);


 /*
  *	Generate an array of indices giving the array position where each successive input
  *	should be placed in the input array.
  *	This array could be stored in ROM and used when inputting data so as to avoid
  *	any subsequent need to spend time rearranging the input array.
  *
  *	Assumes amount of input datais a power of 2.
  *
  *	Inputs:
  *		int * 			mapindex	pointer to an array of ints that will hold the index to be used
  *									in the complex input array for each successive data input.
  *									(the length of this array must be a power of 2
  *		int				lenP2		log to base 2 of array length. 0 <= lenP2 <= 30
  *									(regarding the array as an array of complex numbers)
  *
  *	Result:				mapindex will hold appropriate index in data array to use to hold complex input i
  *
  *	N.B. if the input data is an array of scalars, the length of mapindex[] will be half of the length
  *		 of the array to hold the scalars, and the input index k used to look up the mapindex array should only
  *	     increment after each pair of reals have been added to dataarray[mapindex[k]<< 1] and dataarray[(mapindex[k]<< 1) +1]
  *	     (they will be treated as the real and imaginary parts of the complex number at the index given by mapindex[k])
  *
  *	N.B. Space for the mapindex array of ints must be allocated before calling this function,
  *		 same length as complex data input array, half length of real data input arrays.
  *
  *	N.B. The mapping array is its own inverse
  *
  */
 void fftArrangeIndex(int* mapindex, int lenP2);


 /*
   *	If the 'realonly' option is selected, the input data is treated as having in-phase components only
   *	and any quadrature component is ignored and will be overwritten in processing.
   *	The 'real', i.e. in-phase input components are rearranged into half the number of inputs
   *	with every second in-phase input now in the quadrature position, and zero in remaining half of the input array.
   *	This allows faster FFT processing but requires post-processing.
   *	(Post processing is distributed across all processors and done later in a final processing phase.
   *	 It produces the positive frequency components only, the negative components can be obtained by taking
   *	 the complex conjugates of these.)
   *
   *	Assumes amount of input data is a power of 2.
   *
   *	Inputs:
   *		fft_cpx * 		cin			input array with all 0s in quadrature (imaginary) positions.
   *									Altered so that quadrature positions hold alternate real inputs,
   *									and remaining array entries are made (0.0,0.0)
   *		int				lenP2		log to base 2 of array length. 0 <= lenP2 <= 30.
   *									(regarding the array as an array of complex numbers)
   *
   *	Result:				first half of cin holds rearranged data, rest is (0.0,0.0)
   *
   *						N.B. 		use of this function can usually be avoided by simply casting an input
   *									array of scalars to be an array of complex numbers.
   */
 void fftRealOnlyArrange(fft_cpx * cin, const int lenP2);


 /*
  * TWIDDLES ARRAY
  */

 /*
  * For an input data array of length N the 'twiddle' factors are the N complex roots
  * of unity. Only N/2 of these are calculated, as the others are just minus these,
  * and are provided by using subtraction rather than addition.
  *
  * Getting the complex roots of unity is time consuming, so where the size of
  * array to be processed is known in advance the twiddles array should be stored
  * in ROM. Two arrays should be stored, one for the forward FFT, one for the inverse.
  * Each has length N/2, half that of the complex data array.
  *
  * Time consuming calculations of sines and cosines are involved, but in some cases
  * can be avoided by exploiting various symmetries. These require the complex data array
  * size to be at least 8, but it is restricted here to at least 1 << FFT_MIN_LENP2, or 32 at present.
  *
  *
  * Inputs:
  *		fft_cpx* 		twiddles	pointer to an array of complex numbers
  *									(the length of this array will be half the length
  *									 of the complex data array with which the twiddles
  *									 array will be used)
  *		const int		lenP2		log to base 2 of complex data array length.
  *									(from FFT_MIN_LENP2	to FFT_MAX_LENP2, i.e. 5 to	17)
  *
  *		const int		type 		0 for FFT, 1 for inverse FFT
  *
  * Result:				N/2 of the N complex roots of unity in twiddles array
  *
  * N.B.				Space for the twiddles array of N/2 complex numbers must be allocated
  *						before calling this function, where N is the length of the complex data
  *						array. i.e. 1 << lenP2.
  *
  */
 void fillTwiddles(fft_cpx* twiddles, const int lenP2, const int type);


 /*
  * FFT PROCESSING
  *
  * The (in-phase, quadrature) input data is regarded as an array of complex numbers.
  * The resulting FFT is also an array of complex numbers.
  *
  * If the data is not already in rearranged order, by setting the 'arrange' parameter
  * the input array can be rearranged in parallel, with each CPU doing its own section.
  * This involves copying the original input data into the area at FFT_RAIN,
  * where it will eventually be overwritten by the FFT result, with the original
  * input data left unchanged.
  *
  */

 /*
  * This function calculates all or part of the FFT.
  * It can be used on separate processors or in separate threads.
  *
  * It processes part or all of a complex data array, replacing parts of the
  * array with appropriate FFT data. After any initial rearrangement
  * all processing is done in situ .
  *
  * The length of a sub-array must be a power of 2 and no less than 8.
  * The sub-array can be the whole array.
  *
  * The size of the butterfly blocks with which to start processing a sub-array is
  * input as a power of 2.
  *
  * Where the blocks in the subarray have not already been processed,
  * the starting power of 2 is 1.
  *
  * If parts of a subarray have already been processed, the correct butterfly
  * step must be given. All parts of the subarray must have the same size.
  *
  * The post-processing phrase needed when the original data is in-phase only
  * is done if 'secondpass' is set.
  *
  * Rearrangement of the data for this section is done if the 'arrange' is set.
  *
  *
  * Inputs:
  *		fft_cpx* 		arranged	pointer to the array of complex numbers that
  *									initially contains the rearranged complex data input.
  *									Will hold result.
  *		const int		lenP2		Power of 2 that gives length of input array.
  *
  *		const int		start		position in array at which sub-array starts
  *		const int		nextstart	position one after end of data to process
  *									N.B. (nextstart - start) must be a power of 2
  *
  *		const int		oddadd		Added internally to half of sub-array size.
  *									Used to get distance from even to odd entries.
  *
  *		const int		startP2		power of 2 that gives butterfly at which to start.
  *									Must be at least 1.
  *									1 if no previous processing done on this sub-array.
  *									Otherwise gives size of butterfly to start with.
  *
  *		const fft_cpx*	twiddles	pointer to twiddles array - this must correspond
  *									to the original complete complex input array.
  *		const int		toffset		initial offset into twiddles array
  *		const int		twinc		initial increment to using going through twiddles array
  *
  *		const int		normalise	0 for don't normalize, anything else -> normalize
  *
  *		const int		arrange		0 for don't arrange, anything else -> arrange
  *
  *		const int		realonly	0 if original data complex, 1 if in-phase only
  *									and has been rearranged
  *
  *		const int		secondpass	0 if not ready to do post-processing,
  *									1 to do post-processing only
  *
  * Result:				FFT processed values in situ in this subarray
  *
  *	N.B.				If 'realonly' is selected, assumes that the input array is a compressed
  *						version of something twice as long, and that this additional space
  *						is available.
  *
  *	N.B.				If 'realonly' is selected, FFT result contains the positive frequencies
  *						only. Negative frequencies are easily derived if needed as the complex
  *						conjugates of these.
  *
  */
 void fftProc(fft_cpx* arranged, const int lenP2,
 			 const int start, const int nextstart,
 			 const int oddadd, const int startP2,
     		 const fft_cpx* twiddles,
     		 const int toffset, const int twinc,
     		 const int normalise, const int arrange,
     		 const int realonly, const int secondpass);


 /*
  * Allocates work of doing FFT to one or more CPUs
  * and combines results from them to give FFT.
  *
  * The input is an array of complex numbers representing the
  * (in-phase,quadrature) data.
  * Its length must be a power of 2 - this power rather than the
  * length is input.
  *
  * The complex number array giving the FFT is returned in situ
  * in the input array, unless the 'arrange parameter is set,
  * in which case the result is in the area at FFT_RAIN.
  *
  * At present the number of CPUs to use must be 1, 2 or 4.
  *
  * Inputs:
  *		fft_cpx* 		arranged	pointer to an array of complex numbers that
  *									holds the rearranged complex data input.
  *									Will hold result.
  *		const int		lenP2		power of 2 that gives size of array
  *
  *		const fft_cpx*	twiddles	pointer to twiddles array
  *									- this must correspond to input array
  *									  (it will be half its size)
  *
  *		const int		whichcpus	1 to 0xf, which CPUs to use (binary xxxx, x=1 for use)
  *
  *		const int		normalise	0 for don't normalize, anything else -> normalize
  *
  *		const int		arrange		0 for don't arrange, anything else -> arrange
  *
  *		const int		realonly	0 if original data has both components, 1 if in-phase only
  *									and has been rearranged
  *
  * Result:				Complex valued FFT in situ in arranged
  *
  * N.B.				Assumes START_CPU_INDEX is 0
  *
  */
 void fftMulti(fft_cpx* cin, const int lenP2,
 			  const fft_cpx* twiddles,
     		  const int whichcpus,
     		  const int normalise,
     		  const int arrange,
     		  const int realonly);


 /*
  * Gives one part of the FFT calculation to the selected CPU for execution.
  *
  * The input is part of an array of complex numbers.
  * Its length must be a power of 2 - this power rather than the
  * length is input.
  *
  * The output is placed in the appropriate locations in the input array,
  * unless the 'arrange' parameter is set, in which case the output is
  * in the area at FFT_RAIN.
  *
  * At present the identifier of the CPU to use must be 1, 2, 4 or 8, anything
  * else the function returns without doing anything.
  *
  * Inputs:
  *		const int		which		which CPU to use - must be 1, 2, 4 or 8
  *
  *		fft_cpx* 		arranged	pointer to an array of complex numbers that
  *									holds the rearranged complex data input.
  *									Will hold result.
  *		const int		lenP2		power of 2 that gives size of array
  *
  *		const int		begin		index into arranged, location where this processing is to start
  *
  *		const int		nextbegin	index into arranged, begin <= index < nextbegin
  *
  *		const int		oddadd		Added internally to half of sub-array size.
  *									Used to get distance from even to odd entries.
  * 	const int		startP2		1 << (lenP2 - startP2) gives initial step size in twiddles array
  *
  *		const fft_cpx*	twiddles	pointer to twiddles array
  *									- this must correspond to input array
  *									  (it will be half its size)
  *
  *  	const int		toffset		initial offset into twiddles array
  *		const int		twinc		initial increment to using going through twiddles array
  *
  *		const int		normalise	0 for don't normalize, anything else -> normalize
  *
  *		const int		arrange		0 for don't arrange, anything else -> arrange
  *
  * 	const int		realonly	0 if original data complex, 1 if in-phase only
  *									and has been rearranged
  *
  * Result:				Updated complex values for this stage of FFT processing in arranged
  *
  */
 void giveToCpu(const int which,							// 1,2,4 or 8
 			   fft_cpx* cin, const int lenP2,				// input array
 			   const int begin, const int nextbegin,		// part indices
 			   const int oddadd, const int startP2,
 			   const fft_cpx* twiddles,
 			   const int toffset, const int twinc,
 			   const int normalise, const int arrange,
 			   const int realonly, const int secondpass);


 /*
  * Wait for selected CPUs to have finished.
  * Returns when no CPU has work pending.
  * Can wait indefinitely. TODO put in a time check
  *
  * All stages of the FFT except at most the last two (for the S698PM) proceed independently
  * with no need for communication between CPUs.
  *
  * If using 4 CPUs, each of the last two stages must wait for all CPUs to have finished
  * the preceding stage.
  *
  * If using 2 CPUs, they must wait to process the last stage until both CPUs have finished
  * the earlier stages.
  *
  * Inputs:
  *		const int		whichcpus	which CPUs to check - xxxx, 1 to check, 0 not to check
  *
  *
  * Result:				Wait, only return when no selected CPU has work pending
  */
 void fftWait();


/*
 * UTILITIES
 */

 /*
  * Get the number of CPUs enabled by the selections in whichcpus.
  *
  * 	Input:
  * 		const int 	whichcpus	0 to 0xf  (binary xxxx, 1 = use CPU
  *
  * 	Output:
  * 		int			Count of currently selected CPUs ( 0 to 4)
  *
  */
 int getNumCpus(const int whichcpus);


 /*
  * Validate the input settings and correct if possible
  *
  *	Inputs:
  *		const int		lenP2			must be from FFT_MIN_LENP2	to FFT_MAX_LENP2, (i.e. 5 to 17)
  *
  *		int *			cpusrequest		should be from 0 to 0xf, corresponding to 1, 2 or 4 cpus (but not 3)
  *										(correction attempt made if three CPUs chosen)
  *
  *		const int		numruns			no validation
  *
  *		const int		type			must be 0 or 1
  *
  *		int *			lines			corrected to 1 << lenP2 if greater than this
  *
  *	Returns:
  *		int				0 if all o.k., 1 if a problem (messages given),
  *
  * Inputs:
  *
  */
 int settingsIncorrect(const int lenP2, int *cpusrequest, const int numruns, const int type, int* lines);


 /*
  * Display an array of complex numbers on the console.
  * Can be used to display data input or FFT result.
  *
  * N.B. At present assumes floats being used
  *
  * Inputs:
  * 	fft_cpx* 		in			input array of complex numbers
  *
  *		const int		length		length of the array
  *
  *		const double	tolerance	absolute values less than this are shown as 0.0
  *
  *	Result:				Array is output to the console
  *
  */
 void showC(const fft_cpx* in, const int length, const double tolerance);


 /*
   * Display main parts of FFT result on the console.
   * Shows initial 8 positive frequency components including 0 component,
   * show 8 components before Nyquist frequency, 8 components from Nyquist,
   * and last 8 components.
   * If 'realonly' is selected, shows result up to Nyquist frequency only.
   *
   * N.B. At present assumes floats being used
   *
   * Inputs:
   * 		fft_cpx* 		in			input array of complex numbers
   *
   *		const int		length		length of the array
   *
   *		const double	tolerance	absolute values less than this are shown as 0.0
   *
   *		const int 		realonly	0 to show full FFT, 1 to show non-negative frequencies only.
   *
   *	Result:				Array is output to the console
   *
   */
 void show8C(const fft_cpx* in, const int length, const double tolerance, const int realonly);


#endif // FFT_H_