;+ ; NAME: matrix_dft ; ; PURPOSE: ; Matrix Discrete Fourier Transform, following Soummer et al. 2007. ; This is not as fast as an FFT, but allows you to arbitrarily choose the ; sampling and range covered in the Fourier domain. ; ; NOTES: ; This is an IDL translation of the Python code in SFT.py. ; ; Compute an adjustible-center matrix fourier transform. ; ; For an output array with ODD size n, ; the PSF center will be at the center of pixel (n-1)/2 ; ; For an output array with EVEN size n, ; the PSF center will be in the corner between pixel (n/2-1,n/2-1) and (n/2,n/2) ; ; Those coordinates all assume Python/IDL style pixel coordinates running from ; (0,0) up to (n-1, n-1). ; ; ; This version supports rectangular, non-square arrays, ; in which case nlamD and npix should be 2-element ; tuples or lists, using the usual Pythonic order (Y,X) ; ; INPUTS: ; ; ; KEYWORDS: ; Parameters ; ---------- ; pupil : array ; pupil array (n by n) ; nlamD : float or tuple ; size of focal plane array, in units of lam/D ; (corresponds to 'm' in Soummer et al. 2007 4.2) ; npix : float or tuple ; number of pixels per side side of focal plane array ; (corresponds to 'N_B' in Soummer et al. 2007 4.2) ; offset: tuple ; an offset in pixels relative to the above ; ; OUTPUTS: ; ; HISTORY: ; Began 2011-08-16 18:47:19 by Marshall Perrin based on SFT.py ; Benchmarked against SFT & agreement achieved. ;- PRO test_matrix_dft ; Test matrix DFT, including non-square arrays, in both the ; forward and inverse direction. ; ; This is an IDL equivalent to test_SFT_rect in SFT.py; they should give ; identical results. npupil=156 indices, intarr(npupil, npupil), x, y x -= npupil*0.5 y -= npupil*0.5 r = sqrt(x^2+y^2) pupil = r lt npupil/2 pupil[0:60-1, 0:60-1] = 0 pupil[*, 0:10-1] = 0 pupil /= sqrt(total(pupil)) nlamD=[20,10] sampling=10 npix=nlamD*sampling asf = matrix_dft(pupil, nlamD, npix) !p.multi=[0,4,1] !p.charsize=2 loadct, 33 ; = numpy Jet imdisp, bytscl(pupil,0,8e-3*1.5),/axis, title='Pupil',/noscale imdisp, bytscl(alog10(float(asf)), -8,1),/axis, title='ASF',/noscale psf = asf * conj(asf) imdisp, bytscl(alog10(float(psf)), -8,1),/axis, title='PSF',/noscale pupil2 = matrix_dft(asf, nlamD, npupil,/inverse) pupil2r = float(pupil2 * conj(pupil2)) imdisp, bytscl(pupil2r,0,8e-5*1.5), /axis, title='Pupil (inverse)', /noscale end ;-------- function matrix_dft, pupil, nlamD, npix, offset=offset, inverse=inverse if ~(keyword_set(inverse)) then inverse=0 if ~(keyword_set(offset)) then offset=[0,0] sz = size(pupil) npupY = sz[2] npupX = sz[1] if n_elements(npix) ge 2 then begin npixX = npix[0] npixY = npix[1] endif else begin npixX = npix npixY = npix endelse if n_elements(nlamD) ge 2 then begin nlamDX = nlamD[0] nlamDY = nlamD[1] endif else begin nlamDX = nlamD nlamDY = nlamD endelse if keyword_set(inverse) then begin dX = nlamDX / float(npupX) dY = nlamDY / float(npupY) dU = 1.0/float(npixY) dV = 1.0/float(npixX) ;pupil = transpose(pupil) ;this is necessary because the matrix code assumes Pythonic axes order and ; now we're in IDL. Empirically this works ; here... endif else begin dU = nlamDX / float(npixX) dV = nlamDX / float(npixX) dX = 1.0/float(npupX) dY = 1.0/float(npupY) endelse Xs = (findgen(npupX) - float(npupX)/2.0 - offset[1] + 0.5) * dX Ys = (findgen(npupY) - float(npupY)/2.0 - offset[0] + 0.5) * dY Us = (findgen(npixX) - float(npixX)/2.0 - offset[1] + 0.5) * dU Vs = (findgen(npixY) - float(npixY)/2.0 - offset[0] + 0.5) * dV XU = Xs ## Us ; np.outer(Xs, Us) YV = Ys ## Vs ; np.outer(Ys, Vs) expXU = exp(-2.0 * !pi * complex(0,1) * XU) expYV = exp(-2.0 * !pi * complex(0,1) * YV) expYV = TRANSPOSE(expYV) ;.T.copy() t1 = expYV ## pupil ; np.dot(expYV, pupil) t2 = t1 ## expXU ; np.dot(t1, expXU) if keyword_set(inverse) then begin ;#print expYV.shape, pupil.shape, expXU.shape ;t2 = t2[::-1, ::-1] t2 = reverse(reverse(t2,1),2) endif else begin ;t2 = transpose(t2) ; this is necessary because the above assumes Pythonic axes order and ; now we're in IDL... empirically this works to ; transpose here. endelse norm_coeff = sqrt( float( nlamDY* nlamDX) / (npupY*npupX*npixY*npixX)) return, norm_coeff * t2 end