;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
undef("demod_cmplx")
function demod_cmplx(y:numeric, frqdem[1]:numeric, frqcut[1]:numeric, nwt[1]:integer,  ndim[1]:integer,  opt[1]:logical)
;
; Simple complex demodulation 
; Extract approximately periodic components (ie, slowly varying) from a time series
;
; Source: http://www.uni-muenster.de/ZIV.BennoSueselbeck/s-html/helpfiles/demod.html
;   Complex demodulation is a technique for analyzing nonstationary time series by estimating 
;   the instantaneous amplitude and phase of a given harmonic component. To better understand 
;   the results of complex demodulation several lowpass filters should be tried: the smaller 
;   the pass band, the less instantaneous in time but more specific in frequency is the result.
;
; Nice summaries:
; http://faculty.washington.edu/kessler/generals/complex-demodulation.pdf
; http://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2007/pdfHMANlJLqCE.pdf
;
local dimy, ranky, pi, pi2, zero, one, frqdem25, frqdem50, ny, arg, ARG, yr, ys \
    , delf, ihp, nsigma, fcb, wgt, amp, pha
begin
   if (typeof(y).eq."double" .or. typeof(frqdem).eq."double" \
                             .or. typeof(frqcut).eq."double" ) then
       pi       = 4d0*atan(1d0)
       zero     = 0.0d
       one      = 1.0d
       frqdem25 = 0.25d
       frqdem50 = 0.50d
   else
       pi       = 4.0*atan(1.0)
       zero     = 0.0 
       one      = 1.0 
       frqdem25 = 0.25 
       frqdem50 = 0.50 
   end if
   pi2 = 2*pi

   dimy  = dimsizes(y)
   ranky = dimsizes(dimy) 
   ny    = dimy(ndim) 
   delf  = one/ny

; error check

   if (frqcut.le.delf .or. frqcut.gt.frqdem50 .or. \
       frqdem.le.delf .or. frqdem.gt.frqdem50      ) then
       print("demod_cmplx: illegal valu(s): frqcut="+frqcut+" : frqdem="+frqdem)
       exit
   end if

   if (frqcut.gt.frqdem) then 
       print("demod_cmplx: frqcut > frqdem: frqcut="+frqcut+" : frqdem="+frqdem)
       exit
   end if

; demodulate (primary): Bloomfield, 1976: page 147

   arg = ispan(0,ny-1,1)*frqdem*pi2
   ARG = conform(y, arg, ndim)     ; all dimensions
                                   ; complex demodulate 
   yr  =  y*cos(ARG)*2             ; 'real' series
   yi  = -y*sin(ARG)*2             ; 'imaginary'       

;---Create weights for Lanczos low pass filter

   ihp    = 0                      ; low pass 
   nsigma = one                       
   if (opt .and. isatt(opt, "nsigma")) then
       nsigma := totype(opt@nsigma, typeof(pi))
   end if 

   fcb    = totype(-999, typeof(pi))
   wgt    = filwgts_lanczos (nwt, ihp, frqcut, fcb, nsigma)  

;---Apply filter weights to raw demodulated series (nwt/2 point 'lost; at start/end)

   yr = wgt_runave_n_Wrap ( yr, wgt, 0, 0) 
   yi = wgt_runave_n_Wrap ( yi, wgt, 0, 0) 
  
;---Use the low pass filtered values
   
   amp  = sqrt(yr^2 + yi^2)
   pha  = where(amp.eq.zero, zero, atan2(yi, yr))    ; primary phase 

   pha2 = conform(pha, zero, -1)                     ; secondary phase (Bloomfield)
   pha2 = where(pha.ge.zero, pha-pi2, pha+pi2)   

   amp@long_name  = "amplitude"
   pha@long_name  = "phase"         
   pha2@long_name = "secondary phase"         
   if (isatt(y,"units")) then
       amp@units = y@units
   end if

   copy_VarCoords(y, amp )
   copy_VarCoords(y, pha )
   copy_VarCoords(y, pha2)

   if (opt .and. isatt(opt, "phase_12")) then
       return([/ amp, pha, pha2 /])     ; return as a list
   else
       return([/ amp, pha /])           ; return as a list
   end if 
end