# Fft.h

Jump to navigation
Jump to search

/* * 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_