SIAMR EVIEW SocietyforIndustrialandAppliedMathematics Vol

SIAMR EVIEW SocietyforIndustrialandAppliedMathematics Vol - Description

46No3pp443454 Accelerating the Nonuniform Fast Fourier Transform Leslie Greengard JuneYub Lee Abstract ThenonequispacedFouriertransformarisesinavarietyofapplicationareasfrommedical imagingtoradioastronomytothenumericalsolutionofpartialdi64256erential ID: 24427 Download Pdf

79K - views

SIAMR EVIEW SocietyforIndustrialandAppliedMathematics Vol

46No3pp443454 Accelerating the Nonuniform Fast Fourier Transform Leslie Greengard JuneYub Lee Abstract ThenonequispacedFouriertransformarisesinavarietyofapplicationareasfrommedical imagingtoradioastronomytothenumericalsolutionofpartialdi64256erential

Similar presentations


Download Pdf

SIAMR EVIEW SocietyforIndustrialandAppliedMathematics Vol




Download Pdf - The PPT/PDF document "SIAMR EVIEW SocietyforIndustrialandAppli..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.



Presentation on theme: "SIAMR EVIEW SocietyforIndustrialandAppliedMathematics Vol"— Presentation transcript:


Page 1
SIAMR EVIEW 2004SocietyforIndustrialandAppliedMathematics Vol.46,No.3,pp.443–454 Accelerating the Nonuniform Fast Fourier Transform Leslie Greengard June-Yub Lee Abstract. ThenonequispacedFouriertransformarisesinavarietyofapplicationareas,frommedical imagingtoradioastronomytothenumericalsolutionofpartialdifferentialequations. In atypicalproblem,oneisgivenanirregularsamplingof datainthefrequencydomain andoneisinterestedinreconstructingthecorrespondingfunctioninthephysicaldomain. When the sampling is uniform, the fast Fourier transform (FFT) allows this calculation

tobecomputedin log )operationsratherthan )operations. Unfortunately, when the sampling is nonuniform, the FFT does not apply. Over the last few years, a number of algorithms have been developed to overcome this limitation and are often referredtoas nonuniformFFTs (NUFFTs). Theserelyonamixtureofinterpolationand thejudicioususeoftheFFTonanoversampledgrid[A.Duttand(.)o*hlin, SIAMJ. Sci.Comput. ,14(1993),pp.1368113832. Inthispaper,weobservethatoneofthestandardinterpolationor3gridding4schemes, based on Gaussians, can be accelerated by a signi6cant factor without precomputation and storage of the

interpolation weights. This is of particular value in two7 and three7 dimensional settings, saving either 18 in storage in dimensions or a factor of about 9118inC;Utime(independentofdimension). Key words. nonuniformfastFouriertransform,fastgridding,FFT,imagereconstruction AMS subject classifications. 42A38,44A39,69T98,69)18 DOI. 18.113=/S883614498343288@ 1. Introduction. In this note, we describe an extremely simple and efficient implementation of the nonuniform fast Fourier transform (NUFFT). There are a host of applications of such algorithms, and we refer the reader to the

references [2, 6, 8, 11, 13, 14, 17( for examples. )e restrict our attention here to one: function (or image) reconstruction from Fourier data as discussed in [6, 8, 11, 14(. Let us begin, howe,er, with a more precise description of the computational tas-. In two dimensions, we de.ne the nonuniform discrete Fourier transform of types 1 and 2 according to the formulae ,k )= =0 ,k (1) )eceived by the editors July 23, 2883B accepted for publication (in revised form) December 1, 2883BpublishedelectronicallyJuly38,2884. httpC//www.siam.org/journals/sirev/4673/43288.html Courant Institute of

Mathematical Sciences, New Eor* University, New Eor*, NE 18812 (greengardFcims.nyu.edu). The wor* of this author was supported by the Applied Mathematical Sciences;rogramoftheU.S.DepartmentofGnergyundercontractDGFGO288G)29893. Department of Mathematics, Gwha Womans University, Seoul, 1287=98, Korea (jyleeFmath. ewha.ac.*r). Thewor*ofthisauthorwassupportedbytheKorea)esearchFoundationundergrant 288278197C;8844. 443
Page 2
444 LESLIEGREENGARDANDJ,NE-Y,BLEE )= ,k ,k (2) respecti,ely, where [0 [0 ( and ,k It is, perhaps, con,enient to thin- of (1) as a discreti1ation of the Fourier

integral ,k )= (2 ,k (3) with ser,ing as the discreti1ation points. If we let denote the 2uadrature weight corresponding to , then we obtain (1) by setting . E2uation (2), of course, is simply the e,aluation of a .nite Fourier series )= ,k ,k (4) at an arbitrary set of targets. Nonuniform FFTs of the types discussed here are based, in essence, on combining some interpolation scheme with the standard FFT. 4ddly enough, it was a number of years after their use in applications before a rigorous analysis of such schemes was introduced by Dutt and 6o-hlin [7(. 8ubse2uent papers, such as [1, 3, 9,

10(, described ,ariants based on alternati,e interpolation:approximation approaches. Before discussing the algorithm itself, we would li-e to comment briefly on appli= cations that in,ol,e e,aluation of (3) or the Fourier integral ,s )= ,s (7) when the transform data ) is -nown at a scatter of points rather than on a regular Cartesian mesh. There is some confusion in the literature about the use of NUFFTs in this context (see 6emar- 1 below). There are three separate issues in,ol,ed: ac2uisi= tion of data ) in the Fourier domain, the choice of a 2uadrature scheme ,w and the a,ailability

of a fast algorithm for computing the discrete transform itself. They are often blended together when they should not be. Dutt and 6o-hlin [7( ap= pear to ha,e been the .rst to try to isolate one of these problems; they addressed the algorithmic 2uestion and showed that sums of the form (1) or (2) can be computed in log ) time with complete control of precision. )hile they concentrated on the one=dimensional case, higher dimensional ,ersions ha,e been considered by a ,a= riety of authors [3, 6, 14(. The .rst rigorous two=dimensional ,ersion can be found in a paper by 8train [17(, which uses

the NUFFT to sol,e a class of elliptic partial differential e2uations. Remark 1. Not all schemes for reconstructing Fourier integrals of the type (3) can be represented formally as a quadrature of the type (1) . Different schemes for interpolating to a uniform mesh do give rise to different reconstructed functions . However, if the decision has been made to use a quadrature approach such as (1) then the remaining task is entirely computational. The best algorithm is the one that evaluates the relevant sums as quickly and as accurately as possible. 2. The NUFFT. Let us .rst

consider the one=dimensional analogs of the summa= tion problems (1) and (2). )ith [0 (, the type=1 NUFFT is de.ned by the
Page 3
ACCELERATINGT0EN1N,NI21RM2AST21,RIERTRANS21RM 445 calculation of )= =0 ikx for ,..., (6) It is based on the following set of obser,ations: 1. E2uation (6) describes the exact Fourier coefficients of the function )= =0 (7) ,iewed as a periodic function on [0 (. Here, ) denotes the Dirac delta function. It is clearly not well=resol,ed by a uniform mesh in 2. Let ) denote the one=dimensional periodic heat -ernel on [0 (, gi,en by )= l If we de.ne ) to be

the con,olution )= )= dy, (8) then isa2 =periodic function and can be well=resol,ed by a uniform mesh in whose spacing is determined by (Figure 1). The Fourier coeffi= cients of , namely, )= ikx dx, can be computed with high accuracy using the standard FFT on an o,ersam= pled grid =0 (2 πm/M ik πm/M (9) where (2 πm/M )= =0 (2 πm/M (10) 3. 4nce the ,alues ) are -nown, an elementary calculation shows that )= (11) (This is a direct conse2uence of the con,olution theorem and the fact that the Fourier transform of is )= τe .)
Page 4
446

LESLIEGREENGARDANDJ,NE-Y,BLEE image image Fig. 1 In the version of the NUFFT described here, each delta function source at a point such as in (=) isreplacedbyaGaussian. Thissmearsthesourcestrengthtonearbyregulargrid points. Theregulargridmustbefineenoughtoresolvethesmearedfunction in (8) . Note that we include -periodic images of the sources in the definition of the heat kernel Theydecaysufficientlyrapidlythatallbutthenearestonescanbeignored. 6ecall that a type=2 transformation e,aluates a regular Fourier series at irregular target points. Thus, in one dimension, for [0 (, the

type=2 NUFFT is de.ned by the calculation of )= ikx (12) It is based on a closely related idea: 1. )e .rst decon,ol,e the Fourier coefficients, de.ning )by )= (13) and e,aluate the corresponding function ) on a uniform mesh with points on [0 ( using the FFT )= =0 ikx (14) In the preceding expression, we set )=0for k and ,iew )asan =periodic function )= ). 2. )e then compute the desired ,alues ) from )= )= dx =0 (2 πm/M πm/M (17) This is again a direct conse2uence of the con,olution theorem except that we decon,ol,e the effect of Gaussian smoothing in (13) before we

actually carry out the smoothing in (17)C
Page 5
ACCELERATINGT0EN1N,NI21RM2AST21,RIERTRANS21RM 443 Remark 2. There are a number of details that need to be fixed here, including the selection of , the convolution with a Gaussian in both p rocedures, and the length of the FFTs used. We will not repeat the analysis of [7( , since the relevant results can be summarized very simply: with =2 and =12 /M , Gaussian spreading of each source to the nearest 24 grid points yields about 12 digits of accuracy. With =6 /M , Gaussian spreading of each source to the nearest 12 grid points

yields about digits of accuracy. 3. Fast Gaussian Gridding. The dominant tas- in the NUFFT is the calcula= tion of (2 πm/M ) in (10) and ) in (17). Following standard practice, we will refer to these processes as gridding and the =point mesh as the oversampled mesh. In dimensions, gridding re2uires 12 exponential e,aluations for single precision accuracy and about 24 exponential e,aluations for double precision accuracy. In order to a,oid that cost using existing schemes, one can precompute all the necessary 2uantities, incurring a storage cost of 12 or 24 and a computational cost of 12

or 24 multiplications. This cost (in either storage or CDU time or both) becomes a signi.cant burden in two, three, and higher dimensions. It is sometimes called the curse of dimensionality; in the absence of a separable coordinate system, interpolation=type processes ha,e costs that grow exponentially with dimension. In the remainder of this paper, we donEt o,ercome the curse, but we show that (1 F exponential e,aluations or (1F storage is sufficient, followed by (12 F12 or (24 F24 multiplications, depending on the re2uired accuracy. The net reduction in CDU time is by a factor of 7G10.

By inspection of (10), it is e,ident that we only need ,alues of the function at e2uispaced points on the o,ersampled mesh. For this, we ha,e (2 πm/M )= =0 πm/M l This expression for loo-s much more expensi,e than it actually is. 8ince the Gaussian sources are sharply pea-ed (in a manner dependent on ), each source of strength located at is nonnegligible only at nearby grid points. Hs mentioned abo,e, we only need to compute its effect at the nearest 12 or 24 points (Figure 1) to achie,e either 6= or 12=digit accuracy. Thus, for the purpose of computation, we change our point

of ,iew from the recei,ing point (2 πm/M ) to the source point and consider one Gaussian source at a time. Hn elementary calculation shows that πm/M π/M πm/M / (16) But this means that we can compute and store two exponentials, and π/M for each source point. The third exponential, πm/M / , is independent of The gridding algorithm is straightforward: Let =2 πm/M denote the nearest regular grid point that is less than or e2ual to on the o,ersampled grid. Beginning at , let sp denote the number of grid points to which the spreading will be accounted for in each

direction. Compute the two exponentials: π/M The spreading contribution to (2 /M ) for sp sp is ), where sp = 6 for single precision, sp = 12 for double precision, and )= πm /M /
Page 6
444 LESLIEGREENGARDANDJ,NE-Y,BLEE Careful organi1ation of the loop shows that, for each source point, two exponential e,aluations are re2uired, followed by two multiplications at each of 2 sp regular mesh points. The algorithm for (17) is similar. H nice feature of Gaussian spreading is that the heat -ernel is built as a tensor product: πm/M πn/M π/M π/M πm/M /

πn/M / Thus, one can carry out spreading one dimension at a time. There follows an informal description of the implementation of the type=1 fast gridding algorithm in two dimensions, ,k )= =0 ,k for ,k and [0 [0 (. FastGriddingAlgorithmofType1inTwoDimensions. StepI:Initialization 1. 8et the o,ersampling ratio /M , the spreading parameter sp , and the Gaussian -ernel parameter according to the desired precision 2. Drecompute )= πl/M / for 0 sp and )= τk for | StepC:ConvolutionforEachSourcePoint ,y 1. Find the nearest grid point ( , )= ,m ) with 2. Compute (( +( /M /M and )= )=

for sp ,l sp 3. Con,ol,e the Gaussian spreading function with as follows: for sp F1 ,M sp for sp F1 ,M sp Hdd )to ,m ). StepD:FFTandDeconvolution 1. Compute two=dimensional FFT of ,m ) to obtain ,k ). 2. 8et ,k )= ,k ) for ,k The total cost is that of the o,ersampled two=dimensional FFT (8tep D), ( F1) ex= ponential e,aluations (8tep C2), ( sp ) multiplications (8tep C2), and (2 sp mul= tiplications (8tep C3) per source point. 4. Numerical E(am)les. The NUFFT of types 1 and 2 ha,e been implemented (in Fortran) with fast gridding in one and two dimensions. The following is a selection of

examples illustrating their performance. E(am)le 1 (Verification of Accuracy). In order to chec- the accuracy of the gridding algorithm numerically, we compare the results of the type=1 and =2 transforms
Page 7
ACCELERATINGT0EN1N,NI21RM2AST21,RIERTRANS21RM 445 10 15 10 −15 10 10 10 10 Type 1: 1024 points in 1D Spreading distance Computed Error 10 15 10 15 10 10 10 10 Type 1&2: 64*64 grid Spreading distance 10 15 10 15 10 10 10 10 Type 1&2: 16384 points in 2D Spreading distance Fig. 2 Gridding error in the NUFFT. In each figure, the -axis indicates the spreading

distance definedbyparameter sp andthe -axisindicatesthecomputederror. Theleftmostfigure correspondstoaone-dimensionalexamplewith 1824 points. Themiddleexampleusesauni- formsetofdatapointsintwodimensionsandtherightmostfigureusesarandomdistribution ofdatapointsintwodimensions. (ash-dottedlinesshowtheerrorforanoversampledgrid with I1 ,solidlinesfor I2 ,anddashedlinesfor I3 Table 1 sp I1 I2 I2 I3 I3 I4 9.8G73 1.9G73 8.9G74 9.3G74 3.8G74 3.1G74 8.1G79 3.9G76 =.2G7= 2.8G7= 1.9G7= 1.8G7= =.2G7= 6.9G79 6.2G718 1.9G718 9.8G711 3.8G711 12 6.9G79 1.2G711 9.9G713 8.8G714 2.3G714 9.2G719

(6), (12) with uniformly distributed random data points using direct summation and the NUFFT. Figure 2 shows the errors in the =norm. In the fast gridding algorithm, the accuracy of the type=1 transformation (6) is controlled by three parameters: the Gaussian -ernel parameter , the o,ersampling ratio /M , and the spreading distance sp . )e set according to the formula 7) sp Hs indicated earlier, we refer to the paper [7( for a detailed analysis. Here, we present some numerical experiments ,erifying the precision estimates deri,ed there (see Ta= ble 1). E(am)le 2 (Fast Gridding Com)ared to

Gridding). The nai,e Gaussian grid= ding algorithm and the fast gridding method for (1), (2), (6), and (12) ha,e been tested in ,arious computing en,ironments, including 8olaris 2.7 with gcc=2.97 on a 470MH1 Ultra=60 with 768MB 6HM, cygwin=7.1 with gcc=3.2 on a 1GH1 Dentium=3 with 384MB 6HM, and cygwin=7.1 with gcc=3.2 on a 2.4GH1 Dentium=4 with 1GB 6HM. Figure 3 summari1es the computational costs for =2 sp = 6 yielding six digits of accuracy. Both algorithms use the standard FFT on the o,ersampled mesh, and the time for this step is indicated in Figure 3 by dotted lines. The actual

computation time depends on a number of factors, including compiler options, type of CDU, performance
Page 8
450 LESLIEGREENGARDANDJ,NE-Y,BLEE 0 50,000 100,000 0.2 0.4 0.6 0.8 1.2 1.4 1D: Ultra 60(450Mhz) Number of Data Points Time (Sec) 0 50,000 100,000 0.2 0.4 0.6 0.8 1D: Pentium III(1GHz) Number of Data Points Time (Sec) 0 50,000 100,000 0.1 0.2 0.3 0.4 0.5 0.6 0.7 1D: Pentium IV(2.4Ghz) Number of Data Points Time (Sec) 0 50,000 100,000 10 2D: Ultra 60(450Mhz) Number of Data Points Time (Sec) 0 50,000 100,000 10 2D: Pentium III(1GHz) Number of Data Points Time (Sec) 0 50,000

100,000 2D: Pentium IV(2.4Ghz) Number of Data Points Time (Sec) Fig. 3 Computingtime. (ottedlinesrepresentthecostforFFT,dashedlinesforfastgridding,and dash-dottedlinesfornaivegridding. )eavydots( )areusedfortype- andplusmarks( fortype- transformation. of the math coprocessor, cache si1e, etc. Note that the type=1 and type=2 transforms are ,ery similar in terms of floating point operations; the differences in CDU time are due mainly to memory caching issues. In any case, the speed=up of the fast gridding algorithm is signi.cant in two dimensions and would be e,en more signi.cant in

the three=dimensional case. )e used an optimi1ed ,ersion of the FFT and compiled all codes with the 42 optimi1ation flag, but we did not carry out a complete optimi1ation of the fast gridding algorithm itself. More detailed (but less portable) implementation wor- could probably yield a further factor of 7G10 in performance. )e ha,e not carried out such .ne=tuning. E(am)le 3 (Com)arison with Direct Method). In this example, we compare the computational performance of our fast gridding algorithm with direct summation and the standard FFT. The direct code was implemented and compiled with

the same options: the gcc=2.97 compiler with =42 optimi1ation on a 470MH1 8parc Ultra=60. Figure 4 shows our results. )e set =2 , and set the number of data points in one dimension and in two dimensions. Using the data of Figure 4, we can summari1e the performance of the method in one and two dimensions as follows. First, direct summation re2uires about 0.022 sec) in one dimension and 0.087 sec) in two dimensions. In one dimension, the o,ersampled FFT re2uires about 2.4 sec (using a linear .t of the data o,er the range of tested). Note that this is already twice as expensi,e as an =point


Page 9
ACCELERATINGT0EN1N,NI21RM2AST21,RIERTRANS21RM 456 10 10 10 10 10 10 1D Transformation Cost Number of Data Points Time (Sec) 10 10 10 10 10 10 10 10 2D Transformation Cost Number of Data Points Time (Sec) Fig. 4 C,Ure-uirementsofthedirectcode(solidlineswithdots),thefastgriddingcode(foursolid lineswithtolerance I18 18 18 18 12 frombottomtotop),andthestandardFFT fortheoversampledmeshwith I2 (dottedlines). Table 2 1D I18 I18 I18 I18 12 Time( sec) 3.==N 4.48N 4.=8N 9.44N Krea*even( 1=8 288 228 298 2D I18 I18 I18 I18 12 Time( sec) 8.49N 28.6=N 36.39N 99.38N Krea*even( 18L18 19L19

28L28 26L26 FFT. In two dimensions, the o,ersampled FFT re2uires about 2.7 sec. This is four times as expensi,e as an FFT. Table 2 shows how the computational cost for gridding grows with precision. It also shows the brea-=e,en point with respect to the direct algorithm. In summary, the NUFFT is about 4 times more expensi,e in one dimension than a traditional =point FFT for single precision accuracy. It is about 27 times slower in the current implementation than the traditional two=dimensional FFT—a factor of 21 from gridding and a factor of 4 from the o,ersampled FFT. E(am)le 4 (M.IImage

.econstruction). 4ne of the important applications of the nonuniform FFT is to magnetic resonance imaging (M6I) [6, 8, 11, 12, 13, 14(. The M6I hardware is able to ac2uire the Fourier transform of a particular tissue property at selected points in the fre2uency domain. In most clinical systems, the de,ice is designed to ac2uire data on a uniform Cartesian mesh, from which a standard FFT can be used for image reconstruction. For a ,ariety of technical reasons, howe,er, nonuniform data sampling techni2ues are much better suited for fast data ac2uisition, motion correction, and functional M6I

[4(. In this example, we create simulated M6I
Page 10
452 LESLIEGREENGARDANDJ,NE-Y,BLEE data by using a type=2 transformation in two dimensions: ,s )= ,j ,j ,s (17) followed by a type=1 transformation to reconstruct the image, ,j )= =0 ,j ,s (18) If the function ,s ) were -nown e,erywhere, then the exact reconstruction would ob,iously be the Fourier integral ,j )= r, ,j cos ,r sin rdrdθ, (19) written in polar coordinates. It is probably worth repeating a point made in the introduction: once the decision has been made to use (18) for reconstruction, one still has a number of

degrees of freedom to wor- with. Engineering considerations determine the selection of points ,s , which will certainly affect the image 2uality. 4ne must also select 2uadrature weights so that, in the transform (18), ,s ). There are a number of interesting optimi1ation 2uestions that arise here, which will be addressed in subse2uent wor-. Here, we will simply use a radial grid and truncate the integral (19) at introducing a LringingM artifact that can be seen in the reconstruction. For a .xed NUFFT tolerance , the result computed agrees with the exact sum (18) to within that error. More

precisely, we let ,s )= (cos( sin( )) ,r πj , πi (20) for ,0 j ,0 i< , so that =2 . Thus, the 2uadrature weight for the point indexed by is =( jπ/M (2 π/ π/M ). Figure 7 shows the image reconstructed in this manner from the well=-nown 8heppGLogan phantom sampled on a 276 276 grid. )e set the Fourier transform tolerance to =10 5. Conclusions. The nonuniform FFT (NUFFT) is an important, and relati,ely recent, algorithm. There is, howe,er, some confusion in the literature about its use. It is simply a fast algorithm for computing discrete sums of a certain type. It is

completely independent of considerations ha,ing to do with ac2uisition of data in the Fourier domain or the choice of a 2uadrature scheme in computing Fourier integrals. The use of Gaussian spreading inside the algorithm has no theoretical ad,antage o,er any other properly applied Lspreading function.M It does, howe,er, allow a particularly simple and fast implementation, as described abo,e. )e belie,e it pro,ides the .rst reasonably efficient three=dimensional scheme. )e do not mean to suggest, howe,er, that all schemes for Fourier reconstruction must be based on a 2uadrature approach

and the calculation of sums li-e (18) or (1). H more general linear reconstruction algorithm could ta-e the form ,k )= =0 j,k ,k ,k (21)
Page 11
ACCELERATINGT0EN1N,NI21RM2AST21,RIERTRANS21RM 453 10 20 30 40 50 y=0.250 10 20 30 40 50 y= 0.609 10 20 30 40 50 Fig. 5 .econstructed image from Fourier data sampled on a radial (polar coordinate) grid. The curvesattherightshowthe reconstructedfunctionalongthelines I8 29 and 689 Here the weight depends on both the LsourceM ( ) and the LtargetM ( ,k ). The NUFFT does not apply to this calculation. RE2ERENCES [12 C.AndersonandM.D.Dahleh

.apidcomputationofthediscreteFouriertransform ,SIAM J.Sci.Comput.,1=(1996),pp.9131919. [22 S. Bagchi and S. Mitra TheNonuniform(iscreteFourierTransformandItsApplications inSignal, rocessing ,KluwerAcademic,Koston,1999. [32 G. Beylkin On the fast Fourier transform of functions with singularities , Appl. Comput. MarmonicAnal.,2(1999),pp.3631383. [42 M.Bourgeois,F.Wajer,D.Ormondt,andD.Graveron-Demilly .econstructionofM.I imagesfromnon-uniformsamplingandapplicationtoIntrascanmotioncorrectioninfunc- tionalM.I ,inModernSamplingTheoryC MathematicsandApplications,J.J.Kenedetto

and;.Ferreira,eds.,Appl.Numer.Marmon.Anal.,Kir*hN auser,Koston,2881,pp.3431363. [92 A. Dutt and -. Rokhlin Fast Fourier transforms for none-uispaced data , SIAM J. Sci. Comput.,14(1993),pp.136811393. [62 ..A.FesslerandB.P.Sutton NonuniformfastFouriertransformsusingmin-maxinter- polation ,IGGGTrans.Signal;rocess.,91(2883),pp.96819=4. [=2 L.GreengardandP.Lin Spectralapproximationofthefree-spaceheatkernel ,Appl.Comput. MarmonicAnal.,9(2888),pp.8319=. [82 ..I..ackson,C.H.Meyer,D.G.1ishimura,andA.Macovski Selectionofaconvolution function for Fourier inversion using gridding , IGGG Trans. Med.

Imag., 18 (1991), pp. 4=314=8. [92 2. H. Liu and 1. 1guyen An accurate algorithm for nonuniform fast Fourier transforms (NUFFT) ,IGGGMicrowaveGuidedWaveOett.,8(1998),pp.18128. [182 1.1guyenand2.H.Liu TheregularFouriermatricesandnonuniformfastFouriertrans- forms ,SIAMJ.Sci.Comput.,21(1999),pp.2831293. [112 .. D. O’Sullivan AfastsincfunctiongriddingalgorithmforFourierinversionincomputer tomography ,IGGGTrans.Med.Imag.,MI74(1989),pp.288128=. [122 G.B.Pike MultidimensionalFouriertransformationinmagneticresonanceimaging ,inThe

FourierTransformationinKiomedicalGngineering,T.;etersandJ.Williams,eds.,Appl. Numer.Marmon.Anal.,Kir*hN auser,Koston,1998,pp.891128.
Page 12
454 LESLIEGREENGARDANDJ,NE-Y,BLEE [132 D. Potts, G. Steidl, and M. Tasche Fast Fourier transforms for none-uispaced data: A tutorial ,inModernSamplingTheoryC MathematicsandApplications,J.J.Kenedettoand ;.Ferreira,eds.,Appl.Numer.Marmon.Anal.,Kir*hN auser,Koston,2881,pp.24912=4. [142 G. E. Sarty, R. Bennett, and R. W. Cox (irect reconstruction of non-Cartesian -space datausinganonuniformfastFouriertransform ,Magn.)eson.Med.,49(2881),pp.9881 919.

[192 .. Strain Fastpotentialtheory IIC Layerpotentialsanddiscretesums ,J.Comput.;hys.,99 (1992),pp.29112=8. [162 A. R. Thompson and R. 1. Bracewell InterpolationandFouriertransformationoffringe visibilities ,Astronom.J.,=9(19=4),pp.11124. [1=2 A.F.Ware FastapproximateFouriertransformsforirregularlyspaceddata ,SIAM)ev.,48 (1998),pp.8381896.