finfo = fourier_info(data, nhx, sclPhase) ;返回数组 finfo[3][...][N/2]
;对data的最右边一维做傅里叶变换，返回每个周期对应的振幅、第一个极大值的位置、方差百分比
;其中，振幅即每个周期对应的傅里叶系数的平方和的均方根
;设data最右边一维的长度为N，则分解得到的每个周期中包含的谐波数范围为 (0, N/2]，谐波数值越大，周期越短。
;nhx=0时会返回谐波数分别为 (0, N/2] 时所对应的周期的信息，若 0 < nhx <= N/2，则返回谐波数分别为 (0, nhx] 时所对应的周期信息
;sclPhase其实就是data最右边一维点与点之间的间隔数，若最右边一维是经度，有72个点，则 sclPhase = 360/72 = 5

;注意点：data不能有缺测值，不包括周期点
;此外，该函数也有 fourier_info_n 的形式

cf = ezfftf(data) ;返回数组 [2][...][N/2]
;对data最右边一维做傅里叶变换，设其最右边一维的长度为N，返回傅里叶的实部与虚部系数
;同时以属性 xbar 返回data最右边一维的均值，以 npts 返回data最右边一维的长度，最右边一维不需要是2的倍数
;如果有缺测值，计算出来的傅里叶系数为0

data = ezfftb(cf, xbar) 
;cf 即利用 ezfftf 生成的傅里叶系数，其最左边维的长度必须是2，且 cf(0,...) 为实部，cf(1,...) 为虚部
;xbar 是data最右边维的平均数，可以是一个数，也可以是一维数组，长度为 data 左边维数长度的乘积

;利用该函数及 ezfftf 函数可以提取出 data 前几个波数周期的合成
;或者只提取出其中一个波动的曲线
;设现在有一个 data (ntime, nlev, nlat, nlon)
cf = ezfftf (data) ;得到 cf (2, ntime, nlev, nlat, nlon/2)
cf(:,:,:,:,3:nlon/2) = 0.0
data_wave1_3 = ezfftb (cf, 0)
;重构的data_wave1-3的平均数为0
data_wave1_3 = ezfftb (cf, cf@xbar)
;重构的data_wave1-3的平均数为cf@xbar

(spcx(0)+spcx(N/2-1))*(df/2) + SUM{spcx(1:N/2-2)*df} = total variance of the series 

d   = 0     ; detrending opt: 0=>remove mean 1=>remove mean + detrend
sm  = 1     ; smooth: should be at least 3 and odd; if small than 3, it would do no smoothing
pct = 0.10  ; percent taper: (0.0 <= pct <= 1.0) 0.10 common. If pct =0.0, no tapering will be done. If pct = 1.0, the whole series is affected

  ;************************************************
  ; calculate mean spectrum spectrum and lag1 auto cor
  ;************************************************
  
  ; loop over each segment of length ntim
  
   spcavg = new ( ntim/2, typeof(x))
   spcavg = 0.0
  
   r1zsum = 0.0
  
   do nt=0,nyear-1
      dof    = specx_anal(x(nt,:),d,sm,pct)      ; current segment spc
      spcavg = spcavg + dof@spcx                ; sum spc of each segment
      r1     = dof@xlag1                        ; extract segment lag-1
      r1zsum = r1zsum  + 0.5*log((1+r1)/(1-r1)) ; sum the Fischer Z
   end do
  
   r1z  = r1zsum/nyear                          ; average r1z
   r1   = (exp(2*r1z)-1)/(exp(2*r1z)+1)         ; transform back
                                                ; this is the mean r1
   spcavg  = spcavg/nyear                       ; average spectrum
  
  ;************************************************
  ; Assign mean spectrum to data object
  ;************************************************
  
   df      = 2.0*nyear   ; deg of freedom
   df@spcx = spcavg      ; assign the mean spc
   df@frq  = dof@frq
   df@xlag1= r1          ; assign mean lag-1
  
   splt = specx_ci(df, 0.05, 0.95)   ; confidence interval

;splt(4,ntim/2) including input spectrum (indx=0), Markov Red Noise spectrum (indx=1), 
;lower confidence bound for Markov (indx=2), upper confidence bound for Markov (indx=3)

  xyres = True
  xyres@gsnMaximize    = True
  xyres@gsnFrame       = False
  xyres@tiMainString   = "Band Pass: 20-100 srate=1: sigma = " + sigma
  xyres@tiXAxisString  = "frequency"
  xyres@tiYAxisString  = "response"
  xyres@gsnLeftString  = "fca=" + fca + "; fcb=" + fcb
  xyres@gsnRightString = nWgt
  xyres@trXMaxF        = 0.1
  xyres@trYMaxF        = 1.1
  xyres@trYMinF        = -0.1
  plot = gsn_csm_xy(wks, wgt@freq, wgt@resp, xyres)

  X = (/0.0, fca, fca, fcb, fcb, 0.1/)      ; ideal filter
  Y = (/0.0, 0.0, 1.0, 1.0, 0.0,  0.0 /) 
  resGs = True
  resGs@gsLineThicknessF = 1.0
  gsn_polyline(wks,plot,X,Y,resGs)

    ua    = f->U_anom(:,{LAT},{LON})   ; ua(time), read from one grid point

    ca    = 50.0        ; band start (longer period)
    cb    = 40.0        ; band end （时间间隔不能等于2）

    fca   = 1.0/ca      ; 'left'  frequency
    fcb   = 1.0/cb      ; 'right' frequency

    dims  = 0           ; 'time' dimension of ua  

    opt   = True        ; options to set
    opt@return_envelope = True ; time series of filtered and envelope values

    ua_bf = bw_bandpass_filter (ua,fca,fcb,opt,dims)       ; (ua,fca,fcb,opt,dims)
    copy_VarMeta(ua, ua_bf)
    ua_bf@long_name = "Band Pass: "+cb+"-"+ca+" day"