1
- from pynfft .nfft import NFFT
1
+ #from future import __division__
2
+ from nfft import nfft_adjoint
2
3
from .utils import Summations
3
4
import numpy as np
4
5
from math import floor
@@ -67,20 +68,24 @@ def direct_summations(t, y, w, freqs, nh):
67
68
68
69
69
70
70
- def fast_summations (t , y , w , freqs , nh , eps = 1E-5 ):
71
+ def fast_summations (t , y , w , freqs , nh , sigma = 2 , tol = 1E-7 , m = None ,
72
+ kernel = 'gaussian' , use_fft = True , truncated = True ):
71
73
"""
72
74
Computes C, S, YC, YS, CC, CS, SS using
73
- pyNFFT
75
+ nfft Python implementation by Jake Vanderplas
74
76
"""
77
+ nfft_kwargs = dict (sigma = sigma , tol = tol , m = m ,
78
+ kernel = kernel , use_fft = use_fft ,
79
+ truncated = truncated )
80
+
75
81
nf , df , dnf = inspect_freqs (freqs )
76
82
tmin = min (t )
77
83
78
84
# infer samples per peak
79
85
baseline = max (t ) - tmin
80
86
samples_per_peak = 1. / (baseline * df )
81
87
82
- eps = 1E-5
83
- a = 0.5 - eps
88
+ a = 0.5 - 1E-8
84
89
r = 2 * a / df
85
90
86
91
tshift = a * (2 * (t - tmin ) / r - 1 )
@@ -90,38 +95,23 @@ def fast_summations(t, y, w, freqs, nh, eps=1E-5):
90
95
# nf_nfft_w / 2 - 1 = 2H * (nf - 1 + dnf)
91
96
nf_nfft_u = 2 * ( nh * (nf + dnf - 1 ) + 1 )
92
97
nf_nfft_w = 2 * ( 2 * nh * (nf + dnf - 1 ) + 1 )
93
- n_w0 = int (floor (nf_nfft_w / 2 ))
94
- n_u0 = int (floor (nf_nfft_u / 2 ))
95
98
96
99
# transform y -> w_i * y_i - ybar
97
100
ybar = np .dot (w , y )
98
101
u = np .multiply (w , y - ybar )
99
102
100
- # plan NFFT's and precompute
101
- plan = NFFT (nf_nfft_w , len (tshift ))
102
- plan .x = tshift
103
- plan .precompute ()
104
-
105
- plan2 = NFFT (nf_nfft_u , len (tshift ))
106
- plan2 .x = tshift
107
- plan2 .precompute ()
108
-
109
- # NFFT(weights)
110
- plan .f = w
111
-
112
- f_hat_w = plan .adjoint ()[n_w0 :]
113
-
114
- # NFFT(y - ybar)
115
- plan2 .f = u
116
- f_hat_u = plan2 .adjoint ()[n_u0 :]
103
+
104
+ n_w0 = int (floor (nf_nfft_w / 2 ))
105
+ n_u0 = int (floor (nf_nfft_u / 2 ))
106
+ f_hat_u = nfft_adjoint (tshift , u , nf_nfft_u , ** nfft_kwargs )[n_u0 :]
107
+ f_hat_w = nfft_adjoint (tshift , w , nf_nfft_w , ** nfft_kwargs )[n_w0 :]
117
108
118
109
# now correct for phase shift induced by transforming t -> (-1/2, 1/2)
119
110
beta = - a * (2 * tmin / r + 1 )
120
111
I = 0. + 1j
121
112
twiddles = np .exp (- I * 2 * np .pi * np .arange (0 , n_w0 ) * beta )
122
113
f_hat_u *= twiddles [:len (f_hat_u )]
123
114
f_hat_w *= twiddles [:len (f_hat_w )]
124
-
125
115
all_computed_sums = []
126
116
127
117
# Now compute the summation values at each frequency
0 commit comments