[ncl-talk] problem with wavelet power calculation

[Dennis Shea]

Well, one obvious issue is that the 1st and last two values are [I assume]
missing [_FillValues] values.

197901 9969209968386869046778552952102584320.000000
197902 9969209968386869046778552952102584320.000000
197903 10.643162
197904 37.253376
[SNIP]
201709 10.624535
201710 1.700483
201711 9969209968386869046778552952102584320.000000
201712 9969209968386869046778552952102584320.000000
===

The underlying function code [fortran] uses an FFT. Missing values are not
allowed in FFT functions.  Hence, they must be eliminated.

The attached eliminates those values.
*%> ncl YanYan_wavelet.ncl*


Note that setting the assorted input values may be more art-than-science.
Setting these appropriately is the user's responsibility. It may require
some iteration.
ncl-talk is not an expert on wavelets.

Good luck








On Tue, Feb 4, 2020 at 8:40 PM Yan Yan via ncl-talk <ncl-talk at ucar.edu>
wrote:

> Hi there,
>
> I ran into some problems lately when using wavelet analysis in NCL. Here
> they are.
>
> I was using wavelet to analyze an EOF PC1, which is in monthly resolution.
> Two problems stuck out.
> One is when I dtrend it, the result time series became constant,
> -8.520693e+34. The following wavelet gave me a time series with power being
> all "inf". I couldn't plot the power spectrum. This was obvious the input
> problem, so I changed the input into the PC1 without detrending, which gave
> me the same problem, the power of the wavelet time series still showed
> "inf".
>
> I suspect there is something wrong with this input time series, but I
> don't know enough to identify what the problem is, and how to deal with it.
> PC1 data is attached. Can somebody point me to the right direction on this?
>
> Also, the following is the code I used, which is basically using NCL
> example code. Since data is monthly, I tried 1. and 1./12 for dt, both gave
> me the same "inf" power.
>
> wts = eof_ts(0,:)
> ; wts = dtrend(wts,False)
>   mother  = 0
>   param   = 6.0
>   dt      = 1.
>   s0      = dt
>   dj      = 0.25
>   jtot    = 1 + floattointeger(((log10(N*dt/s0))/dj)/log10(2.))
>   npad    = N
>   nadof   = new(2,float)
>   noise   = 1
>   siglvl  = .05
>   isigtest= 1
>
>   w =
> wavelet(wts,mother,dt,param,s0,dj,jtot,npad,noise,isigtest,siglvl,nadof)
> ;************************************
> ; create coodinate arrays for plot
> ;************************************
>   power            = onedtond(w at power,(/jtot,N/))
>   power!0          = "period"                        ; Y axis
>   power&period     = w at period
>
>   power!1          = "time"                          ; X axis
>   power&time       = time
>
>   power at long_name  = "Power Spectrum"
>   power at units      = "unit^2"
>
> ; compute significance ( >= 1 is significant)
>   SIG              = power                            ; transfer meta data
>   SIG              = power/conform (power,w at signif,0)
>   SIG at long_name    = "Significance"
>   SIG at units        = " "
>
> ;************************************
> ; initial resource settings
> ;************************************
>   wks = gsn_open_wks("x11","wavelet")             ; send graphics to PNG
> file
>   plot = new(2,graphic)                           ; create graphical array
>
>   res                     = True                  ; plot mods desired
>   res at gsnFrame            = False                 ; don't advance frame
> yet
>   res at gsnDraw             = False                 ; don't draw yet
>
>   res at cnFillOn            = True                  ; turn on color
>   res at cnFillPalette       = "BlAqGrYeOrReVi200"   ; set color map
>   res at cnFillMode          = "RasterFill"          ; turn on raster mode
>   res at cnRasterSmoothingOn = True                  ; turn on raster
> smoothing
>   res at cnLinesOn           = False                 ; turn off contour lines
>   res at lbOrientation       = "Vertical"            ; vertical label bar
>
>   res at trYReverse          = True                  ; reverse y-axis
>   res at tmLabelAutoStride   = True
>
>   res at vpHeightF           = .4                    ; height and width of
> plot
>   res at vpWidthF            = .85
>
>
>   res at cnLevelSelectionMode = "ManualLevels"       ; set manual contour
> levels
>   res at cnMinLevelValF       = 0.0                  ; set min contour level
>   res at cnMaxLevelValF       = 20.0                 ; set max contour level
>   res at cnLevelSpacingF      = 2.0                  ; set contour spacing
>
> ; power
>   res at gsnLeftString       = "Wavelet Power"
>   plot(0) = gsn_csm_contour(wks,power({0:35},0:N-1),res)
>
> ; significance
>   res at cnMinLevelValF      = 1.0                ; set min contour level
>   res at cnMaxLevelValF      = 4.0                ; set max contour level
>   res at cnLevelSpacingF     = 3.0                ; set contour spacing
>   res at gsnLeftString       = "Wavelet Significance"
>   plot(1) = gsn_csm_contour(wks,SIG({0:35},0:N-1),res)
>
> Any help is very much appreciated!
>
> Best,
> Yan
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20200205/7d5d14ab/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: YanYan_wavelet.ncl
Type: application/octet-stream
Size: 2447 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20200205/7d5d14ab/attachment.obj>