undef("SPECX_ANAL")
function SPECX_ANAL(X[*]:numeric, iopt[1]:integer, jave[1]:integer, pct[1]:numeric)
;
; all-NCL version of NCL's specx_anal
;
; The input 'x' array may return anomalies sepending upon 'iopt'
;
local x, N, xVari, xVaro, cf, cr, ci, px, spcx, total_area, df, frq
    , wgts, sclVar, sdof, acr, tapcf
begin
   N     = dimsizes(X)
   xVari = variance(X)*(N-1.)/N        ; input variance
  
   if (iopt.ge.0) then                 ; remove mean
       x = X - avg(X)
   else
       x = X
   end if

   if (iopt.ge.1) then
       x = dtrend(x,False)             ; step 1: detrend
   end if
                                       ; step 2: Calculate input sample variance 
   xVaro = variance(x)*(N-1.)/N        ; population variance after detrending

   x     = taper (x,pct,0)             ; step 3: taper the series
                                       ; tapering correction factor
   tapcf = 0.5*(128-93*pct)/(8-5*pct)^2; See: Bloomfield: Intro to Time Series
  
   cf    = ezfftf(x)                   ; step 4: perform forward FFT
   cr    = cf(0,:)                     ; clarity ... real coef
   ci    = cf(1,:)                     ;             imag coef
   delete( cf )

   px    = cr^2 + ci^2                 ; step 5: periodogram of "x"
  
   wgts         = fspan(1,1,jave)
   wgts(0)      = 0.5                  ; not sure why I did this
   wgts(jave-1) = 0.5
   wgts         = wgts/sum(wgts)       ; normalize sum of wgts to one 

   spcx  = wgt_runave(px,wgts,1)       ; weighted run ave

   delete(px)
  
  ;step 7: normalize the area under the curve [6] to the variance.
  ;       Remember the first and last bandwidths are only 0.5*df
   
   df         = 1./N 
   total_area = (spcx(0) + spcx(N/2-1))*(df/2) + sum(spcx(1:N/2-2))*df
   sclVar     = xVaro/total_area
   spcx       = spcx*sclVar
  
   frq        = fspan(df,0.5,N/2)
  
   sdof       = 2/(tapcf*sum(wgts^2))  ; sum squares of normalized wgts 
                                       ; to calculate dof 
   sdof@bw    = 0.5*sdof*df            ; bw: jones pg 202: murphy+katz
   sdof@xvari = xVari
   sdof@xvaro = xVaro
   sdof@frq   = frq
   sdof@spcx  = spcx
   acr        = esacr(x,1)             ; lag-1 autocorrelation                 
   sdof@xlag1 = acr(1) 
  
   return( sdof )
end