(* ::Package:: *) (* :Title: nfft *) (* :Author: Sseziwa Mukasa *) (* :Summary: nfft: nonequispaced Fourier transform This package exports seven functions: nfft, ndft, nfftTranspose, ndftTranspose, precompute\[Psi], precomputec and precomputen. The functions nfft and ndft approximate or calculate respectively the nonequispaced Fourier transform as described in D. Potts, G. Steidl, and M. Tasche, "Fast Fourier transforms for nonequispaced data: A tutorial" in Modern Sampling Theory: Mathematics and Applications, J.J. Benedetto and P. Ferreira (Eds.), Chapter 12, pages 249-274. The functions nfftTranspose and ndftTranspose approximate or compute the transposed nonequispaced Fourier transform described therein. precompute\[Psi], precomputec and precompute n compute data structures used by nfft and nfftTransposed. The structures are functions of some of the parameters of nfft and nfftTransposed and can be stored for reuse in future calls. See also: http://www.math.mu-luebeck.de/potts/nfft/. Comments, bug reports, etc. can be forwarded to Sseziwa Mukasa, smukasa@earthlink.net *) (* :Package Version: .10 *) (* :Context: nffft` *) (* :Copyright: Copyright 2005 by Sseziwa Mukasa The author, makes no representations, express or implied, with respond to this documentation, of the software it describes and contains, including without limitations, any implied warranties of mechantability or fitness for a particular purpose, all of which are expressly disclaimed. The author shall in no event be liable for any indirect, incidental, or consequential damages. *) (* :History: Version 0.10 by Sseziwa Mukasa, 2005 *) (* :Keywords: nonequispaced, Fourier, transform, NFFT *) BeginPackage["nfft`",{"GUIKit`"}] approximationWidth::usage = "Option to nfft package functions. Half the size of the approximation \ window - 1, equivalent to the parameter m described in the tutorial at \ http://www.math.mu-luebeck.de/potts/nfft/." oversamplingFactor::usage = "Option to nfft package functions. Amount the data is oversampled to form \ the regular grid, equivalent to the parameter \[Alpha] described in the \ tutorial at http://www.math.mu-luebeck.de/potts/nfft/." compressed::usage = "Option to nfft package functions. A boolean value indicating the form of the \ frequency list. False indicates the frequency list is in the form \ {\!\(\*SubscriptBox[\"\[Nu]\", \"11\"]\),\!\(\*SubscriptBox[\"\[Nu]\", \"12\"]\),...,\!\(\*SubscriptBox[\"\[Nu]\", RowBox[{\"1\", \"d\"}]]\)},{\!\(\*SubscriptBox[\"\[Nu]\", \"21\"]\),\!\(\*SubscriptBox[\"\[Nu]\", \"22\"]\),...,\!\(\*SubscriptBox[\"\[Nu]\", RowBox[{\"2\", \"d\"}]]\)},...}, True indicates the list is of the form \ {{\!\(\*SubscriptBox[\"\[Nu]\", \"11\"]\),\!\(\*SubscriptBox[\"\[Nu]\", \"21\"]\),...},{\!\(\*SubscriptBox[\"\[Nu]\", \"12\"]\),\!\(\*SubscriptBox[\"\[Nu]\", \"22\"]\),...},...} and the complete set of frequencies is \ Outer[List,Sequence@@\[Nu]]." parallel::usage = "Option to nfft package functions. A boolean value indicating whether to \ perform the operation in parallel, default is False." showProgress::usage = "Option to nfft package functions. A boolean value indicating whether to \ display a progress indicator, default is False." Evaluate[ToExpression["precomputed\[Psi]"]]::usage = "Option to nfft package functions. Table of \ \!\(\[Psi] \((\[Nu]\_j - l\/n)\)\) and the upper and lower bounds of l as \ described in the tutorial at http://www.math.mu-luebeck.de/potts/nfft/." precomputedc::usage = "Option to nfft package functions. Table of \ \!\(\(c\_k\) \((\(\[CurlyPhi]\&~\))\)\) at the nodes \ \!\(k \[Element] I\_N\%d\) as described in the tutorial at http://www.math.mu-luebeck.de/potts/nfft/." windowFunction::usage = "Option to nfft package functions. The function \[Psi] as described in the tutorial at http://www.math.mu-luebeck.de/potts/nfft/." windowFunctionTransformed::usage = "Option to nfft package functions. The function \ \!\(\(c\_k\) \((\(\[CurlyPhi]\&~\))\)\) as described in the tutorial at http://www.math.mu-luebeck.de/potts/nfft/." nfft::usage = "nfft[f,\[Nu]] approximates the nonequispaced Fourier transform (NFFT) \ (see tutorial at http://www.math.mu-luebeck.de/potts/nfft/) of data f at \ frequencies \[Nu]. Options are: approximationWidth; half the size of \ the approximating window - 1 (default 6), oversamplingFactor the ratio of \ the size of the approximating grid to the dimensions of the input data \ (default 2), precomputed\[Psi]; the approximation window \ evaluated at points \!\(\[Nu] - l\/n\ \((l \[Element] I\_\(n, m\))\)\) \ as computed by precompute\[Psi], precomputedc; \ as computed by precomputec and compressed; the form of the frequency list. \ The default window function is a Kaiser-Bessel function indicated by a \ precomputed\[Psi] of None. The default transformed window function is a modified \ Bessel function indicated by a precomputedc of None." nfftTranspose::usage = "nfftTranspose[f,\[Nu],n] approximates the transposed Nonequispaced \ Fourier Transform (NFFT) \ (see tutorial at http://www.math.mu-luebeck.de/potts/nfft/) of data f at \ frequencies \[Nu]. n is a list of the desired size of the result in \ each dimension. Options are: approximationWidth; half the size of \ the approximating window - 1 (default 6), oversamplingFactor the ratio of \ the size of the approximating grid to the dimensions of the input data \ (default 2), precomputed\[Psi]; the approximation window \ evaluated at points \!\(\[Nu] - l\/n\ \((l \[Element] I\_\(n, m\))\)\) \ as computed by precomputed\[Psi], precomputedc; \ as computed by precomputec and compressed; the form of the frequency list. \ The default window function is a Kaiser-Bessel function indicated by a \ precomputed\[Psi] of None. The default transformed window function is a modified \ Bessel function indicated by a precomputedc of None." ndft::usage = "ndft[f,\[Nu]] computes the nonequispaced Fourier transform (NFFT) \ (see tutorial at http://www.math.mu-luebeck.de/potts/nfft/) of data f at \ frequencies \[Nu] using an \[ScriptCapitalO](m\!\(8\+\(i = 1\)\%d n\_i\)\ ) algorithm. d is the dimensionality of the f, m is the number of \ frequencies \[Nu], \!\(n\_i\) is the size of f in dimension i." ndftTranspose::usage = "ndft[f,\[Nu],n] computes the transposed nonequispaced Fourier \ transform (NFFT) (see tutorial \ at http://www.math.mu-luebeck.de/potts/nfft/) of data f at \ frequencies \[Nu] using an \[ScriptCapitalO](m\!\(8\+\(i = 1\)\%d n\_i\)\ ) algorithm. d is the dimensionality of the f, m is the number of \ frequencies \[Nu], \!\(n\_i\) is the size of f in dimension i." precomputen::usage = "precomputen[s_] returns the dimensions of the oversampled data set \ (Max[oversamplingFactor*n,2*approximatingWidth+2]) for a data set of \ dimensions s. Options are approximationWidth; half the size of the \ approximating window - 1(default 6) , oversamplingFactor; the ratio of n \ to the dimensions of the input data (default 2)" Evaluate[ToExpression["precompute\[Psi]"]]::usage = "precompute\[Psi][\[Nu],n] computes the approximation function \ \!\(\[Psi] \((\[Nu]\_j - l\/n)\)\) at the nodes \!\(\[Nu]\_j\) where \ \[Nu]*n-approximationWidth <= l <= \[Nu]*n+approximationWidth (see the \ tutorial at http://www.math.mu-luebeck.de/potts/nfft/ for details). \ \[Nu] is the list of nodes at which to compute the Fourier sum, and n the \ size of the oversampled data set (see precomputen). Options are: \ approximationWidth; half the size of the approximating window - 1 (default \ 6), oversamplingFactor; the ratio of n to the dimensions of the input data \ (default 2), windowFunction; the function used as an approximation \ window, and compressed; a boolean value indicating the form of the frequency \ list. The default window function is a Kaiser-Bessel function indicated by a \ windowFunction of None. The result is suitable for use as the precomputed\[Psi] \ option." precomputec::usage = "precomputec[n] computes the Fourier transformed window \ function \!\(\(c\_k\) \((\(\[CurlyPhi]\&~\))\)\) at the nodes \ \!\(k \[Element] I\_N\%d\) (see the tutorial at \ http://www.math.mu-luebeck.de/potts/nfft/ for details). n is the \ size of the oversampled data set (see precomputen). Options are: \ approximationWidth; half the size of the approximating window - 1 (default \ 6), oversamplingFactor; the ratio of n to the dimensions of the input data \ (default 2), and windowFunctionTransformed; the Fourier transform of the \ function used as an approximation window. The default transformed window \ function is a Modified Bessel Function indicated by a \ windowFunctionTransformed of None. The result is an array of with the \ number of rows equal to the length of n and the number of columns equal to \ the length of the data in that dimension." Options[precomputen] = {approximationWidth -> 6, oversamplingFactor -> 2} Options[precompute\[Psi]] = {approximationWidth -> 6, oversamplingFactor -> 2, windowFunction -> None, compressed -> False, parallel -> False} Options[precomputec] = {approximationWidth -> 6, oversamplingFactor -> 2, windowFunctionTransformed -> None, parallel -> False} Options[ndft] = showProgress -> False Options[nfft] = {approximationWidth -> 6, oversamplingFactor -> 2, precomputed\[Psi] -> None, precomputedc -> None, compressed -> False, parallel -> False, showProgress -> False} Options[ndftTranspose] = showProgress -> False Options[nfftTranspose]={approximationWidth->6,oversamplingFactor->2, precomputed\[Psi]->None, precomputedc->None, compressed -> False, showProgress -> False} Begin["`Private`"] kaiserBessel[v_,n_,m_,a_] := With[{b=Pi*(2-1/a)}, If[Abs[v] < m/n, Sinh[b*Sqrt[m^2-n^2*v^2]] / Sqrt[m^2-n^2*v^2], If[Abs[v]==m/n,b, Sin[b*Sqrt[n^2*v^2-m^2]] / Sqrt[n^2*v^2-m^2]]] / Pi] kaiserBesselTransformed[k_,n_,m_,a_] := With[{b=Pi*(2-1/a)}, If[-n*(1-1/(2*a)) <= k <= n*(1-1/(2*a)), BesselI[0,m*Sqrt[b^2 - (2*Pi*k/n)^2]], 0]] precomputen[s_,opts___?OptionsQ] := Module[{m,a}, {m,a}={nfft`approximationWidth, nfft`oversamplingFactor} /. {opts} /. Options[precomputen]; If[IntegerQ[m],m=Table[m,{Length[s]}]]; If[NumberQ[a],a=Table[a,{Length[s]}]]; Max/@Transpose[{a*s,2*m+2}]] precompute\[Psi][w_/;ArrayDepth[w]==1,n_,opts___?OptionQ] := Module[{m,a,p,f}, {m,a,p} = {nfft`approximationWidth, nfft`oversamplingFactor, nfft`windowFunction} /. {opts} /. Options[precompute\[Psi]]; If[IntegerQ[m], m = Table[m, {Length[w]}]]; If[NumberQ[a], a = Table[a, {Length[w]}]]; f = If[p === None, kaiserBessel, p]; Table[ Module[{u = Floor[n[[d]]*w[[d]]-m[[d]]], o = Ceiling[n[[d]]*w[[d]]+m[[d]]]}, If[o - u < 2*m[[d]]+1, If[w[[d]] < 0, u = u - 1, o = o + 1]]; {Table[f[w[[d]]-l/n[[d]],n[[d]],m[[d]],a[[d]]], {l,u,o}], u, o}], {d, Length[w]}]] toTree[{}] := {} toTree[{x_}] := {x,{},{}} toTree[{x_,y_}] := {x,{},{y}} toTree[lst_] := With[{p=Quotient[Length[lst],2]}, {lst[[p]],toTree[lst[[;;p-1]]],toTree[lst[[p+1;;]]]}]; DistributeDefinitions[toTree]; precompute\[Psi][w_,n_,opts___?OptionQ] := Module[{m,a,p,f,iscompressed,doparallel}, {m,a,p,iscompressed,doparallel} = {nfft`approximationWidth, nfft`oversamplingFactor, nfft`windowFunction, nfft`compressed, nfft`parallel} /. {opts} /. Options[precompute\[Psi]]; If[IntegerQ[m], m = Table[m, {If[iscompressed,Length[w],Length[w[[1]]]]}]]; If[NumberQ[a], a = Table[a, {If[iscompressed,Length[w],Length[w[[1]]]]}]]; f = If[p === None, kaiserBessel, p]; If[doparallel, DistributeDefinitions[w,n]; DistributeDefinitions["nfft`Private`"]; Parallelize@@#, ReleaseHold[#]]&[Hold[ MapThread[ Function[{wd,nd,md,ad}, toTree[ Module[{u = Floor[nd*#-md],o = Ceiling[nd*#+md]}, If[o-u < 2*md+1, If[# < 0, u = u - 1, o = o + 1]]; {#,Table[f[#-l/nd,nd,md,ad], {l,u,o}],u,o}]&/@wd]], {If[iscompressed,Sort/@w,Union/@Transpose[w]],n,m,a}]]]] DistributeDefinitions[precompute\[Psi]]; precomputec[n_,opts___?OptionQ] := Module[{m,a,phat,f,doparallel}, {m,a,phat,doparallel} = {nfft`approximationWidth, nfft`oversamplingFactor, nfft`windowFunctionTransformed, nfft`parallel} /. {opts} /. Options[precomputec]; If[IntegerQ[m], m = Table[m, {Length[n]}]]; If[NumberQ[a], a = Table[a, {Length[n]}]]; f = If[phat === None, kaiserBesselTransformed, phat]; If[doparallel, DistributeDefinitions[n]; DistributeDefinitions["nfft`Private`"]; Parallelize@@#, ReleaseHold[#]]&[ Hold[ Table[f[k,n[[d]],m[[d]],a[[d]]], {d,Length[n]}, {k,Quotient[-n[[d]]/a[[d]],2]+Mod[-n[[d]]/a[[d]],2], Quotient[n[[d]]/a[[d]],2]-1+Mod[n[[d]]/a[[d]],2]}]]]] ndft[f_,w_,opts___?OptionQ] := Module[{k,dims=Dimensions[f],d,lims,showprogress,prog=0}, showprogress = nfft`showProgress /.{opts} /. Options[ndft]; d=Length[dims]; lims=Table[{k[i],Quotient[-dims[[i]],2]+Mod[dims[[i]],2], If[EvenQ[dims[[i]]],dims[[i]]/2-1,Quotient[dims[[i]],2]]}, {i,d}]; Monitor[ (prog++; Sum[f[[Sequence@@Table[k[i]+Quotient[dims[[i]],2]+1,{i,d}]]]* Product[Exp[-2*Pi*I*k[i]*#[[i]]],{i,d}], Evaluate[Sequence@@lims]])&/@w, If[showprogress,ProgressIndicator[prog,{0,Length[w]}],"Calculating"]]] ndftTranspose[f_,w_,s_,opts___?OptionQ] := Module[{k,d=Length[s],lims,showprogress,prog=0}, showprogress = nfft`showProgress /.{opts} /. Options[ndftTranspose]; lims=Table[{k[i],Quotient[-s[[i]],2]+Mod[s[[i]],2], If[EvenQ[s[[i]]],s[[i]]/2-1,Quotient[s[[i]],2]]}, {i,d}]; Monitor[ Table[prog++;Sum[f[[j]]*Exp[2*Pi*I*Table[k[i],{i,d}].w[[j]]],{j,Length[w]}], Evaluate[Sequence@@lims]], If[showprogress, ProgressIndicator[prog,{0,Times@@(#[[3]]-#[[2]]+1&/@lims)}], "Calculating"]]] treeFind[tree_,item_] := If[First[First[tree]]===item,Rest[First[tree]], If[itemm, nfft`oversamplingFactor->a], MapThread[treeFind[#1,#2]&,{psi,w}]]; Total[ Flatten[ g[[Sequence@@ ((l=If[NonNegative[#[[1]]],#[[1]]+1,#[[1]]]; u=If[NonNegative[#[[2]]],#[[2]]+1,#[[2]]]; If[l*u<0, Join[Range[#[[3]]+l+1,#[[3]]],Range[u]], Span[l,u]])&/@ Join[\[Psi][[All,{2,3}]],Transpose[{n}],2])]]* Outer[Times,Sequence@@\[Psi][[All,1]]]]]] DistributeDefinitions[computeTerm]; nfft[f_,w_,opts___?OptionQ] := Module[{m,a,psi,c,iscompressed,doparallel,n,g,l,ret,showprogress,prog=0}, {m,a,psi,c,iscompressed,doparallel}={nfft`approximationWidth, nfft`oversamplingFactor, nfft`precomputed\[Psi], nfft`precomputedc, nfft`compressed, nfft`parallel} /. {opts} /. Options[nfft]; showprogress = nfft`showProgress /.{opts} /. Options[nfft]; If[IntegerQ[m],m=Table[m,{If[iscompressed,Length[w],Length[w[[1]]]]}]]; If[NumberQ[a],a=Table[a,{If[iscompressed,Length[w],Length[w[[1]]]]}]]; n = Max/@Transpose[{a*Dimensions[f],2*m+2}]; a = n/Dimensions[f]; (*Step 2: Compute g = FFT (ghat)*) g = Fourier[ RotateLeft[ PadLeft[ PadRight[ (*Step 1: Compute ghat=f/c*) Fold[Transpose[#1/#2,Prepend[Range[1,ArrayDepth[#]-1], ArrayDepth[#]]]&,f, If[c===None,precomputec[n,nfft`approximationWidth->m, nfft`oversamplingFactor->a],c]], Quotient[n,2]+Quotient[Dimensions[f],2]+Mod[n,2]],n], Quotient[n,2]], FourierParameters->{1,-1}]; (*Step 3: Compute result_j=Sum[g_l psi (v_j-l/n),{l,I_n,m (v_j)}]*) If[doparallel, DistributeDefinitions["nfft`Private`"]]; ret= If[!iscompressed, If[doparallel, SetSharedVariable[prog]; Monitor[ ParallelMap[(prog++;computeTerm[g,#,n,m,a,psi])&,w], If[showprogress, ProgressIndicator[prog,{0,Length[w]}],"Computing"]], Monitor[ (prog++;computeTerm[g,#,n,m,a,psi])&/@w, If[showprogress, ProgressIndicator[prog,{0,Length[w]}],"Computing"]]], Flatten[If[doparallel, SetSharedVariable[prog]; Monitor[ Parallelize[ Outer[(prog++; computeTerm[g,{##},n,m,a,psi])&, ##]]&[Sequence@@w], If[showprogress, ProgressIndicator[prog,{0,Times@@Dimensions[w]}], "Computing"]], Monitor[ Outer[(prog++;computeTerm[g,{##},n,m,a,psi])&, Sequence@@w], If[showprogress, ProgressIndicator[prog,{0,Times@@Dimensions[w]}], "Computing"]]], Length[w]-1]]; Remove[g];Return[ret]] computeTransposedTerm[f_,w_,n_,m_,a_,psi_] := With[{\[Psi] = If[psi===None, precompute\[Psi][w,n,nfft`approximationWidth->m, nfft`oversamplingFactor->a], MapThread[treeFind[#1,#2]&,{psi,w}]]}, f*SparseArray[Flatten[Outer[List,Sequence@@ Transpose[Mod[#,n]+1&/@ Transpose[Range@@@\[Psi][[All,{2,3}]]]]],1]-> Flatten[Outer[Times,Sequence@@\[Psi][[All,1]]]],n]] nfftTranspose[f_,w_,s_,opts___?OptionQ] := Module[{m,a,psi,c,iscompressed,n,showprogress,prog=0}, {m,a,psi,c,iscompressed}={nfft`approximationWidth, nfft`oversamplingFactor, nfft`precomputed\[Psi], nfft`precomputedc, nfft`compressed} /. {opts} /. Options[nfftTranspose]; showprogress = nfft`showProgress /.{opts} /. Options[nfftTranspose]; If[IntegerQ[m],m=Table[m,{If[iscompressed,Length[w],Length[w[[1]]]]}]]; If[NumberQ[a],a=Table[a,{If[iscompressed,Length[w],Length[w[[1]]]]}]]; n=Max/@Transpose[{a*s,2*m+2}]; a=n/s; (*Step 3: Compute result = ghat/c*) Fold[ Transpose[#1/#2,Prepend[Range[1,ArrayDepth[#]-1],ArrayDepth[#]]]&, (*Step 2: Compute ghat = IFFT (g)*) RotateRight[InverseFourier[ (*Step 1: Compute g=Sum[f_j psihat (v_j-l/n),{j,I_n,m}]*) N[Normal[ Monitor[ Fold[(prog++; #1+computeTransposedTerm[#2[[1]],#2[[2]],n,m,a,psi])&, 0, Transpose[{f,If[iscompressed, Flatten[Outer[List,Sequence@@w], Length[w]-1],w]}]], If[showprogress, ProgressIndicator[prog, If[iscompressed, Times@@Dimensions[w], Length[w]]], "Computing"]]]], FourierParameters->{-1,-1}],Quotient[s,2]][[Sequence@@Range[s]]], If[c===None,precomputec[n,nfft`approximateWidth->m, nfft`oversamplingFactor->a],c]]] End[] (*'Private'*) EndPackage[]