pfafft.c

该程序是用vc开发的对动态数组进行管理的DLL

/* Copyright (c) Colorado School of Mines, 2003.*/
/* All rights reserved.                       */

/*********************** self documentation **********************/
/*****************************************************************************
PFAFFT - Functions to perform Prime Factor (PFA) FFT's, in place

npfa		return valid n for fcomplex-to-fcomplex PFA
npfar		return valid n for real-to-fcomplex/fcomplex-to-real PFA
npfao		return optimal n for fcomplex-to-fcomplex PFA
npfaro		return optimal n for real-to-fcomplex/fcomplex-to-real PFA
pfacc		1D PFA fcomplex to fcomplex
pfacr		1D PFA fcomplex to real
pfarc		1D PFA real to fcomplex
pfamcc		multiple PFA fcomplex to real
pfa2cc		2D PFA fcomplex to fcomplex
pfa2cr		2D PFA fcomplex to real
pfa2rc		2D PFA real to fcomplex

*****************************************************************************
Function Prototypes:
int npfa (int nmin);
int npfao (int nmin, int nmax);
int npfar (int nmin);
int npfaro (int nmin, int nmax);
void pfacc (int isign, int n, fcomplex z[]);
void pfacr (int isign, int n, fcomplex cz[], float rz[]);
void pfarc (int isign, int n, float rz[], fcomplex cz[]);
void pfamcc (int isign, int n, int nt, int k, int kt, fcomplex z[]);
void pfa2cc (int isign, int idim, int n1, int n2, fcomplex z[]);
void pfa2cr (int isign, int idim, int n1, int n2, fcomplex cz[], float rz[]);
void pfa2rc (int isign, int idim, int n1, int n2, float rz[], fcomplex cz[]);

*****************************************************************************
npfa:
Input:
nmin		lower bound on returned value (see notes below)

Returned:	valid n for prime factor fft

******************************************************************************
npfao
Input:
nmin		lower bound on returned value (see notes below)
nmax		desired (but not guaranteed) upper bound on returned value

Returned:	valid n for prime factor fft

******************************************************************************
npfar
Input:
nmin		lower bound on returned value

Returned:	valid n for real-to-fcomplex/fcomplex-to-real prime factor fft

*****************************************************************************
npfaro:
Input:
nmin		lower bound on returned value
nmax		desired (but not guaranteed) upper bound on returned value

Returned:	valid n for real-to-fcomplex/fcomplex-to-real prime factor fft

******************************************************************************
pfacc:
Input:
isign		sign of isign is the sign of exponent in fourier kernel
n		length of transform (see notes below)
z		array[n] of fcomplex numbers to be transformed in place

Output:
z		array[n] of fcomplex numbers transformed

******************************************************************************
pfacr:
Input:
isign       sign of isign is the sign of exponent in fourier kernel
n           length of transform (see notes below)
cz          array[n/2+1] of fcomplex values (may be equivalenced to rz)

Output:
rz          array[n] of real values (may be equivalenced to cz)

******************************************************************************
pfarc:
Input:
isign       sign of isign is the sign of exponent in fourier kernel
n           length of transform; must be even (see notes below)
rz          array[n] of real values (may be equivalenced to cz)

Output:
cz          array[n/2+1] of fcomplex values (may be equivalenced to rz)

******************************************************************************
pfamcc:
Input:
isign       	sign of isign is the sign of exponent in fourier kernel
n           	number of fcomplex elements per transform (see notes below)
nt          	number of transforms
k           	stride in fcomplex elements within transforms
kt          	stride in fcomplex elements between transforms
z           	array of fcomplex elements to be transformed in place

Output:
z		array of fcomplex elements transformed

******************************************************************************
pfa2cc:
Input:
isign       	sign of isign is the sign of exponent in fourier kernel
idim        	dimension to transform, either 1 or 2 (see notes)
n1          	1st (fast) dimension of array to be transformed (see notes)
n2          	2nd (slow) dimension of array to be transformed (see notes)
z           	array[n2][n1] of fcomplex elements to be transformed in place

Output:
z		array[n2][n1] of fcomplex elements transformed

******************************************************************************
pfa2cr:
Input:
isign       sign of isign is the sign of exponent in fourier kernel
idim        dimension to transform, which must be either 1 or 2 (see notes)
n1          1st (fast) dimension of array to be transformed (see notes)
n2          2nd (slow) dimension of array to be transformed (see notes)
cz          array of fcomplex values (may be equivalenced to rz)

Output:
rz          array of real values (may be equivalenced to cz)

******************************************************************************
pfa2rc:
Input:
isign       sign of isign is the sign of exponent in fourier kernel
idim        dimension to transform, which must be either 1 or 2 (see notes)
n1          1st (fast) dimension of array to be transformed (see notes)
n2          2nd (slow) dimension of array to be transformed (see notes)
rz          array of real values (may be equivalenced to cz)

Output:
cz          array of fcomplex values (may be equivalenced to rz)

******************************************************************************
Notes:
Table of valid n and cost for prime factor fft.  For each n, cost
was estimated to be the inverse of the number of ffts done in 1 sec
on an IBM RISC System/6000 Model 320H, by Dave Hale, 08/04/91.
(Redone by Jack Cohen for 15 sec to rebuild NTAB table on advice of
David and Gregory Chudnovsky, 05/03/94).
Cost estimates are least accurate for very small n.  An alternative method
for estimating cost would be to count multiplies and adds, but this
method fails to account for the overlapping of multiplies and adds
that is possible on some computers, such as the IBM RS/6000 family.

npfa:
The returned n will be composed of mutually prime factors from
the set {2,3,4,5,7,8,9,11,13,16}.  Because n cannot exceed
720720 = 5*7*9*11*13*16, 720720 is returned if nmin exceeds 720720.

npfao:
The returned n will be composed of mutually prime factors from
the set {2,3,4,5,7,8,9,11,13,16}.  Because n cannot exceed
720720 = 5*7*9*11*13*16, 720720 is returned if nmin exceeds 720720.
If nmin does not exceed 720720, then the returned n will not be 
less than nmin.  The optimal n is chosen to minimize the estimated
cost of performing the fft, while satisfying the constraint, if
possible, that n not exceed nmax.

npfar and npfaro:
Current implemenations of real-to-fcomplex and fcomplex-to-real prime 
factor ffts require that the transform length n be even and that n/2 
be a valid length for a fcomplex-to-fcomplex prime factor fft.  The 
value returned by npfar satisfies these conditions.  Also, see notes 
for npfa.

pfacc:
n must be factorable into mutually prime factors taken 
from the set {2,3,4,5,7,8,9,11,13,16}.  in other words,
	n = 2**p * 3**q * 5**r * 7**s * 11**t * 13**u
where
	0 <= p <= 4,  0 <= q <= 2,  0 <= r,s,t,u <= 1
is required for pfa to yield meaningful results.  this
restriction implies that n is restricted to the range
	1 <= n <= 720720 (= 5*7*9*11*13*16)

pfacr:
Because pfacr uses pfacc to do most of the work, n must be even 
and n/2 must be a valid length for pfacc.  The simplest way to
obtain a valid n is via n = npfar(nmin).  A more optimal n can be 
obtained with npfaro.

pfarc:
Because pfarc uses pfacc to do most of the work, n must be even 
and n/2 must be a valid length for pfacc.  The simplest way to
obtain a valid n is via n = npfar(nmin).  A more optimal n can be 
obtained with npfaro.

pfamcc:
To perform a two-dimensional transform of an n1 by n2 fcomplex array 
(assuming that both n1 and n2 are valid "n"), stored with n1 fast 
and n2 slow:
    pfamcc(isign,n1,n2,1,n1,z); (to transform 1st dimension)
    pfamcc(isign,n2,n1,n1,1,z); (to transform 2nd dimension)

pfa2cc:
Only one (either the 1st or 2nd) dimension of the 2-D array is transformed.

If idim equals 1, then n2 transforms of n1 fcomplex elements are performed; 
else, if idim equals 2, then n1 transforms of n2 fcomplex elements are 
performed.

Although z appears in the argument list as a one-dimensional array,
z may be viewed as an n1 by n2 two-dimensional array:  z[n2][n1].

Valid n is computed via the "np" subroutines.

To perform a two-dimensional transform of an n1 by n2 fcomplex array 
(assuming that both n1 and n2 are valid "n"), stored with n1 fast 
and n2 slow:  pfa2cc(isign,1,n1,n2,z);  pfa2cc(isign,2,n1,n2,z);

pfa2cr:
If idim equals 1, then n2 transforms of n1/2+1 fcomplex elements to n1 real 
elements are performed; else, if idim equals 2, then n1 transforms of n2/2+1 
fcomplex elements to n2 real elements are performed.

Although rz appears in the argument list as a one-dimensional array,
rz may be viewed as an n1 by n2 two-dimensional array:  rz[n2][n1].  
Likewise, depending on idim, cz may be viewed as either an n1/2+1 by 
n2 or an n1 by n2/2+1 two-dimensional array of fcomplex elements.

Let n denote the transform length, either n1 or n2, depending on idim.
Because pfa2rc uses pfa2cc to do most of the work, n must be even 
and n/2 must be a valid length for pfa2cc.  The simplest way to
obtain a valid n is via n = npfar(nmin).  A more optimal n can be 
obtained with npfaro.

pfa2rc:
If idim equals 1, then n2 transforms of n1 real elements to n1/2+1 fcomplex 
elements are performed; else, if idim equals 2, then n1 transforms of n2 
real elements to n2/2+1 fcomplex elements are performed.

Although rz appears in the argument list as a one-dimensional array,
rz may be viewed as an n1 by n2 two-dimensional array:  rz[n2][n1].  
Likewise, depending on idim, cz may be viewed as either an n1/2+1 by 
n2 or an n1 by n2/2+1 two-dimensional array of fcomplex elements.

Let n denote the transform length, either n1 or n2, depending on idim.
Because pfa2rc uses pfa2cc to do most of the work, n must be even 
and n/2 must be a valid length for pfa2cc.  The simplest way to
obtain a valid n is via n = npfar(nmin).  A more optimal n can be 
obtained with npfaro.

******************************************************************************
References:  
Temperton, C., 1985, Implementation of a self-sorting
in-place prime factor fft algorithm:  Journal of
Computational Physics, v. 58, p. 283-299.

Temperton, C., 1988, A new set of minimum-add rotated
rotated dft modules: Journal of Computational Physics,
v. 75, p. 190-198.

Differ significantly from Press et al, 1988, Numerical Recipes in C, p. 417.

******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 04/27/89
*****************************************************************************/
/**************** end self doc ********************************/

#include "cwp.h"

#define NTAB 240
static struct {
int n;  float c;
} nctab[NTAB] = {
{       1, 0.000052f },
{       2, 0.000061f },
{       3, 0.000030f },
{       4, 0.000053f },
{       5, 0.000066f },
{       6, 0.000067f },
{       7, 0.000071f },
{       8, 0.000062f },
{       9, 0.000079f },
{      10, 0.000080f },
{      11, 0.000052f },
{      12, 0.000069f },
{      13, 0.000103f },
{      14, 0.000123f },
{      15, 0.000050f },
{      16, 0.000086f },
{      18, 0.000108f },
{      20, 0.000101f },
{      21, 0.000098f },
{      22, 0.000135f },
{      24, 0.000090f },
{      26, 0.000165f },
{      28, 0.000084f },
{      30, 0.000132f },
{      33, 0.000158f },
{      35, 0.000138f },
{      36, 0.000147f },
{      39, 0.000207f },
{      40, 0.000156f },
{      42, 0.000158f },
{      44, 0.000176f },
{      45, 0.000171f },
{      48, 0.000185f },
{      52, 0.000227f },
{      55, 0.000242f },
{      56, 0.000194f },
{      60, 0.000215f },
{      63, 0.000233f },
{      65, 0.000288f },
{      66, 0.000271f },
{      70, 0.000248f },
{      72, 0.000247f },
{      77, 0.000285f },
{      78, 0.000395f },
{      80, 0.000285f },
{      84, 0.000209f },
{      88, 0.000332f },
{      90, 0.000321f },
{      91, 0.000372f },
{      99, 0.000400f },
{     104, 0.000391f },
{     105, 0.000358f },
{     110, 0.000440f },
{     112, 0.000367f },
{     117, 0.000494f },
{     120, 0.000413f },
{     126, 0.000424f },
{     130, 0.000549f },
{     132, 0.000480f },
{     140, 0.000450f },
{     143, 0.000637f },
{     144, 0.000497f },
{     154, 0.000590f },
{     156, 0.000626f },
{     165, 0.000654f },
{     168, 0.000536f },
{     176, 0.000656f },
{     180, 0.000611f },
{     182, 0.000730f },
{     195, 0.000839f },
{     198, 0.000786f },
{     208, 0.000835f },
{     210, 0.000751f },
{     220, 0.000826f },
{     231, 0.000926f },
{     234, 0.000991f },
{     240, 0.000852f },
{     252, 0.000820f },
{     260, 0.001053f },
{     264, 0.000987f },
{     273, 0.001152f },
{     280, 0.000952f },
{     286, 0.001299f },
{     308, 0.001155f },
{     312, 0.001270f },
{     315, 0.001156f },
{     330, 0.001397f },
{     336, 0.001173f },
{     360, 0.001259f },
{     364, 0.001471f },
{     385, 0.001569f },
{     390, 0.001767f },
{     396, 0.001552f },
{     420, 0.001516f },
{     429, 0.002015f },
{     440, 0.001748f },
{     455, 0.001988f },
{     462, 0.001921f },
{     468, 0.001956f },
{     495, 0.002106f },
{     504, 0.001769f },
{     520, 0.002196f },
{     528, 0.002127f },
{     546, 0.002454f },
{     560, 0.002099f },
{     572, 0.002632f },
{     585, 0.002665f },
{     616, 0.002397f },
{     624, 0.002711f },
{     630, 0.002496f },
{     660, 0.002812f },
{     693, 0.002949f },
{     715, 0.003571f },
{     720, 0.002783f },
{     728, 0.003060f },
{     770, 0.003392f },
{     780, 0.003553f },
{     792, 0.003198f },
{     819, 0.003726f },
{     840, 0.003234f },
{     858, 0.004354f },
{     880, 0.003800f },
{     910, 0.004304f },
{     924, 0.003975f },
{     936, 0.004123f },
{     990, 0.004517f },
{    1001, 0.005066f },
{    1008, 0.003902f },
{    1040, 0.004785f },
{    1092, 0.005017f },
{    1144, 0.005599f },
{    1155, 0.005380f },
{    1170, 0.005730f },
{    1232, 0.005323f },
{    1260, 0.005112f },
{    1287, 0.006658f },
{    1320, 0.005974f },
{    1365, 0.006781f },
{    1386, 0.006413f },
{    1430, 0.007622f },
{    1456, 0.006679f },
{    1540, 0.007032f },
{    1560, 0.007538f },
{    1584, 0.007126f },
{    1638, 0.007979f },
{    1680, 0.007225f },
{    1716, 0.008961f },
{    1820, 0.008818f },
{    1848, 0.008427f },
{    1872, 0.009004f },
{    1980, 0.009398f },
{    2002, 0.010830f },
{    2145, 0.012010f },
{    2184, 0.010586f },
{    2288, 0.012058f },
{    2310, 0.011673f },
{    2340, 0.011700f },
{    2520, 0.011062f },
{    2574, 0.014313f },
{    2640, 0.013021f },
{    2730, 0.014606f },
{    2772, 0.013216f },
{    2860, 0.015789f },
{    3003, 0.016988f },
{    3080, 0.014911f },
{    3120, 0.016393f },
{    3276, 0.016741f },
{    3432, 0.018821f },
{    3465, 0.018138f },
{    3640, 0.018892f },
{    3696, 0.018634f },
{    3960, 0.020216f },
{    4004, 0.022455f },
{    4095, 0.022523f },
{    4290, 0.026087f },
{    4368, 0.023474f },
{    4620, 0.024590f },
{    4680, 0.025641f },
{    5005, 0.030303f },
{    5040, 0.025253f },
{    5148, 0.030364f },
{    5460, 0.031250f },
{    5544, 0.029412f },
{    5720, 0.034404f },
{    6006, 0.037500f },
{    6160, 0.034091f },
{    6435, 0.040214f },
{    6552, 0.037221f },
{    6864, 0.042735f },
{    6930, 0.040214f },
{    7280, 0.042980f },
{    7920, 0.045872f },
{    8008, 0.049505f },
{    8190, 0.049834f },
{    8580, 0.055762f },
{    9009, 0.057034f },
{    9240, 0.054945f },
{    9360, 0.056818f },
{   10010, 0.066667f },