Radially averaged Fourier Spectrum -- 2d

It is often necessary, when dealing with multidimensional data to compute the power spectrum as a function of |K|, the magnitude of the wavenumber vector. This is accomplished by computing the FFT and averaging power over circles of diameter |K|. Though fairly straightforward in theory, writing such a routine is often cumbersome because of details related to the arrangement of Fourier coefficients and the discrete representation of the circular average. Currently, such a routine is not part of IDL's standard math library.

Using the DIST and SHIFT functions together with FFT, this can be programmed extremely concisely in IDL. The basic code is given below.

IDL> array = dist(200,200)
IDL> pspec = powerspectrum(array)
IDL> plot, pspec, /ylog

;
; Function PowerSpectrum accepts a square 2D array as input
; and returns a 1D array that includes the Fourier power spectrum
; as a function of the waveVector magnitude, |K| = sqrt(kX*kX + kY*kY)
;
function PowerSpectrum, inputArray
sZ = SIZE(inputArray)
if (sZ[0] ne 2 or (sZ[0] eq 2 and sZ[1] ne sZ[2])) then begin
   void = DIALOG_MESSAGE(/Error, ['Error in Power_Spectrum.',$
                         'Input needs to be a square 2D array.',$
                         'Calling sequence:',$
                         'outputSpectrum = PowerSpectrum(inputArray)'])                      
   RETURN, inputArray
endif
nX = sZ[1] & nWaveNumber = nX / 2
waveNumberDist = SHIFT(DIST(nX), nWaveNumber, nWaveNumber)
fftArray = SHIFT(FFT(inputArray, -1), nWaveNumber, nWaveNumber)
fftArray = FLOAT(fftArray*CONJ(fftArray))
outputSpectrum = FLTARR(nWaveNumber)
for waveNumber = 0, nWaveNumber-1 do outputSpectrum[waveNumber] = $
    TOTAL(fftArray[WHERE((waveNumberDist ge (waveNumber-0.5)) and $
                         (waveNumberDist lt (waveNumber+0.5)))])
RETURN, outputSpectrum
end