The following IDL code is provided for estimating 
noise power spectra (NPS) without any guarantee or claim.  
The author (KMH) and LANL do not accept any responsibility 
for its correctness.  It may be freely distributed.   
In this file, the individual IDL procedures/functions 
are separated by a string of semi-colons.  
You need to break out each procedure and place it in 
a separate file in the appropriate IDL subdirectory. 
To calculate the radial dependence of the NPS for an 
input (square) noise image, fn, which has been previously 
flattened by subtracting a slowly varying fit to the 
original noise image, you should enter:spec_calc, fn, 
dpixwhere dpix is pixel spacing (in mm if you want the 
NPS in image units**2 times mm**2).  You may need to adjust 
the value of NFRQ set by the statementnfrq = 7   ; 
number of freqsto get more or fewer frequency points. 
Reference for algorithm:
"A simplified method of estimating noise power spectra," 
K. M. Hanson, Proc. SPIE, vol. 3336,  in Physics of Medical Imaging}, 
J. T. Dobbins III and J. M. Boone, eds., 1998.
=================================================IDL CODE FOR NPS ESTIMATION           LA-CC-98-15;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro nps_calc, fn, dpixin, sp, freq, sperr

; improved version - runs faster thanks to Todd Kauppila

;SCOTT WATSON's Notes:  To obtain the correct NPS using this tool use a 
;linear noise flat field with no backgrounds or fixed pattern noise.  
;Subtract the mean value of the noise field from the noise image (fn).   
;i.e. use zero mean data. Finally, DIVIDE fn BY THE MEAN VALUE SUBTRACTED
;ABOVE (to check use an image of the form i = randomu(s,512,512)*10+100 
;i.e. an image with 100 x-ray quanta should yeild a DQE of 1.0.
;Be sure to include the image plane pixel size in mm!
;This process will yeild nps where DQE ~ 1/(NPS x Nquanta) where nquanta
;is the number of x-ray quanta per square mm in the detector.

; calculates radial nps using Gaussian pyramid
; using kernels convolved with decimated versions of image
; fn = input noise image with pixel spacing dpix (mm)
; sp = output spectrum at radial frequency values freq
; with estimated rms dev of sperr

time0 = systime(1)  ; for timing

dpix = dpixin

nfrq = 7   ; number of freqs
sp = make_array(nfrq)
sperr = make_array(nfrq)
freq = make_array(nfrq)

f = fn                     ; store input image in working array
res=size(f)
if res[0] ne 2 then print,"SPEC_CALC: not an image"
ni = res[1]
nj = res[2]
;print, "size of image: ni,nj",ni,nj
; assume square image, i.e., nj tacitly assumed to be same as ni
print, "ni, pixel spacing",ni,dpix

i = 0
gen_b2,kern,ki,kj    ; generate kernel for highest freq.

kie=(ki+1)/2  ; useless edge width + 1 after convolution
kje=(kj+1)/2
g0f = convol(f,kern,/edge_truncate)
;print,size(kf)
;fconv = ff[kie-1:ni-kie,kje-1:nj-kje]
f = f - g0f  ; original image - convolution of orignal with b2
varf = variance(f[kie-1:ni-kie,kje-1:nj-kje])

sp[i] = varf*dpix*dpix  ; full-sized image, l2n
freq[i] = 0.92/(2*dpix)
; approximate rms uncertainty
gen_l2,kdif    ; generate kernel for normalization and uncertainty
sp[i] = sp[i]/total(kdif*kdif)
fact=sqrt(total(conv2d(kdif,kdif)^2))/total(kdif^2)
sperr[i] = (2*sp[i]*fact)/sqrt(2.*(ni-2*kie+2)*(nj-2*kje+2))

print, "ni, freq, nps, nps_err", ni, freq[i], sp[i], sperr[i]
if i ge nfrq-1 then goto, ending

i = 1
gen_b2,kern
; two convolutions with b2 same as one with b4
g1f = convol(g0f,kern,/edge_truncate)
l1f = g0f - g1f
print,"mean g1f",mean(g1f)
print,"mean l1f",mean(l1f)
kie = 2*kie - 1
kje = 2*kje - 1
varf = variance(l1f[kie-1:ni-kie,kje-1:nj-kje])
sp[i] = varf*dpix*dpix  ; full-sized image, l2n
freq[i] = 0.56/(2*dpix)
; approximate rms uncertainty
gen_l4,kdif    ; generate kernel for uncertainty
print,total(kdif*kdif)
sp[i] = sp[i]/total(kdif*kdif)
fact=sqrt(total(conv2d(kdif,kdif)^2))/total(kdif^2)
sperr[i] = (2*sp[i]*fact)/sqrt(2.*(ni-2*kie+2)*(nj-2*kje+2))

print, "ni, freq, nps, nps_err", ni, freq[i], sp[i], sperr[i]
if i ge nfrq-1 then goto, ending

gen_b2,kern,ki,kj    ; generate b2
kie=(ki+1)/2  ; useless edge width + 1 after convolution
kje=(kj+1)/2
kie = 2*kie - 1  ; for b4 & l4 (two colvolutions by b2)
kje = 2*kje - 1
spnorm=0.95/(0.152262)  ; nps normalization for remainder
spfact=1.37  ; factor used to account for higher level filters
gfprev = g1f
for i = 2,nfrq-1 do begin
  ni=ni/2                    ; decimate image by facter of two
  nj=nj/2
  dpix = 2*dpix
  if ni ge ki then begin    ; check that image is larger than kernel
  gftmp = rebin(gfprev,ni,nj,sample=1)  ; decimate previous image
  print,"rebin gfprev, continue",mean(gftmp)
  gftmp2 = convol(gftmp,kern,/edge_truncate) ; effectively convolved with b2
  gfprev = convol(gftmp2,kern,/edge_truncate) ; effectively convolved with b4
  gftmp = gftmp - gfprev  ; difference unconv., convov b4
  varf = variance(gftmp[kie-1:ni-kie,kje-1:nj-kje])
  sp[i] = varf*dpix*dpix  ; full-sized image, l2n
  freq[i] = 0.68/(2*dpix)
  ; approximate rms uncertainty
  print,"spnorm",spnorm
  sp[i] = sp[i]*spnorm
  spnorm=spnorm*spfact
  spfact=(spfact)^0.23
  fact=sqrt(total(conv2d(kdif,kdif)^2))/total(kdif^2)
  sperr[i] = (2*sp[i]*fact)/sqrt(2.*(ni-2*kie+2)*(nj-2*kje+2))

  print, "ni, freq, nps, nps_err", ni, freq[i], sp[i], sperr[i]

  endif
endfor
ending:
print, "time (sec) = ",systime(1) - time0
print, " "

; display spectrum
; display black on white
print, "nps",sp
print, "nps_err",sperr
print, "freq",freq

plot, freq, sp, psym = 7,/xlog,/ylog, $
   xtitle = 'SPATIAL FREQUENCY (mm^-1)', $
   ytitle = 'NOISE POWER SPECTRUM (mm^2)'
oploterr, freq, sp, sperr

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function conv2d, a, b

; calculates two 2D arrays a and b

res = size(a)
nax = res(1)
nay = res(2)
res = size(b)
nbx = res(1)
nby = res(2)
;print, "convolve: input array sizes", $
;  nax," X",nay," and ",nbx," X",nby
ncx = nax + nbx - 1
ncy = nay + nby - 1
c = make_array(ncx,ncy)      ; result

for ja = 0,nay-1 do for ia = 0,nax-1 do begin
  for jb = 0,nby-1 do for ib = 0,nbx-1 do begin
    jc = ja + jb
    ic = ia + ib
    c(ic,jc) = c(ic,jc) + a(ia,ja)*b(ib,jb)
  endfor
endfor
;print, "output array ",ncx," X",ncy
;print,c
return, c
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro gen_b2, kern, ki, kj

; create kernel: 3X3 binomial, normalized to unit sum
; output is kernel image and its size

b1 = [[1,1],[1,1]]
b1 = b1/total(b1)
b2 = conv2d(b1,b1)

kern = b2
res=size(kern)
ki = res[1]
kj = res[2]

end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function mean, x

; returns mean value of x

mom = moment(x)
return, mom[0]

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function variance, x

; returns variance of x

mom = moment(x)
return, mom[1]

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro gen_l2, kern

; create kernel 3X3 Laplacian

b1 = [[1,1],[1,1]]
b1 = b1/total(b1)
;print, b1

b2 = conv2d(b1,b1)
;print, b2

l2 = [[0,0,0],[0,1,0],[0,0,0]]
l2 = l2 - b2
;print,"L2"
;print,l2

kern = l2

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro gen_l4, kern

; create kernel 5X5 Laplacian

b1 = [[1,1],[1,1]]
b1 = b1/total(b1)
;print, b1

b2 = conv2d(b1,b1)
;print, b2

b4 = conv2d(b2,b2)
;print, b4

l4 = b4
l4[1:3,1:3] = l4[1:3,1:3] - b2
l4=-l4
;print,"L4"
;print,l4

kern = l4

end