calculation was for nyquist freq not fully correct.

See Notes on FFT-based differentiation
Steven G. Johnson, MIT Applied Mathematics
Created April, 2011, updated May 4, 2011:
This commit is contained in:
Martin Diehl 2015-09-05 11:51:36 +00:00
parent ba418a5c97
commit cf90b764f8
2 changed files with 28 additions and 8 deletions

View File

@ -29,13 +29,23 @@ def curlFFT(geomdim,field):
TWOPIIMG = (0.0+2.0j*math.pi)
for i in xrange(grid[0]):
k_s[0] = i
if(i > grid[0]/2 ): k_s[0] = k_s[0] - grid[0]
if(grid[0]%2==0 and i == grid[0]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
k_s[0]=0
elif (i > grid[0]//2):
k_s[0] = k_s[0] - grid[0]
for j in xrange(grid[1]):
k_s[1] = j
if(j > grid[1]/2 ): k_s[1] = k_s[1] - grid[1]
for k in xrange(grid[2]/2+1):
if(grid[1]%2==0 and j == grid[1]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
k_s[1]=0
elif (j > grid[1]//2):
k_s[1] = k_s[1] - grid[1]
for k in xrange(grid[2]//2+1):
k_s[2] = k
if(k > grid[2]/2 ): k_s[2] = k_s[2] - grid[2]
if(grid[2]%2==0 and k == grid[2]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
k_s[2]=0
xi = np.array([k_s[2]/geomdim[2]+0.0j,k_s[1]/geomdim[1]+0.j,k_s[0]/geomdim[0]+0.j],'c16')
if dataType == 'tensor':
for l in xrange(3):

View File

@ -24,13 +24,23 @@ def divFFT(geomdim,field):
TWOPIIMG = (0.0+2.0j*math.pi)
for i in xrange(grid[0]):
k_s[0] = i
if(i > grid[0]/2 ): k_s[0] = k_s[0] - grid[0]
if(grid[0]%2==0 and i == grid[0]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
k_s[0]=0
elif (i > grid[0]//2):
k_s[0] = k_s[0] - grid[0]
for j in xrange(grid[1]):
k_s[1] = j
if(j > grid[1]/2 ): k_s[1] = k_s[1] - grid[1]
for k in xrange(grid[2]/2+1):
if(grid[1]%2==0 and j == grid[1]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
k_s[1]=0
elif (j > grid[1]//2):
k_s[1] = k_s[1] - grid[1]
for k in xrange(grid[2]//2+1):
k_s[2] = k
if(k > grid[2]/2 ): k_s[2] = k_s[2] - grid[2]
if(grid[2]%2==0 and k == grid[2]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
k_s[2]=0
xi=np.array([k_s[2]/geomdim[2]+0.0j,k_s[1]/geomdim[1]+0.j,k_s[0]/geomdim[0]+0.j],'c16')
if n == 9: # tensor, 3x3 -> 3
for l in xrange(3):