ncl做带通滤波求助

Yoger

最近在学习用NCL做带通滤波，在官网上下载了filter_3.ncl，用2019年逐日位势高度场资料，做10-90天的带通滤波，但是过滤后得到的滤波结果数值量级和原始场相差甚远，请问各位前辈，这是怎么回事？附上ncl代码：
begin

;  diri    = "/Users/shea/Data/AMWG/"
;  diri    = "./"
  vName   = "hgt"        ; name of variable on the file
  fili    = "hgt.2019.nc"
  f       = addfile(fili, "r")
  hgt        =        f->$vName$
  printVarSummary(hgt)
  x       = f->$vName$(:,{500},{0},{120})
  print(x)

; ***********************************************
; create the filter weights and apply
; ***********************************************

  ihp     = 2                             ; band pass
  sigma   = 1.0                           ; Lanczos sigma

  nWgt    = 201                           ; loose 100 each end                           
  fca     = 1./90.                       ; start freq
  fcb     = 1./10.                       ; last  freq
  wgt     = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma )
  printVarSummary(wgt)
  print(wgt)
  xBPF    = wgt_runave ( x, wgt, 0 )      ; 10-90 day
  printVarSummary(xBPF)
  print(xBPF)

; ***********************************************
; create new date array for use on the plot
; ***********************************************

  date    = cd_calendar(x&time,-2)        ; yyyymmdd
  yrfrac  = yyyymmdd_to_yyyyfrac (date, 0)
  delete(yrfrac@long_name)
  delete(x@long_name)

  pStrt = 20190101      ; 4 years: winter 96-97 MJO gold standard
  pLast = 20191231
  iStrt = ind(date.eq.pStrt)        ; user specified dates
  iLast = ind(date.eq.pLast)
  delete(date)

  pltType = "png"                    ; send graphics to PNG file
  pltName = "filters"
  plot    = new ( 2, "graphic")
  wks     = gsn_open_wks (pltType,pltName)

  res              = True            ; plot mods desired
  res@gsnDraw      = False           ; don't draw
  res@gsnFrame     = False           ; don't advance frame yet      

  res@vpHeightF    = 0.4             ; change aspect ratio of plot
  res@vpWidthF     = 0.8                  
  res@vpXF         = 0.1             ; start plot at x ndc coord
  res@gsnYRefLine  = 0.0             ; create a reference line   

  res@gsnCenterString = "500hPa Height [0, 120E]"
  res@tmXBFormat = "f"
  plot(0) = gsn_csm_xy (wks,yrfrac(iStrt:iLast),x(iStrt:iLast),res)      

  res@gsnCenterString  = "Band Pass Filtered: 10-90 day"     
  plot(1) = gsn_csm_xy (wks,yrfrac(iStrt:iLast),xBPF(iStrt:iLast),res)     

  resP                 = True   
  resP@gsnMaximize     = True
  gsn_panel(wks,plot,(/2,1/),resP)

  ; Create XY plot of frequency versus response

  xyres = True

  xyres@gsnMaximize    = True
  xyres@gsnFrame       = False
  xyres@tiMainString   = "Band Pass: 10-90 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)

  frame(wks)
end

oyy.

没毛病，带通滤波就相当于对它求距平。
图中黑线是原始数据，红色是滤波，绿线是减去多年平均的距平值。

喵鳕鱼

学习了，最近也遇到这个问题

Yoger

oyy. 发表于 2020-2-25 18:47
没毛病，带通滤波就相当于对它求距平。
图中黑线是原始数据，红色是滤波，绿线是减去多年平均的距平值。

多谢前辈，已经懂了

niuyingli

想请教一下这个10-90天是怎么选的，也就是fca,fcb

wsx

niuyingli 发表于 2021-8-25 10:02
想请教一下这个10-90天是怎么选的，也就是fca,fcb

一般是先进行功率谱或者小波分析得到周期，然后针对这个周期来滤波

什夜0.0

你好，请问输入这个程序就直接出图片了吗，我也用了这个程序，但为什么没出图片呢，也没有报错