Under the Hood

We strongly recommend everyone to follow the procedure as instructed here when using WaLSAtools — a user-friendly tool — which gives you all information you need to do your analysis. However, for experts who want to make themselves familiar with the techniques and codes under the hood, inspect them and modify/develop/improve them, some of the main codes are also provided below. Please note that all codes and their dependencies are available in the GitHub repository.

Spectral Analyzer

WaLSA_speclizer

This code computes power spectrum and its statistical significance level for a 1D signal (or all pixels of an image sequence, i.e., a 3D cube) using FFT (Fast Fourier Transform), Lomb-Scargle, Wavelet, and HHT (Hilbert-Huang Transform) analysis techniques. In addition, the code can output mean power spectrum (averaged over power spectra of several pixels) as well as dominant frequency and power using the above-mentioned analysis methods.

WaLSA_speclizer.pro

; -----------------------------------------------------------------------------------------------------
; WaLSAtools: Wave analysis tools
; Copyright (C) 2025 WaLSA Team - Shahin Jafarzadeh et al.
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
; 
; http://www.apache.org/licenses/LICENSE-2.0
; 
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
; 
; Note: If you use WaLSAtools for research, please consider citing:
; Jafarzadeh, S., Jess, D. B., Stangalini, M. et al. 2025, Nature Reviews Methods Primers, in press.
; -----------------------------------------------------------------------------------------------------


;+
; NAME: WaLSA_speclizer: WaLSA Spectral Analyzer
;       part of -- WaLSAtools --
;
; PURPOSE:
;   Compute power spectrum and its statistical significance level for a 1D signal 
;           (or all pixels of an image sequence, i.e., a 3D cube) using 
;           FFT (Fast Fourier Transform), Lomb-Scargle, Wavelet, or 
;           HHT (Hilbert-Huang Transform) analyses.
;           -- Signals are detrended (linearly or using higher-order polynomial fits) and 
;              apodized (using a Tukey window) prior to the spectral analysis (unless otherwise it is omitted).
;           -- Power (and significance levels) are returned in DN^2/mHz, frequencies in mHz.
;
; CALLING SEQUENCE:
;   EXAMPLES:
;   power = walsa_speclizer(cube,time,mode=1,/fft,frequency=frequency,significance=significance,siglevel=0.01)
;   power = walsa_speclizer(cube,time,mode=1,/wavelet,/global,frequency=frequency,significance=significance)
;
; + INPUTS:
;   data:           1D time series, or (x,y,t) datacube, any type
;                   (an ordered sequence of data points, typically some variable quantity measured at successive times.)
;   time:           observing times of the time series in seconds
;
; + OPTIONAL KEYWORDS:
; ----type of analysis----
;   fft:            if set, Fast Fourier Transform (FFT) power spectrum is computed: for regular (evenly sampled) time series.
;   lombscargle:    if set, Lomb-Scargle power spectrum is computed: for irregular (unevenly sampled) time series.
;   hht:            if set, a power spectrum associated to EMD + Hilbert Transform is computed: for regular (evenly sampled) time series.
;   wavelet:        if set, Wavelet power spectrum is computed (default: Morlet function with omega=6): for regular (evenly sampled) time series.
;   welch:          if set, Welch power spectrum is computed
; ----padding, detrending, and apodization parameters----
;   padding:        oversampling factor: zero padding (increasing timespan) to increase frequency resolution (NOTE: doesn't add information)
;   apod:           extent of apodization edges (of a Tukey window); default 0.1
;   nodetrendapod:  if set, neither detrending nor apodization is performed!
;   pxdetrend:      subtract linear trend with time per pixel. options: 1=simple, 2=advanced; default: 2
;   polyfit:        the degree of polynomial fit to the data to detrend it.
;                   if set, instead of linear fit this polynomial fit is performed.
;   meantemporal:   if set, only a very simple temporal detrending is performed by subtracting the mean signal from the signal.
;                   i.e., the fitting procedure (linear or higher polynomial degrees) is omitted.
;   meandetrend:    if set, subtract linear trend with time for the image means (i.e., spatial detrending)
;   recon:          optional keyword that will Fourier reconstruct the input timeseries.
;                   note: this does not preserve the amplitudes and is only useful when attempting 
;                   to examine frequencies that are far away from the 'untrustworthy' low frequencies.
;   resample         if recon is set, then by setting resample, amplitudes are scaled to approximate actual values.
; ----significance-level parameters----
;   siglevel:       significance level (default: 0.05 = 5% significance level = 95% confidence level)
;   nperm:          number of random permutations for the significance test -- the larger the better (default: 1000)
;   nosignificance: if set, no significance level is calculated.
; ----power calibration----
;   mode:           outputted power mode: 0 = log(power), 1 = linear power (default), 2 = sqrt(power) = amplitude
; ----wavelet parameters/options----
;   mother:         wavelet function, providing different localisation/resolution in frequency and in time (also depends on param, m).
;                   currently, 'Morlet','Paul','DOG' (derivative of Gaussian) are available. default: 'Morlet'.
;   param:          optional mother wavelet parameter. 
;                   For 'Morlet' this is k0 (wavenumber), default is 6. 
;                   For 'Paul' this is m (order), default is 4.
;                   For 'DOG' this is m (m-th derivative), default is 2 (i.e., the real-valued Mexican-hat wavelet)
;   dj:             spacing between discrete scales. default: 0.025
;   global:         only if wavelet is set: returns global wavelet spectrum (averaged over time domain)
;   oglobal:        global wavelet spectrum excluding regions influenced by cone-of-influence (CoI; regions subject to edge effect)
;   rgws:       time-integral of wavelet power excluding regions influenced by cone-of-influence and only for those above the confidence level
;                   this returns power-weighted frequency distribution (with significant power & unaffected by CoI)
;                   Note: this is likely the most correct spectrum!
;   colornoise:     if set, noise background is based on Auchère et al. 2017, ApJ, 838, 166 / 2016, ApJ, 825, 110
; ----HHT parameters/options----
;   stdlimit:       standard deviation to be achieved before accepting an IMF (recommended value between 0.2 and 0.3; perhaps even smaller); default: 0.2
;   nfilter:        Hanning window width for two dimensional smoothing of the Hilbert spectrum. default: 3 
;                   (an odd integer, prefrabely equal to or larger than 3; equal to 0 to avoid the windowing)
;   emd:            if set, intrinsic mode functions (IMFs) and their associated frequencies (i.e., instantaneous frequencies) are outputted
; ----dominant frequency----
;   nodominantfreq: if set, dominant frequency and dominant power are not calculated (to, e.g., save computational time for large datasets)
;
; + OUTPUTS:
;   power:          a 1D array of power (or a 3D array if the input is a 3D cube).
;                   the only exception is for wavelet (where global is not set).
;                   power is divided by the first (non-zero) frequency. unit: DN^2/mHz
;   significance:   significance levels, with the same dimension as the power. unit: DN^2/mHz
;   frequency:      an array of frequencies, with the same size as the power. unit: mHz
;   period:         1D array of periods (in seconds)
;   coicube:        cone-of-influence cube, only when wavelet analysis is performed --> if wavelet is set
;   imf:            the intrinsic mode functions (IMFs) from EMD alalysis within the HHT --> if hht and emd are set
;   instantfreq:    instantaneous frequencies of each component time series --> if hht and emd are set
;   dominantfreq:   dominant frequency, i.e., frequency corresponding to the maximum power (in mHz): same saptial size as input data (i.e., 1D or 2D)
;                   note: if there are multiple peaks with the exact same power, the lowest dominant frequency is returned!
;   dominantpower:  power (in DN^2/mHz) corresponding to the dominant frequency: same saptial size as input data (i.e., 1D or 2D)
;   rangefreq:      frequency range over which the dominant frequency is computed. default: full frequency range
;   averagedpower:  spatially averaged power spectrum (of multiple 1D power spectra). unit: DN^2/mHz
;   amplitude:      a 1D array of oscillation amplitude (or a 3D array if the input is a 3D cube).
;
; MODIFICATION HISTORY
;
;  2010 plotpowermap: Rob Rutten, assembly of Alfred de Wijn's routines 
;                     (https://webspace.science.uu.nl/~rutte101/rridl/cubelib/plotpowermap.pro)
;  2014-2021 Extensively modified/extended by Shahin Jafarzadeh, with contributions from Marco Stangalini and David B. Jess
;-

;---------------------------- HHT (EMD + Hilbert) -----------------------------
function getpowerHHT,cube,cadence,stdlimit,nfilter=nfilter,significance=significance,siglevel=siglevel,nperm=nperm,padding=padding,$
                     frequencies=frequencies,nosignificance=nosignificance,emd=emd,imf=imf,instantfreq=instantfreq,averagedpower=averagedpower,$
                     dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower,amplitude=amplitude,$
                     originalcube=originalcube,apod=apod,nodetrendapod=nodetrendapod,pxdetrend=pxdetrend,meandetrend=meandetrend,polyfit=polyfit,$
                     meantemporal=meantemporal,recon=recon,resample_original=resample_original,silent=silent
  ; Hilbert-Huang Transform (HHT) power spectra

  if padding gt 1 then begin ; zero padding (optional): to increase frequency resolution
      if silent eq 0 then begin
          print, ' '
          print, ' -- Zero Padding (oversampling factor: '+strtrim(padding,2)+') .....'
          print, ' '
      endif
      nx = n_elements(cube[*,0,0])
      ny = n_elements(cube[0,*,0])
      nt = n_elements(cube[0,0,*])      
      padded=fltarr(nx,ny,padding*nt)
      mid_point = ROUND((padding*nt)/2.)
      lower_point = mid_point-nt/2.
      upper_point = mid_point+nt/2.-1
      padded[*,*,mid_point-nt/2.:mid_point+nt/2.-1] = cube
      cube = padded
  endif
  sizecube=size(cube)
  nx=sizecube[1]
  ny=sizecube[2]
  nt=sizecube[3]

  if silent eq 0 then begin
      print, ' '
      print, ' ...... output Marginal HHT Spectra '
      print, ' '
  endif

  dt = cadence
  IMFcal = walsa_emd_function(reform(cube[0,0,*]),stdlimit,dt=dt)
  hhs = walsa_hilbert_spec(IMFcal,dt,freq=frequencies, marginal=pm, nfilter=nfilter)
  nff= n_elements(frequencies)
  if frequencies[0] eq 0 then begin
    frequencies = frequencies[1:nff-1]
    pm = pm[1:nff-1]
    f0 = 1
  endif else f0 = 0
  nf= n_elements(frequencies)
  frequencies = frequencies*1000. ; in mHz

  powermap=fltarr(nx,ny,nf) ; HHT power spectra
  amplitude=fltarr(nx,ny,nf)

  if emd then begin
      imf = fltarr(nx,ny,nt,20) ; 20: maximum number of IMFs that can be created
      instantfreq = fltarr(nx,ny,nt,20)
  endif
  if nosignificance eq 0 then significance=fltarr(nx,ny,nf) ; significancec cube
  if nodominantfreq eq 0 then begin
      dominantfreq=fltarr(nx,ny) ; dominant-frequency map
      dominantpower=fltarr(nx,ny) ; dominant-power map (i.e., powers corresponding to dominant frequencies)
  endif
  averagedpower = fltarr(nf)
  for ix=0, nx-1 do begin
      for iy=0, ny-1 do begin
          signal=reform(cube[ix,iy,*])
          IMFcal = walsa_emd_function(signal,stdlimit,dt=dt)
          hhs = walsa_hilbert_spec(IMFcal,dt, marginal=pm, nfilter=nfilter, instfreq=instfreq, amplitudemarginal=amplitudemarginal)
          nimimf = n_elements(IMFcal[0,*])
          if emd then begin
              imf[ix,iy,*,0:nimimf-1] = IMFcal
              instantfreq[ix,iy,*,0:nimimf-1] = instfreq*1000. ; instantaneous frequencies in mHz
          end
          if f0 then powermap[ix,iy,*] = (pm[1:nff-1]*padding)/frequencies[0] else powermap[ix,iy,*] = (pm*padding)/frequencies[0] ; in DN^2/mHz
          amplitude[ix,iy,*] = amplitudemarginal[1:nff-1]*padding

          if nodominantfreq eq 0 then begin
              dominantfreq[ix,iy] = walsa_dominant_frequency(reform(powermap[ix,iy,*]),frequencies,rangefreq,dominantpower=dompm)
              dominantpower[ix,iy] = dompm
          endif

          if nosignificance eq 0 then begin
              Nsig = n_elements(signal)
              ps_perm = fltarr(nf,nperm)
              for ip=0L, nperm-1 do begin
                  permutation = walsa_randperm(Nsig)
                  signalo=reform(originalcube[ix,iy,*]) 
                  y_perm = signalo(permutation)
                  if nodetrendapod eq 0 then $
                     y_perm=walsa_detrend_apod(y_perm,apod,meandetrend,pxdetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original,cadence=cadence,/silent) 
                  IMFcal = walsa_emd_function(y_perm,stdlimit,dt=dt)
                  hhs = walsa_hilbert_spec(IMFcal,dt, marginal=pstmp, nfilter=nfilter)
                  ps_perm[*,ip] = pstmp[1:nff-1]
                  if silent eq 0 then print,string(13b)+' >>> % Running Monte Carlo (significance test): ',(ip*100.)/(nperm-1),format='(a,f4.0,$)'
              endfor
              signif = walsa_confidencelevel(ps_perm, siglevel=siglevel, nf=nf)
              significance[ix,iy,*] = (signif*padding)/frequencies[0] ; in DN^2/mHz
          endif
          averagedpower = averagedpower+reform(powermap[ix,iy,*])
      endfor
      if long(nx) gt 1 then $ 
         writeu,-1,string(format='(%"\r == HHT next row... ",i5,"/",i5)',ix,nx)
  endfor
  powermap = reform(powermap)
  frequencies = reform(frequencies)
  amplitude = reform(amplitude)
  averagedpower = reform(averagedpower/float(nx)/float(ny))
  if emd then begin
      imf = reform(imf)
      instantfreq = reform(instantfreq)
  endif
  if nodominantfreq eq 0 then begin
      dominantfreq = reform(dominantfreq)
      dominantpower = reform(dominantpower)
  endif
  if nosignificance eq 0 then significance = reform(significance)

  return, powermap
end
;--------------------------------- Lomb-Scargle -------------------------------
function getpowerLS,cube,time,OFAC=OFAC,siglevel=siglevel,frequencies=frequencies,significance=significance,nperm=nperm,$
                    nosignificance=nosignificance,averagedpower=averagedpower,amplitude=amplitude,originalcube=originalcube,$
                    dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower,$
                    apod=apod,nodetrendapod=nodetrendapod,pxdetrend=pxdetrend,meandetrend=meandetrend,polyfit=polyfit,$
                    meantemporal=meantemporal,recon=recon,resample_original=resample_original,silent=silent
  ; Lomb-scargle power spectra
  ; The periodogram values (from LNP_TEST) are converted to power (comparable to FFT values) by myltiplying
  ; with 2.*variance(signal,/double)/nt (see Numerical Recipes in C: The Art of Scientific Computing; Press at al. 2007)

  sizecube=size(cube)
  nx=sizecube[1]
  ny=sizecube[2]
  nt=sizecube[3]

  if OFAC gt 1 then begin
      if silent eq 0 then begin
          print, ' '
          print, ' -- Zero Padding (oversampling factor: '+strtrim(OFAC,2)+') .....'
          print, ' '
      endif
  endif

  r = LNP_TEST(reform(time), reform(cube[0,0,*]), /DOUBLE, WK1=frequencies, WK2=pm, OFAC=OFAC)
  frequencies = frequencies*1000. ; in mHz
  nf=n_elements(frequencies)
  powermap=fltarr(nx,ny,nf) ; Lomb-scargle power spectra
  amplitude=fltarr(nx,ny,nf)
  if nosignificance eq 0 then significance=fltarr(nx,ny,nf) ; significancec cube
  if nodominantfreq eq 0 then begin
      dominantfreq=fltarr(nx,ny) ; dominant-frequency map
      dominantpower=fltarr(nx,ny) ; dominant-power map (i.e., powers corresponding to dominant frequencies)
  endif
  averagedpower = fltarr(nf)
  for ix=0, nx-1 do begin
      for iy=0, ny-1 do begin
          signal = reform(cube[ix,iy,*])
          r = LNP_TEST(reform(time), signal, /DOUBLE, WK1=freq, WK2=pm, OFAC=OFAC)
          powermap[ix,iy,*] = ((pm*(2.*variance(signal,/double)/nt))*OFAC)/frequencies[0] ; in DN^2/mHz

          amplitude[ix,iy,*] = sqrt(((2.*pm*(2.*variance(signal,/double)/nt))*OFAC)) ; K. Hocke 1998, Ann. Geophysics, 16, 356

          if nodominantfreq eq 0 then begin
              dominantfreq[ix,iy] = walsa_dominant_frequency(reform(powermap[ix,iy,*]),frequencies,rangefreq,dominantpower=dompm)
              dominantpower[ix,iy] = dompm
          endif

          if nosignificance eq 0 then begin
              Nsig = n_elements(signal)
              ps_perm = fltarr(nf,nperm)
              for ip=0L, nperm-1 do begin
                  permutation = walsa_randperm(Nsig)
                  signalo=reform(originalcube[ix,iy,*]) 
                  y_perm = signalo(permutation)
                  if nodetrendapod eq 0 then $
                     y_perm=walsa_detrend_apod(y_perm,apod,meandetrend,pxdetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original,cadence=cadence,/silent)
                  results = LNP_TEST(time, y_perm, /DOUBLE, WK1=freq, WK2=psp, OFAC=OFAC)
                  ps_perm[*,ip] = psp*(2.*variance(y_perm,/double)/nt)
                  if silent eq 0 then print,string(13b)+' >>> % Running Monte Carlo (significance test): ',(ip*100.)/(nperm-1),format='(a,f4.0,$)'
              endfor
              signif = walsa_confidencelevel(ps_perm, siglevel=siglevel, nf=nf)
              significance[ix,iy,*] = (signif*OFAC)/frequencies[0] ; in DN^2/mHz
          endif
          averagedpower = averagedpower+reform(powermap[ix,iy,*])
      endfor
      if long(nx) gt 1 then $ 
         writeu,-1,string(format='(%"\r == Lomb-Scargle next row... ",i5,"/",i5)',ix,nx)
  endfor
  powermap = reform(powermap)
  frequencies = reform(frequencies)
  amplitude = reform(amplitude)
  averagedpower = reform(averagedpower/float(nx)/float(ny))
  if nodominantfreq eq 0 then begin
      dominantfreq = reform(dominantfreq)
      dominantpower = reform(dominantpower)
  endif
  if nosignificance eq 0 then significance = reform(significance)

  return, powermap
end
;----------------------------------- Welch ---------------------------------
function welch_psd,cube,cadence,frequencies=frequencies,window_size=window_size,overlap=overlap,wfft_size=wfft_size, significance=significance,$
                     nperm=nperm,nosignificance=nosignificance,averagedpower=averagedpower, originalcube=originalcube,siglevel=siglevel, $
                     dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower,apod=apod,silent=silent,$
                     nodetrendapod=nodetrendapod,pxdetrend=pxdetrend,meandetrend=meandetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original

    ; Initialize variables
    sizecube=size(cube)
    nx=sizecube[1]
    ny=sizecube[2]
    nt=sizecube[3]

    ; Assume wfft_size (nfft) is equal to window_size for simplicity and clarity
    wfft_size = window_size

    step_size = window_size-overlap
    num_segments = fix((nt-overlap)/step_size)
    if num_segments le 0 then begin
        PRINT, ' Number of segments: '+strtrim(num_segments, 2)+' (!)'
        PRINT, ' Error: Overlap or window size too large.'
        stop
    endif

    ; Create a Hann window (this should be a changable option in future versions)
    window = (1.0-cos(2 *!pi*findgen(window_size)/(window_size-1)))/2.0

    frequencies = 1./(cadence*2)*findgen(window_size/2+1)/(window_size/2)
    nff=n_elements(frequencies)
    frequencies = frequencies[1:nff-1]
    frequencies = frequencies*1000. ; in mHz
    nf=n_elements(frequencies)
    powermap=fltarr(nx,ny,nf) ; Welch power spectra
    if nosignificance eq 0 then significance=fltarr(nx,ny,nf) ; significance cube
    if nodominantfreq eq 0 then begin
        dominantfreq=fltarr(nx,ny) ; dominant-frequency map
        dominantpower=fltarr(nx,ny) ; dominant-power map (i.e., powers corresponding to dominant frequencies)
    endif
    averagedpower = fltarr(nf)

    for ix=0, nx-1 do begin
        for iy=0, ny-1 do begin
            signal = reform(cube[ix,iy,*])
            ; Process each segment
            psd = FLTARR(window_size / 2 + 1)
            for segment = 0L, num_segments-1 do begin
                start_index = segment * step_size
                end_index = start_index + window_size
                if end_index gt nt then continue
                ; Extract the segment and apply the window
                segment_data = signal[start_index:end_index]*window
                ; Compute the FFT
                segment_fft = FFT(segment_data, wfft_size)
                ; Compute power spectral density
                segment_psd = (ABS(segment_fft))^2/(window_size*wfft_size)
                psd = psd + segment_psd[0:window_size / 2]

            endfor
            ; Normalize the averaged PSD
            psd = psd / num_segments
            powermap[ix,iy,*] = psd[1:nff-1] / frequencies[0] ; in DN^2/mHz

            if nodominantfreq eq 0 then begin
                dominantfreq[ix,iy] = walsa_dominant_frequency(reform(powermap[ix,iy,*]),frequencies,rangefreq,dominantpower=dompm)
                dominantpower[ix,iy] = dompm
            endif

            if nosignificance eq 0 then begin
                Nsig = n_elements(signal)
                ps_perm = fltarr(nf,nperm)
                for ip=0L, nperm-1 do begin
                    permutation = walsa_randperm(Nsig)
                    signalo=reform(originalcube[ix,iy,*]) 
                    y_perm = signalo(permutation)
                    if nodetrendapod eq 0 then $
                       y_perm=walsa_detrend_apod(y_perm,apod,meandetrend,pxdetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original,cadence=cadence,/silent)

                    ; Process each segment
                    psd = FLTARR(window_size / 2 + 1)
                    for segment = 0L, num_segments-1 do begin
                        start_index = segment * step_size
                        end_index = start_index + window_size
                        if end_index gt nt then continue
                        ; Extract the segment and apply the window
                        segment_data = y_perm[start_index:end_index]*window
                        ; Compute the FFT
                        segment_fft = FFT(segment_data, wfft_size)
                        ; Compute power spectral density
                        segment_psd = (ABS(segment_fft))^2/(window_size*wfft_size)
                        psd = psd + segment_psd[0:window_size / 2]

                    endfor
                    ; Normalize the averaged PSD
                    psd = psd / num_segments
                    ps_perm[*,ip] = psd[1:nff-1]
                    if silent eq 0 then print,string(13b)+' >>> % Running Monte Carlo (significance test): ',(ip*100.)/(nperm-1),format='(a,f4.0,$)'
                endfor
                signif = walsa_confidencelevel(ps_perm, siglevel=siglevel, nf=nf)
                significance[ix,iy,*] = signif/frequencies[0] ; in DN^2/mHz
            endif
            averagedpower = averagedpower+reform(powermap[ix,iy,*])
        endfor
        if long(nx) gt 1 then $ 
           writeu,-1,string(format='(%"\r == FFT next row... ",i5,"/",i5)',ix,nx)
    endfor

    powermap = reform(powermap)
    frequencies = reform(frequencies)
    averagedpower = reform(averagedpower/float(nx)/float(ny))
    if nodominantfreq eq 0 then begin
        dominantfreq = reform(dominantfreq)
        dominantpower = reform(dominantpower)
    endif
    if nosignificance eq 0 then significance = reform(significance)

    return, powermap
end
;----------------------------------- FOURIER ----------------------------------
function getpowerFFT,cube,cadence,siglevel=siglevel,padding=padding,frequencies=frequencies,significance=significance,$
                     nperm=nperm,nosignificance=nosignificance,averagedpower=averagedpower,amplitude=amplitude,originalcube=originalcube,$
                     dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower,apod=apod,silent=silent,$
                     nodetrendapod=nodetrendapod,pxdetrend=pxdetrend,meandetrend=meandetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original
  ; Fast Fourier Transform (FFT) power spectra

  if padding gt 1 then begin ; zero padding (optional): to increase frequency resolution
      if silent eq 0 then begin
          print, ' '
          print, ' -- Zero Padding (oversampling factor: '+strtrim(padding,2)+') .....'
          print, ' '
      endif
      nx = n_elements(cube[*,0,0])
      ny = n_elements(cube[0,*,0])
      nt = n_elements(cube[0,0,*])        
      padded=fltarr(nx,ny,padding*nt)
      mid_point = ROUND((padding*nt)/2.)
      lower_point = mid_point-nt/2.
      upper_point = mid_point+nt/2.-1
      padded[*,*,mid_point-nt/2.:mid_point+nt/2.-1] = cube
      cube = padded
  endif
  sizecube=size(cube)
  nx=sizecube[1]
  ny=sizecube[2]
  nt=sizecube[3]

  frequencies = 1./(cadence*2)*findgen(nt/2+1)/(nt/2)
  nff=n_elements(frequencies)
  frequencies = frequencies[1:nff-1]
  frequencies = frequencies*1000. ; in mHz
  nf=n_elements(frequencies)
  powermap=fltarr(nx,ny,nf) ; FFT power spectra
  amplitude=complexarr(nx,ny,nf) ; FFT amplitudes
  if nosignificance eq 0 then significance=fltarr(nx,ny,nf) ; significance cube
  if nodominantfreq eq 0 then begin
      dominantfreq=fltarr(nx,ny) ; dominant-frequency map
      dominantpower=fltarr(nx,ny) ; dominant-power map (i.e., powers corresponding to dominant frequencies)
  endif
  averagedpower = fltarr(nf)
  for ix=0, nx-1 do begin
      for iy=0, ny-1 do begin
          signal = reform(cube[ix,iy,*])
          ; single-sided power is doubled (compared to double-sided power), assuming P(−f) = P(f):
          spec = (fft(signal,-1,/double))[0:nt/2.]
          pm = 2.*(ABS(spec)^2)
          powermap[ix,iy,*] = (pm[1:nff-1]*padding)/frequencies[0] ; in DN^2/mHz
          amplitude[ix,iy,*] = spec[1:nff-1]*padding

          if nodominantfreq eq 0 then begin
              dominantfreq[ix,iy] = walsa_dominant_frequency(reform(powermap[ix,iy,*]),frequencies,rangefreq,dominantpower=dompm)
              dominantpower[ix,iy] = dompm
          endif

          if nosignificance eq 0 then begin
              Nsig = n_elements(signal)
              ps_perm = fltarr(nf,nperm)
              for ip=0L, nperm-1 do begin
                  permutation = walsa_randperm(Nsig)
                  signalo=reform(originalcube[ix,iy,*])
                  y_perm = signalo(permutation)
                  if nodetrendapod eq 0 then $
                     y_perm=walsa_detrend_apod(y_perm,apod,meandetrend,pxdetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original,cadence=cadence,/silent)
                  pstmp = 2.*(ABS((fft(y_perm,-1,/double))[0:nt/2.])^2.) 
                  ps_perm[*,ip] = pstmp[1:nff-1]
                  if silent eq 0 then print,string(13b)+' >>> % Running Monte Carlo (significance test): ',(ip*100.)/(nperm-1),format='(a,f4.0,$)'
              endfor
              signif = walsa_confidencelevel(ps_perm, siglevel=siglevel, nf=nf)
              significance[ix,iy,*] = (signif*padding)/frequencies[0] ; in DN^2/mHz
          endif
          averagedpower = averagedpower+reform(powermap[ix,iy,*])
      endfor
      if long(nx) gt 1 then $ 
         writeu,-1,string(format='(%"\r == FFT next row... ",i5,"/",i5)',ix,nx)
  endfor
  powermap = reform(powermap)
  frequencies = reform(frequencies)
  averagedpower = reform(averagedpower/float(nx)/float(ny))
  amplitude = reform(amplitude)
  if nodominantfreq eq 0 then begin
      dominantfreq = reform(dominantfreq)
      dominantpower = reform(dominantpower)
  endif
  if nosignificance eq 0 then significance = reform(significance)

  return, powermap
end
;------------------------------------ WAVELET ---------------------------------
function getpowerWAVELET,cube,cadence,dj=dj,mother=mother,siglevel=siglevel,global=global,frequencies=frequencies,$
                         significance=significance,coicube=coicube,oglobal=oglobal,colornoise=colornoise,param=param,$
                         padding=padding,nosignificance=nosignificance,averagedpower=averagedpower,psd=psd,$
                         dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower,$
                         originalcube=originalcube,apod=apod,nodetrendapod=nodetrendapod,pxdetrend=pxdetrend,meandetrend=meandetrend,$
                         polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original, nperm=nperm, rgws=rgws, silent=silent
  ; Wavelet power spectra: either wavelet spectra, or global wavelet spectra (traditional or improved versions)

  input_data = cube ; to prevent the input data be modified by functions unintentionally

  if padding gt 1 then begin ; zero padding (optional): to increase frequency resolution
      if silent eq 0 then begin
          print, ' '
          print, ' -- Zero Padding (oversampling factor: '+strtrim(padding,2)+') .....'
          print, ' '
      endif
      nx = n_elements(input_data[*,0,0])
      ny = n_elements(input_data[0,*,0])
      nt = n_elements(input_data[0,0,*])        
      padded=fltarr(nx,ny,padding*nt)
      mid_point = ROUND((padding*nt)/2.)
      lower_point = mid_point-nt/2.
      upper_point = mid_point+nt/2.-1
      padded[*,*,mid_point-nt/2.:mid_point+nt/2.-1] = input_data
      input_data = padded
  endif

  siglevel=1.0-siglevel ; different convention

  sizecube=size(input_data)
  nx=sizecube[1]
  ny=sizecube[2]
  nt=sizecube[3]
  col=reform(input_data[0,0,*])
  col = (col - TOTAL(col)/nt)
  ; lag1 = (A_CORRELATE(col,1) + SQRT(A_CORRELATE(col,2)))/2.
  ; Wavelet transform:
  wave = walsa_wavelet(reform(col),cadence,PERIOD=period,PAD=1,COI=coi,MOTHER=mother,/RECON,dj=dj,param=param,J=J,/nodetrendapod)

  frequencies = 1./period
  frequencies = frequencies*1000. ; in mHz

  nf = n_elements(frequencies)

  if (global+oglobal+rgws) gt 1 then begin
      print
      print, ' --- [!] Only one of the /global, /oglobal, or /rgws can be flagged at a time!'
      print
      stop
  endif

  if silent eq 0 then print, ' '
  if global eq 1 then begin
      if silent eq 0 then print, ' ...... global: output "Traditional" Global Wavelet Spectrum'
      ftcube=fltarr(nx,ny,nf)
  endif
  if oglobal eq 1 then begin
      if silent eq 0 then print, ' ...... oglobal: output Global Wavelet Spectrum excluding CoI regions'
      ftcube = fltarr(nx,ny,nf) ; power
  endif
  if rgws eq 1 then begin
      if silent eq 0 then print, ' ...... rgws: output rgws Wavelet Spectrum '
      if silent eq 0 then print, ' ...... (power-weighted significant frequency distribution, unaffected by CoI)'
      ftcube = fltarr(nx,ny,nf) ; power
  endif
  if global eq 0 and oglobal eq 0 and rgws eq 0 then begin
      if silent eq 0 then print, ' ...... output Wavelet Spectra '
      ftcube=fltarr(nx,ny,nt,nf)
  endif

  if silent eq 0 then begin
      print, ' '
      print, ' Wavelet (mother) function: '+mother
      print, ' dj: '+strtrim(dj,2)
      print, ' '
  endif

  iloop = 0
  numm = nx*float(ny)

  if nosignificance eq 0 then significance=fltarr(nx,ny,nf)
  if nodominantfreq eq 0 then begin
      dominantfreq=fltarr(nx,ny) ; dominant-frequency map
      dominantpower=fltarr(nx,ny) ; dominant-power map (i.e., powers corresponding to dominant frequencies)
  endif
  coicube = fltarr(nx,ny,nt)
  averagedpower = fltarr(nf)
  for ix=0, nx-1 do begin
      for iy=0, ny-1 do begin
        col=reform(input_data[ix,iy,*])
        col = (col - TOTAL(col)/nt)
        col_copy = col
        ; lag1 = (A_CORRELATE(col,1) + SQRT(A_CORRELATE(col,2)))/2.
        ; Wavelet transform:
        wave = walsa_wavelet(col_copy,cadence,PERIOD=period,PAD=1,COI=coi,MOTHER=mother,param=param,/RECON,dj=dj,scale=scale,$
                                                    SIGNIF=SIGNIF,SIGLVL=siglevel,/nodetrendapod,colornoise=colornoise,power=power)

        if global eq 1 then begin
            global_ws = TOTAL(power,1,/nan)/nt  ; global wavelet spectrum (GWS)
            ftcube[ix,iy,*] = (reform(global_ws)*padding);/frequencies[nf-1]

            global_amp = TOTAL(wave,1,/nan)/nt

            if nodominantfreq eq 0 then begin
                dominantfreq[ix,iy] = walsa_dominant_frequency(reform(ftcube[ix,iy,*]),frequencies,rangefreq,dominantpower=dompm)
                dominantpower[ix,iy] = dompm
            endif

            ; GWS significance levels:
            if nosignificance eq 0 then begin
                dof = nt - scale   ; the -scale corrects for padding at edges
                global_signif = walsa_wave_signif(col,cadence,scale,1, LAG1=0.0,DOF=dof,MOTHER=mother,CDELTA=Cdelta,PSI0=psi0,siglvl=siglevel)
                significance[ix,iy,*] = (reform(global_signif)*padding);/frequencies[nf-1]
            endif
            averagedpower = averagedpower+reform(ftcube[ix,iy,*])
        endif

        if global eq 0 and oglobal eq 0 and rgws eq 0 then begin
            ftcube[ix,iy,*,*] = (reform(power)*padding); /frequencies[nf-1]
            if nosignificance eq 0 then significance[ix,iy,*] = (reform(SIGNIF)*padding); /frequencies[nf-1] ; in DN^2/mHz
            coicube[ix,iy,*] = reform(coi)
        endif

        ; oglobal: time-average wavelet power only over the areas not affected by CoI
        if oglobal eq 1 then begin
            opower = fltarr(nt,nf)+!VALUES.F_NAN
            for i=0L, nt-1 do begin
                pcol = reform(power[i,*])
                ii = where(reform(period) lt coi[i], pnum)
                if pnum gt 0 then opower(i,ii) = pcol(ii)
            endfor
            opower = mean(opower,dimension=1,/nan)
            ftcube[ix,iy,*] = (opower*padding); /frequencies[nf-1] ; in DN^2

            if nodominantfreq eq 0 then begin
                dominantfreq[ix,iy] = walsa_dominant_frequency(reform(ftcube[ix,iy,*]),frequencies,rangefreq,dominantpower=dompm)
                dominantpower[ix,iy] = dompm
            endif

            ; GWS significance levels:
            if nosignificance eq 0 then begin
                ; dof = nt - scale   ; the -scale corrects for padding at edges
                ; global_signif = walsa_wave_signif(col,cadence,scale,1, LAG1=0.0,DOF=dof,MOTHER=mother,CDELTA=Cdelta,PSI0=psi0,siglvl=siglevel)
                ; significance[ix,iy,*] = (reform(global_signif)*padding); /frequencies[nf-1] ; in DN^2
                Nsig = n_elements(col)
                ps_perm = fltarr(nf,nperm)
                for ip=0L, nperm-1 do begin
                    permutation = walsa_randperm(Nsig)
                    signalo=reform(originalcube[ix,iy,*]) 
                    y_perm = signalo(permutation)
                    if nodetrendapod eq 0 then $
                       y_perm=walsa_detrend_apod(y_perm,apod,meandetrend,pxdetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original,cadence=cadence,/silent)

                    y_perm = (y_perm - TOTAL(y_perm)/nt)
                    wave = walsa_wavelet(y_perm,cadence,PERIOD=period,PAD=1,COI=coi,MOTHER=mother,param=param,/RECON,dj=dj,scale=scale,$
                                         /nodetrendapod,colornoise=colornoise,power=powertemp)

                    opower = fltarr(nt,nf)+!VALUES.F_NAN
                    for i=0L, nt-1 do begin
                        pcol = reform(powertemp[i,*])
                        ii = where(reform(period) lt coi[i], pnum)
                        if pnum gt 0 then opower(i,ii) = pcol(ii)
                    endfor
                    ps_perm = mean(opower,dimension=1,/nan) ; time average only over the areas not affected by CoI
                    if silent eq 0 then print,string(13b)+' >>> % Running Monte Carlo (significance): ',(ip*100.)/(nperm-1),format='(a,f4.0,$)'
                endfor
                signif = walsa_confidencelevel(ps_perm, siglevel=siglevel, nf=nf)
                significance[ix,iy,*] = (signif*padding); /frequencies[nf-1] ; in DN^2
            endif
            averagedpower = averagedpower+reform(ftcube[ix,iy,*])
        endif

        ; rgws: time-integral of wavelet power only over the areas not affected by CoI and only for those above the significance level.
        ; i.e., 'distribution' of significant frequencies (unaffected by CoI) weighted by power.
        if rgws eq 1 then begin
            isig = REBIN(TRANSPOSE(signif),nt,nf)
            istest = where(power/isig lt 1.0, numtest)
            if numtest gt 0 then power[istest]=!VALUES.F_NAN
            ipower = fltarr(nt,nf)+!VALUES.F_NAN
            for i=0L, nt-1 do begin
                pcol = reform(power[i,*])
                ii = where(reform(period) lt coi[i], pnum)
                if pnum gt 0 then ipower[i,ii] = pcol[ii]
            endfor
            ;ipower = mean(ipower,dimension=1,/nan)
            ipower = total(ipower,1,/nan)

            ftcube[ix,iy,*] = (ipower*padding); /frequencies[nf-1] ; in DN^2

            if nodominantfreq eq 0 then begin
                dominantfreq[ix,iy] = walsa_dominant_frequency(reform(ftcube[ix,iy,*]),frequencies,rangefreq,dominantpower=dompm)
                dominantpower[ix,iy] = dompm
            endif
            averagedpower = averagedpower+reform(ftcube[ix,iy,*])
        endif

      endfor
      iloop = long(iloop+1.)
      if silent eq 0 then if nx gt 1 or ny gt 1 then print,string(13b)+' >>> % finished: ',(iloop*100.)/nx,format='(a,f4.0,$)'
  endfor

  powermap = reform(ftcube)
  frequencies = reform(frequencies)
  averagedpower = reform(averagedpower/float(nx)/float(ny))
  if nodominantfreq eq 0 then begin
      dominantfreq = reform(dominantfreq)
      dominantpower = reform(dominantpower)
  endif
  if nosignificance eq 0 then significance = reform(significance)
  coicube = reform(coicube)

  ; if (global+oglobal+rgws) gt 0 AND psd eq 1 then begin
  ;       ; Interpolate the power spectrum to a uniform frequency array (Wavelet's frequency resolution changes with frequency)
  ;       uniform_freqs = findgen(n_elements(frequencies)) * (max(frequencies) - min(frequencies)) / (n_elements(frequencies) - 1) + min(frequencies)
  ;
  ;       powermap = interpol(powermap, frequencies, uniform_freqs, /SPLINE)
  ;
  ;       frequencies = uniform_freqs
  ;       delta_freq = ABS(frequencies[1] - frequencies[0])
  ;
  ;       powermap = powermap / delta_freq
  ;
  ;       averagedpower = averagedpower / delta_freq
  ;       if nodominantfreq eq 0 then dominantpower = dominantpower / delta_freq
  ;       significance = significance / delta_freq
  ; endif

  return, powermap
end
;==================================================== MAIN ROUTINE ====================================================
function walsa_speclizer,data,time,$ ; main inputs
                        frequencies=frequencies, significance=significance, coicube=coicube, imf=imf, instantfreq=instantfreq,$ ; main (additional) outputs
                        averagedpower=averagedpower,amplitude=amplitude, period=period, $
                        fft=fft, lombscargle=lombscargle, wavelet=wavelet, hht=hht, welch=welch,$ ; type of analysis
                        padding=padding, apod=apod, nodetrendapod=nodetrendapod, pxdetrend=pxdetrend, meandetrend=meandetrend,$ ; padding and apodization parameters
                        polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original, psd=psd,$
                        siglevel=siglevel, nperm=nperm, nosignificance=nosignificance,$ ; significance-level parameters
                        mother=mother, param=param, dj=dj, global=global, oglobal=oglobal, rgws=rgws, colornoise=colornoise,$ ; Wavelet parameters/options
                        stdlimit=stdlimit, nfilter=nfilter, emd=emd,$ ; HHT parameters/options
                        mode=mode, silent=silent, $
                        dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower, $ ; dominant frequency
                        window_size=window_size, overlap=overlap, wfft_size=wfft_size ; Welch parameters

cube = reform(data)
sizecube = size(cube)
if sizecube[0] ne 3 then begin
    if sizecube[0] eq 1 then begin
        blablacube = fltarr(1,1,sizecube[1])
        blablacube[0,0,*] = cube
        cube = blablacube
    endif else begin
        print, ' '
        print, ' [!] The datacube must have either 1 or 3 dimension(s).'
        print, ' '
        stop
    endelse
endif

ii = where(~finite(cube),/null,cnull)
if cnull gt 0 then cube(ii) = median(cube)

cadence=walsa_mode(walsa_diff(time))
temporal_Nyquist = 1. / (cadence * 2.)

if n_elements(silent) eq 0 then silent=0

if silent eq 0 then begin
    print,' '
    if sizecube[0] eq 1 then print,'The input datacube is of size: ['+ARR2STR(sizecube[1],/trim)+']'
    if sizecube[0] eq 3 then $
    print,'The input datacube is of size: ['+ARR2STR(sizecube[1],/trim)+', '+ARR2STR(sizecube[2],/trim)+', '+ARR2STR(sizecube[3],/trim)+']'
    print,' '
    print,'Temporally, the important values are:'
    print,'    2-element duration (Nyquist period) = '+ARR2STR((cadence * 2.),/trim)+' seconds'
    if sizecube[0] eq 1 then print,'    Time series duration = '+ARR2STR(cadence*sizecube[1],/trim)+' seconds'
    if sizecube[0] eq 3 then print,'    Time series duration = '+ARR2STR(cadence*sizecube[3],/trim)+' seconds'
    print,'    Nyquist frequency = '+ARR2STR(temporal_Nyquist*1000.,/trim)+' mHz'
    print, ' '
endif

if n_elements(global) eq 0 then global=0 ; if set, global wavelet will be returned; otherwsie wavelet 2D power spectra (default)
if n_elements(dj) eq 0 then dj=0.025

if n_elements(fft) eq 0 then fft=0
if n_elements(psd) eq 0 then psd=0
if n_elements(lombscargle) eq 0 then lombscargle=0
if n_elements(hht) eq 0 then hht=0
if n_elements(welch) eq 0 then welch=0
if n_elements(wavelet) eq 0 then wavelet=0
if n_elements(nosignificance) eq 0 then nosignificance=0
if n_elements(nodominantfreq) eq 0 then nodominantfreq=0

if n_elements(siglevel) eq 0 then siglevel=0.05 ; 5% significance level = 95% confidence level
if n_elements(nperm) eq 0 then nperm=1000
if n_elements(oglobal) eq 0 then oglobal=0
if n_elements(rgws) eq 0 then rgws=0

if n_elements(padding) eq 0 then padding=1
if n_elements(stdlimit) eq 0 then stdlimit=0.2
if n_elements(nfilter) eq 0 then nfilter=3
if n_elements(emd) eq 0 then emd=0
if n_elements(colornoise) eq 0 then colornoise=0 ; if set, noise background based on Auchère et al. 2017, ApJ, 838, 166 / 2016, ApJ, 825, 110

if n_elements(mode) eq 0 then mode=1

if n_elements(apod) eq 0 then apod=0.1 ; width of Tukey window (for apodization)
if n_elements(nodetrendapod) eq 0 then nodetrendapod=0 ; detrend the signal and apodize it
if n_elements(pxdetrend) eq 0 then pxdetrend=2 ; do proper linear detrending
if n_elements(meandetrend) eq 0 then meandetrend=0 ; no spatial detrending
if n_elements(mother) eq 0 then mother='Morlet' ; possible functions: 'Morlet', 'DOG', 'Paul'

if n_elements(NFFT) eq 0 then NFFT=256

; detrend and apodize the cube
if nodetrendapod eq 0 then begin
    apocube=walsa_detrend_apod(cube,apod,meandetrend,pxdetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,cadence=cadence,resample_original=resample_original,silent=silent) 
endif else apocube=cube

sizecube = size(apocube)
if sizecube[0] ne 3 then begin
    if sizecube[0] eq 1 then begin
        blablacube = fltarr(1,1,sizecube[1])
        blablacube[0,0,*] = apocube
        apocube = blablacube
    endif
endif

if fft then begin
    if silent eq 0 then begin
        print, ' '
        print, ' -- Perform FFT (Fast Fourier Transform) .....'
        print, ' '
    endif
    power=getpowerFFT(apocube,cadence,siglevel=siglevel,padding=padding,frequencies=frequencies,significance=significance,averagedpower=averagedpower,$
        nperm=nperm,nosignificance=nosignificance,dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower,$
        amplitude=amplitude,originalcube=cube,apod=apod,nodetrendapod=nodetrendapod,pxdetrend=pxdetrend,meandetrend=meandetrend,polyfit=polyfit,$
        meantemporal=meantemporal,recon=recon,resample_original=resample_original,silent=silent)
endif

if lombscargle then begin
    if silent eq 0 then begin
        print, ' '
        print, ' -- Perform Lomb-Scargle Transform .....'
        print, ' '
    endif
    power=getpowerLS(apocube,time,OFAC=padding,siglevel=siglevel,frequencies=frequencies,significance=significance,averagedpower=averagedpower,amplitude=amplitude,$
        nperm=nperm,nosignificance=nosignificance,dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower,originalcube=cube,$
        apod=apod,nodetrendapod=nodetrendapod,pxdetrend=pxdetrend,meandetrend=meandetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original,silent=silent)
endif

if welch then begin
    if silent eq 0 then begin
        print, ' '
        print, ' -- Perform Welch method .....'
        print, ' '
    endif
    power=welch_psd(apocube, cadence, frequencies=frequencies, window_size=window_size, overlap=overlap, wfft_size=wfft_size, significance=significance,$
                     nperm=nperm,nosignificance=nosignificance,averagedpower=averagedpower, originalcube=cube,siglevel=siglevel,$
                     dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower,apod=apod,silent=silent,$
                     nodetrendapod=nodetrendapod,pxdetrend=pxdetrend,meandetrend=meandetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original)
endif


if wavelet then begin
    if silent eq 0 then begin
        print, ' '
        print, ' -- Perform Wavelet Transform .....'
        print, ' '
    endif
    power=getpowerWAVELET(apocube,cadence,dj=dj,mother=mother,siglevel=siglevel,global=global,frequencies=frequencies,averagedpower=averagedpower,rgws=rgws,$
        significance=significance,coicube=coicube,oglobal=oglobal,colornoise=colornoise,padding=padding,nosignificance=nosignificance,originalcube=cube,$
        dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower,param=param,nperm=nperm,psd=psd,$
        apod=apod,nodetrendapod=nodetrendapod,pxdetrend=pxdetrend,meandetrend=meandetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original,silent=silent)
endif

if hht then begin
    if silent eq 0 then begin
        print, ' '
        print, ' -- Perform HHT (Hilbert-Huang Transform) .....'
        print, ' '
    endif
    power=getpowerHHT(apocube,cadence,stdlimit,nfilter=nfilter,significance=significance,siglevel=siglevel,nperm=nperm,originalcube=cube,$
        padding=padding,frequencies=frequencies,nosignificance=nosignificance,emd=emd,imf=imf,instantfreq=instantfreq,amplitude=amplitude,$
        dominantfreq=dominantfreq,rangefreq=rangefreq,nodominantfreq=nodominantfreq,dominantpower=dominantpower,averagedpower=averagedpower,$
        apod=apod,nodetrendapod=nodetrendapod,pxdetrend=pxdetrend,meandetrend=meandetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,$
        resample_original=resample_original,silent=silent)
endif

period = 1000./frequencies

sizecube=size(power)
dim=sizecube[0]
nx=sizecube[1]
ny=sizecube[2]

if dim le 2 then begin ; wavelet power spectrum
    if (mode eq 0) then begin
        power=alog10(power)
        if nosignificance eq 0 then significance=alog10(significance)
    endif
    if (mode eq 2) then begin
        power=sqrt(power)
        if nosignificance eq 0 then significance=sqrt(significance)
    endif
endif

if dim eq 3 then begin ; power spectra at multiple pixels
    nf=sizecube[3]
    if (mode eq 0) then begin 
        for jnn=0L, nf-1 do power[*,*,jnn]=alog10(power[*,*,jnn])
        if nosignificance eq 0 then for jnn=0L, nf-1 do significance[*,*,jnn]=alog10(significance[*,*,jnn])
    endif
    if (mode eq 2) then begin
        for jnn=0L, nf-1 do power[*,*,jnn]=sqrt(power[*,*,jnn])
        if nosignificance eq 0 then for jnn=0L, nf-1 do significance[*,*,jnn]=sqrt(significance[*,*,jnn])
    endif
endif

if dim eq 4 then begin ; wavelet power spectra at multiple pixels
    if (mode eq 0) then begin
        for jx=0L, nx-1 do for jy=0L, ny-1 do power[jx,jy,*,*]=alog10(power[jx,jy,*,*])
        if nosignificance eq 0 then for jx=0L, nx-1 do for jy=0L, ny-1 do significance[jx,jy,*,*]=alog10(significance[jx,jy,*,*])
    endif
    if (mode eq 2) then begin
        for jx=0L, nx-1 do for jy=0L, ny-1 do power[jx,jy,*,*]=sqrt(power[jx,jy,*,*])
        if nosignificance eq 0 then for jx=0L, nx-1 do for jy=0L, ny-1 do significance[jx,jy,*,*]=sqrt(significance[jx,jy,*,*])
    endif
endif

if silent eq 0 then begin
    PRINT
    if (mode eq 0) then print,' mode = 0: log(power)'
    if (mode eq 1) then print,' mode = 1: linear power' 
    if (mode eq 2) then print,' mode = 2: sqrt(power)'

    print, ''
    print, 'COMPLETED!'
    print,''
endif

return, power
end

k-ω Analysis and Fourier Filtering

WaLSA_QUB_QUEEFF

A variant of the QUEEns Fourier Filtering (QUEEFF) code, to compute k-ω diagram and perform Fourier filtering in the k-ω space.

WaLSA_qub_queeff.pro

; -----------------------------------------------------------------------------------------------------
; WaLSAtools: Wave analysis tools
; Copyright (C) 2025 WaLSA Team - Shahin Jafarzadeh et al.
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
; 
; http://www.apache.org/licenses/LICENSE-2.0
; 
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
; 
; Note: If you use WaLSAtools for research, please consider citing:
; Jafarzadeh, S., Jess, D. B., Stangalini, M. et al. 2025, Nature Reviews Methods Primers, in press.
; -----------------------------------------------------------------------------------------------------


;+
; NAME: WaLSA_QUB_QUEEFF
;       part of -- WaLSAtools --
;
; ORIGINAL CODE: QUEEns Fourier Filtering (QUEEFF) code
; WRITTEN, ANNOTATED, TESTED AND UPDATED BY:
; (1) Dr. David B. Jess
; (2) Dr. Samuel D. T. Grant
; The original code along with its manual can be downloaded at: https://bit.ly/37mx9ic
;
; WaLSA_QUB_QUEEFF: Slightly modified (i.e., a few additional keywords added) by Shahin Jafarzadeh
;
; CHECK DEPENDENCIES (MAKE SURE ALL REQUIRED PROGRAMMES ARE INSTALLED):
; NOTE
; @/Users/dbj/ARC/IDL_programmes/Fourier_filtering/QUEEFF_code/QUEEFF_dependencies.bat
;
; CALLING SEQUENCE:
;   EXAMPLES:
;   walsa_qub_queeff, datacube, arcsecpx, time=time, power=power, wavenumber=wavenumber, frequencies=frequencies, koclt=1
;   walsa_qub_queeff, datacube, arcsecpx, cadence, /filtering, power=power, wavenumber=wavenumber, frequencies=frequencies, filtered_cube=filtered_cube
;
; + INPUTS:
;   datacube   input datacube, normally in the form of [x, y, t] 
;              [note - at present the input datacube needs to have identical x and y dimensions. if not supplied like this the datacube will be cropped accordingly!]
;   cadence    delta time between sucessive frames - given in seconds. if not set, time must be provided (see optional inputs)
;   arcsecpx   spatial sampling of the input datacube - given in arcseconds per pixel
;
; + OPTIONAL INPUTS:
; (if optional inputs not supplied, the user will need to interact with the displayed k-omega diagram to define these values)
;   time             observing times in seconds (1d array). it is ignored if cadence is provided
;   filtering        if set, filterring is proceeded
;   f1               optional lower (temporal) frequency to filter - given in mhz
;   f2               optional upper (temporal) frequency to filter - given in mhz
;   k1               optional lower (spatial) wavenumber to filter - given in arcsec^-1 (where k = (2*!pi)/wavelength)
;   k2               optional upper (spatial) wavenumber to filter - given in arcsec^-1 (where k = (2*!pi)/wavelength)
;   spatial_torus    makes the annulus used for spatial filtering have a gaussian-shaped profile (useful for preventing aliasing). default: 1
;                    if equal to 0, it is not applied.
;   temporal_torus   makes the temporal filter have a gaussian-shaped profile (useful for preventing aliasing). default: 1
;                    if equal to 0, it is not applied.
;   no_spatial_filt  optional keyword that ensures no spatial filtering is performed on the dataset (i.e., only temporal filtering)
;   no_temporal_filt optional keyword that ensures no temporal filtering is performed on the dataset (i.e., only spatial filtering)
;   silent:          if set, the k-ω diagram is not plotted
;   clt:             color table number (idl ctload)
;   koclt:           custom color tables for k-ω diagram (currently available: 1 and 2)
;   threemin:        if set, a horizontal line marks the three-minute periodicity
;   fivemin:         if set, a horizontal line marks the five-minute periodicity
;   xlog:            if set, x-axis (wavenumber) is plotted in logarithmic scale (base 10)
;   ylog:            if set, y-axis (frequency) is plotted in logarithmic scale (base 10)
;   xrange:          x-axis (wavenumber) range
;   yrange:          y-axis (frequency) range
;   nox2:            if set, 2nd x-axis (spatial size, in arcsec) is not plotted
;                    (spatial size (i.e., wavelength) = (2*!pi)/wavenumber)
;   noy2:            if set, 2nd y-axis (period, in sec) is not plotted
;                    (p = 1000/frequency)
;   smooth:          if set, power is smoothed
;   epsfilename:     if provided (as a string), an eps file of the k-ω diagram is made
;   mode:            outputted power mode: 0 = log(power) (default), 1 = linear power, 2 = sqrt(power) = amplitude
;
; + OUTPUTS:
;   power:           2d array of power (see mode for the scale)
;                    (in dn^2/mhz, i.e., normalized to frequency resolution)
;   frequencies:     1d array of frequencies (in mhz)
;   wavenumber:      1d array of wavenumber (in arcsec^-1)
;   filtered_cube:   3d array of filtered datacube (if filtering is set)
;
;
; IF YOU USE THIS CODE, THEN PLEASE CITE THE ORIGINAL PUBLICATION WHERE IT WAS USED:
; Jess et al. 2017, ApJ, 842, 59 (http://adsabs.harvard.edu/abs/2017ApJ...842...59J)
;-

pro walsa_qub_queeff, datacube, arcsecpx, cadence=cadence, time=time, $ ; main inputs
                      power=power, wavenumber=wavenumber, frequencies=frequencies, filtered_cube=filtered_cube, $ ; main (additional) outputs
                      filtering=filtering, f1=f1, f2=f2, k1=k1, k2=k2, spatial_torus=spatial_torus, temporal_torus=temporal_torus, $ ; filtering options
                      no_spatial_filt=no_spatial_filt, no_temporal_filt=no_temporal_filt, $
                      clt=clt, koclt=koclt, threemin=threemin, fivemin=fivemin, xlog=xlog, ylog=ylog, xrange=xrange, yrange=yrange, $ ; plotting keywords
                      xtitle=xtitle, ytitle=ytitle, x2ndaxistitle=x2ndaxistitle, y2ndaxistitle=y2ndaxistitle, $
                      epsfilename=epsfilename, noy2=noy2, nox2=nox2, smooth=smooth, silent=silent, mode=mode

if n_elements(xtitle) eq 0 then xtitle='Wavenumber (arcsec!U-1!N)'
if n_elements(ytitle) eq 0 then ytitle='Frequency (mHz)'
if n_elements(x2ndaxistitle) eq 0 then x2ndaxistitle='Spatial size (arcsec)!C'
if n_elements(y2ndaxistitle) eq 0 then y2ndaxistitle='Period (s)'
if n_elements(cadence) eq 0 then cadence=walsa_mode(walsa_diff(time))
; DEFINE THE SCREEN RESOLUTION TO ENSURE THE PLOTS DO NOT SPILL OVER THE EDGES OF THE SCREEN
dimensions = GET_SCREEN_SIZE(RESOLUTION=resolution)
xscreensize = dimensions[0]
yscreensize = dimensions[1]
IF (xscreensize le yscreensize) THEN smallest_screensize = xscreensize
IF (yscreensize le xscreensize) THEN smallest_screensize = yscreensize

xsize_cube = N_ELEMENTS(datacube[*,0,0])
ysize_cube = N_ELEMENTS(datacube[0,*,0])
zsize_cube = N_ELEMENTS(datacube[0,0,*])

; FORCE THE CUBES TO HAVE THE SAME SPATIAL DIMENSIONS
IF xsize_cube gt ysize_cube THEN datacube = TEMPORARY(datacube[0:(ysize_cube-1), *, *])
IF xsize_cube gt ysize_cube THEN xsize_cube = ysize_cube
IF ysize_cube gt xsize_cube THEN datacube = TEMPORARY(datacube[*, 0:(xsize_cube-1), *])
IF ysize_cube gt xsize_cube THEN ysize_cube = xsize_cube

if n_elements(spatial_torus) eq 0 then spatial_torus = 1
if n_elements(temporal_torus) eq 0 then temporal_torus = 1

if n_elements(xlog) eq 0 then xlog = 0
if n_elements(ylog) eq 0 then ylog = 0
if n_elements(nox2) eq 0 then nox2 = 0
if n_elements(noy2) eq 0 then noy2 = 0
if not keyword_set(mode) then mode=0
if n_elements(epsfilename) eq 0 then eps = 0 else eps = 1

if n_elements(silent) eq 0 then silent = 0
if n_elements(filtering) eq 0 then filtering = 0 else silent = 0

; CALCULATE THE NYQUIST FREQUENCIES
spatial_Nyquist  = (2.*!pi) / (arcsecpx * 2.)
temporal_Nyquist = 1. / (cadence * 2.)

print,''
print,'The input datacube is of size: ['+strtrim(xsize_cube,2)+', '+strtrim(ysize_cube,2)+', '+strtrim(zsize_cube,2)+']'
print,''
print,'Spatially, the important values are:'
print,'    2-pixel size = '+strtrim((arcsecpx * 2.),2)+' arcsec'
print,'    Field of view size = '+strtrim((arcsecpx * xsize_cube),2)+' arcsec'
print,'    Nyquist wavenumber = '+strtrim(spatial_Nyquist,2)+' arcsec^-1'
IF KEYWORD_SET(no_spatial_filt) THEN print, '***NO SPATIAL FILTERING WILL BE PERFORMED***'
print,''
print,'Temporally, the important values are:'
print,'    2-element duration (Nyquist period) = '+strtrim((cadence * 2.),2)+' seconds'
print,'    Time series duration = '+strtrim(cadence*zsize_cube,2)+' seconds'
print,'    Nyquist frequency = '+strtrim(temporal_Nyquist*1000.,2)+' mHz'
IF KEYWORD_SET(no_temporal_filt) THEN print, '***NO TEMPORAL FILTERING WILL BE PERFORMED***'

; MAKE A k-omega DIAGRAM
sp_out = DBLARR(xsize_cube/2,zsize_cube/2)
print,''
print,'Constructing a k-omega diagram of the input datacube..........'
print,''
; MAKE THE k-omega DIAGRAM USING THE PROVEN METHOD OF ROB RUTTEN
kopower = walsa_plotkopower_funct(datacube, sp_out, arcsecpx, cadence, apod=0.1,  kmax=1., fmax=1.)
   ; X SIZE STUFF
   xsize_kopower  = N_ELEMENTS(kopower[*,0])
   dxsize_kopower = spatial_Nyquist / FLOAT(xsize_kopower-1.)
   kopower_xscale = (FINDGEN(xsize_kopower)*dxsize_kopower) ; IN arcsec^-1
      ; Y SIZE STUFF
      ysize_kopower  = N_ELEMENTS(kopower[0,*])
      dysize_kopower = temporal_Nyquist / FLOAT(ysize_kopower-1.)
      kopower_yscale = (FINDGEN(ysize_kopower)*dysize_kopower) * 1000. ; IN mHz
Gaussian_kernel = GAUSSIAN_FUNCTION([0.65,0.65], WIDTH=3, MAXIMUM=1, /double)
Gaussian_kernel_norm = TOTAL(Gaussian_kernel,/nan)
kopower_plot = kopower
kopower_plot[*,1:*] = CONVOL(kopower[*,1:*],  Gaussian_kernel, Gaussian_kernel_norm, /edge_truncate)

; normalise to frequency resolution (in mHz)
freq = kopower_yscale[1:*]
if freq[0] eq 0 then freq0 = freq[1] else freq0 = freq[0]
kopower_plot = kopower_plot/freq0

if mode eq 0 then kopower_plot = ALOG10(kopower_plot)
if mode eq 2 then kopower_plot = SQRT(kopower_plot)

LOADCT, 0, /silent
!p.background = 255.
!p.color = 0.
x1 = 0.12 
x2 = 0.86 
y1 = 0.10
y2 = 0.80

!p.background = 255.
!p.color = 0.

; WHEN PLOTTING WE NEED TO IGNORE THE ZERO'TH ELEMENT (I.E., THE MEAN f=0) SINCE THIS WILL MESS UP THE LOG PLOT!
komegamap = (kopower_plot)[1:*,1:*]>MIN((kopower_plot)[1:*,1:*],/nan)<MAX((kopower_plot)[1:*,1:*],/nan)

IF silent EQ 0 THEN BEGIN

    if n_elements(komega) eq 0 then komega = 0 else komega = 1 
    if n_elements(clt) eq 0 then clt = 13 else clt=clt 
    ctload, clt, /silent 
    if n_elements(koclt) ne 0 then walsa_powercolor, koclt

    !p.background = 255.
    !p.color = 0.

    positioncb=[x1,y2+0.11,x2,y2+0.13]

    IF EPS eq 1 THEN BEGIN
        walsa_eps, size=[20,22]
        !p.font=0
        device,set_font='Times-Roman'
        !p.charsize=1.3
        !x.thick=4.
        !y.thick=4.
        !x.ticklen=-0.025
        !y.ticklen=-0.025
        positioncb=[x1,y2+0.12,x2,y2+0.14]
    ENDIF ELSE BEGIN
        IF (xscreensize ge 1000) AND (yscreensize ge 1000) THEN BEGIN 
            WINDOW, 0, xsize=1000, ysize=1000, title='QUEEFF: k-omega diagram'
            !p.charsize=1.7
            !p.charthick=1
            !x.thick=2
            !y.thick=2
            !x.ticklen=-0.025
            !y.ticklen=-0.025
        ENDIF
        IF (xscreensize lt 1000) OR  (yscreensize lt 1000) THEN BEGIN 
            WINDOW, 0, xsize=FIX(smallest_screensize*0.9), ysize=FIX(smallest_screensize*0.9), title='QUEEFF: k-omega diagram'
            !p.charsize=1
            !p.charthick=1
            !x.thick=2
            !y.thick=2
            !x.ticklen=-0.025
            !y.ticklen=-0.025       
        ENDIF
    ENDELSE

    walsa_pg_plotimage_komega, komegamap, kopower_xscale[1:*], kopower_yscale[1:*], noy2=noy2, nox2=nox2, smooth=smooth, $
        xtitle='Wavenumber (arcsec!U-1!N)', ytitle='Frequency (mHz)', xst=8, yst=8, xlog=xlog, ylog=ylog, position=[x1, y1, x2, y2], $
        xrange=xrange, yrange=yrange, threemin=threemin, fivemin=fivemin, eps=eps

    tickmarknames = STRARR(4)
    tickmarknames[0] = STRING(MIN(kopower_plot[1:*,1:*],/nan), FORMAT='(F5.1)')
    tickmarknames[1] = STRING(((MAX(kopower_plot[1:*,1:*],/nan)-MIN(kopower_plot[1:*,1:*],/nan)) * 0.33) + MIN(kopower_plot[1:*,1:*],/nan), FORMAT='(F5.1)')
    tickmarknames[2] = STRING(((MAX(kopower_plot[1:*,1:*],/nan)-MIN(kopower_plot[1:*,1:*],/nan)) * 0.67) + MIN(kopower_plot[1:*,1:*],/nan), FORMAT='(F4.1)')
    tickmarknames[3] = STRING(MAX(kopower_plot[1:*,1:*],/nan), FORMAT='(F4.1)')

    cgcolorbar, bottom=0, ncolors=255, divisions=3, minrange=MIN(kopower_plot[1:*,1:*],/nan), maxrange=MAX(kopower_plot[1:*,1:*],/nan), $
        position=positioncb, /top, ticknames=tickmarknames, title='Log!d10!n(Oscillation Power)', yticklen=0.00001

ENDIF

IF EPS eq 1 THEN walsa_endeps, filename=epsfilename, /noboundingbox

power = komegamap
wavenumber = kopower_xscale[1:*]
frequencies = kopower_yscale[1:*]

print, ' '
if filtering then print, ' ..... start filtering (in k-ω space)' else return
print, ' '

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; STEPS USED TO MAKE SURE THE FREQUENCIES ARE CHOSEN
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; NEED f1 AND k1
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, ''
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, 'Please click on the LOWEST frequency/wavenumber value you wish to preserve.....'
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN CURSOR, k1, f1, /data
WAIT, 1.0
; NEED f2 AND k2
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, ''
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, 'Please click on the HIGHEST frequency/wavenumber value you wish to preserve.....'
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN CURSOR, k2, f2, /data
WAIT, 1.0
; NEED ONLY f1 (spatial filtering ON)
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymin = 10^MIN(!y.crange,/nan)
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymax = 10^MAX(!y.crange,/nan)
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [k1, k1], [kopower_plot_ymin, kopower_plot_ymax], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [k2, k2], [kopower_plot_ymin, kopower_plot_ymax], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, ''
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, 'Please click on the LOWEST frequency value you wish to preserve inside the dotted lines.....'
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN CURSOR, nonsense, f1, /data
WAIT, 1.0
; NEED ONLY f2 (spatial filtering ON)
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymin = 10^MIN(!y.crange,/nan)
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymax = 10^MAX(!y.crange,/nan)
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [k1, k1], [kopower_plot_ymin, kopower_plot_ymax], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [k2, k2], [kopower_plot_ymin, kopower_plot_ymax], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, ''
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, 'Please click on the HIGHEST frequency value you wish to preserve inside the dotted lines.....'
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN CURSOR, nonsense, f2, /data
WAIT, 1.0
; NEED ONLY f1 (spatial filtering OFF)
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymin = 10^MIN(!y.crange,/nan)
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymax = 10^MAX(!y.crange,/nan)
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmin = 10^MIN(!x.crange,/nan)
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmax = 10^MAX(!x.crange,/nan)
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmin, kopower_plot_xmin], [kopower_plot_ymin, kopower_plot_ymax], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmax, kopower_plot_xmax], [kopower_plot_ymin, kopower_plot_ymax], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, ''
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, 'Please click on the LOWEST frequency value you wish to preserve inside the dotted lines.....'
IF NOT KEYWORD_SET(f1) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN CURSOR, nonsense, f1, /data
WAIT, 1.0
; NEED ONLY f2 (spatial filtering OFF)
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymin = 10^MIN(!y.crange,/nan)
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymax = 10^MAX(!y.crange,/nan)
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmin = 10^MIN(!x.crange,/nan)
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmax = 10^MAX(!x.crange,/nan)
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmin, kopower_plot_xmin], [kopower_plot_ymin, kopower_plot_ymax], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmax, kopower_plot_xmax], [kopower_plot_ymin, kopower_plot_ymax], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, ''
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, 'Please click on the HIGHEST frequency value you wish to preserve inside the dotted lines.....'
IF NOT KEYWORD_SET(f2) AND KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN CURSOR, nonsense, f2, /data
WAIT, 1.0
; NEED ONLY k1 (temporal filtering ON)
IF KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmin = 10^MIN(!x.crange,/nan)
IF KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmax = 10^MAX(!x.crange,/nan)
IF KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmin, kopower_plot_xmax], [f1, f1], line=1, thick=3, color=255
IF KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmin, kopower_plot_xmax], [f2, f2], line=1, thick=3, color=255
IF KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, ''
IF KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, 'Please click on the LOWEST wavenumber value you wish to preserve inside the dotted lines.....'
IF KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN CURSOR, k1, nonsense, /data
WAIT, 1.0
; NEED ONLY k2 (temporal filtering ON)
IF KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmin = 10^MIN(!x.crange,/nan)
IF KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmax = 10^MAX(!x.crange,/nan)
IF KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmin, kopower_plot_xmax], [f1, f1], line=1, thick=3, color=255
IF KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmin, kopower_plot_xmax], [f2, f2], line=1, thick=3, color=255
IF KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, ''
IF KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN print, 'Please click on the HIGHEST wavenumber value you wish to preserve inside the dotted lines.....'
IF KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND NOT KEYWORD_SET(no_temporal_filt) THEN CURSOR, k2, nonsense, /data
WAIT, 1.0
; NEED ONLY k1 (temporal filtering ON)
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymin = 10^MIN(!y.crange,/nan)
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymax = 10^MAX(!y.crange,/nan)
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmin = 10^MIN(!x.crange,/nan)
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmax = 10^MAX(!x.crange,/nan)
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmin, kopower_plot_xmax], [kopower_plot_ymin, kopower_plot_ymin], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmin, kopower_plot_xmax], [kopower_plot_ymax, kopower_plot_ymax], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN print, ''
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN print, 'Please click on the LOWEST wavenumber value you wish to preserve inside the dotted lines.....'
IF NOT KEYWORD_SET(f1) AND NOT KEYWORD_SET(k1) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN CURSOR, k1, nonsense, /data
WAIT, 1.0
; NEED ONLY k2 (temporal filtering ON)
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmin = 10^MIN(!x.crange,/nan)
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN kopower_plot_xmax = 10^MAX(!x.crange,/nan)
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymin = 10^MIN(!y.crange,/nan)
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN kopower_plot_ymax = 10^MAX(!y.crange,/nan)
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmin, kopower_plot_xmax], [kopower_plot_ymin, kopower_plot_ymin], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN PLOTS, [kopower_plot_xmin, kopower_plot_xmax], [kopower_plot_ymax, kopower_plot_ymax], line=1, thick=3, color=255
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN print, ''
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN print, 'Please click on the HIGHEST wavenumber value you wish to preserve inside the dotted lines.....'
IF NOT KEYWORD_SET(f2) AND NOT KEYWORD_SET(k2) AND NOT KEYWORD_SET(no_spatial_filt) AND KEYWORD_SET(no_temporal_filt) THEN CURSOR, k2, nonsense, /data
WAIT, 1.0
IF KEYWORD_SET(no_spatial_filt) THEN k1 = kopower_xscale[1]
IF KEYWORD_SET(no_spatial_filt) THEN k2 = MAX(kopower_xscale,/nan)
IF KEYWORD_SET(no_temporal_filt) THEN f1 = kopower_yscale[1]
IF KEYWORD_SET(no_temporal_filt) THEN f2 = MAX(kopower_yscale,/nan)
IF (k1 le 0.0) THEN k1 = kopower_xscale[1]
IF (k2 gt MAX(kopower_xscale,/nan)) THEN k2 = MAX(kopower_xscale,/nan)
IF (f1 le 0.0) THEN f1 = kopower_yscale[1]
IF (f2 gt MAX(kopower_yscale,/nan)) THEN f2 = MAX(kopower_yscale,/nan)
IF NOT KEYWORD_SET(no_spatial_filt) THEN BEGIN
    PLOTS, [k1, k2], [f1, f1], line=2, thick=2, color=255
    PLOTS, [k1, k2], [f2, f2], line=2, thick=2, color=255
    PLOTS, [k1, k1], [f1, f2], line=2, thick=2, color=255
    PLOTS, [k2, k2], [f1, f2], line=2, thick=2, color=255
ENDIF
IF KEYWORD_SET(no_spatial_filt) THEN BEGIN
    k1 = kopower_xscale[1]
    k2 = MAX(kopower_xscale,/nan)
    PLOTS, [k1, k2], [f1, f1], line=2, thick=2, color=255
    PLOTS, [k1, k2], [f2, f2], line=2, thick=2, color=255
    PLOTS, [k1, k1], [f1, f2], line=2, thick=2, color=255
    PLOTS, [k2, k2], [f1, f2], line=2, thick=2, color=255
ENDIF
print, ''
print, 'The preserved wavenumbers are ['+strtrim(k1, 2)+', '+strtrim(k2, 2)+'] arcsec^-1'
print, 'The preserved spatial sizes are ['+strtrim((2.*!pi)/k2, 2)+', '+strtrim((2.*!pi)/k1,2)+'] arcsec'
print,''
print, 'The preserved frequencies are ['+strtrim(f1, 2)+', '+strtrim(f2, 2)+'] mHz'
print, 'The preserved periods are ['+strtrim(FIX(1./(f2/1000.)), 2)+', '+strtrim(FIX(1./(f1/1000.)),2)+'] seconds'

pwavenumber = [k1,k2]
pspatialsize = [(2.*!pi)/k2,(2.*!pi)/k1]
pfrequency = [f1,f2]
pperiod = [FIX(1./(f2/1000.)),FIX(1./(f1/1000.))]

print,''
print,'Making a 3D Fourier transform of the input datacube..........'
threedft = FFT(datacube, -1, /double, /center)

; CALCULATE THE FREQUENCY AXES FOR THE 3D FFT
temp_x = FINDGEN((xsize_cube - 1)/2) + 1
is_N_even = (xsize_cube MOD 2) EQ 0
IF (is_N_even) THEN $
    spatial_frequencies_orig = ([0.0, temp_x, xsize_cube/2, -xsize_cube/2 + temp_x]/(xsize_cube*arcsecpx)) * (2.*!pi) $
    ELSE $
    spatial_frequencies_orig = ([0.0, temp_x, -(xsize_cube/2 + 1) + temp_x]/(xsize_cube*arcsecpx)) * (2.*!pi)
temp_x = FINDGEN((zsize_cube - 1)/2) + 1
is_N_even = (zsize_cube MOD 2) EQ 0
IF (is_N_even) THEN $
    temporal_frequencies_orig = [0.0, temp_x, zsize_cube/2, -zsize_cube/2 + temp_x]/(zsize_cube*cadence) $
    ELSE $
    temporal_frequencies_orig = [0.0, temp_x, -(zsize_cube/2 + 1) + temp_x]/(zsize_cube*cadence)

; NOW COMPENSATE THESE FREQUENCY AXES DUE TO THE FACT THE /center KEYWORD IS USED FOR THE FFT TRANSFORM
spatial_positive_frequencies = N_ELEMENTS(WHERE(spatial_frequencies_orig ge 0.))
IF N_ELEMENTS(spatial_frequencies_orig) MOD 2 EQ 0 THEN spatial_frequencies = SHIFT(spatial_frequencies_orig, (spatial_positive_frequencies-2))
IF N_ELEMENTS(spatial_frequencies_orig) MOD 2 NE 0 THEN spatial_frequencies = SHIFT(spatial_frequencies_orig, (spatial_positive_frequencies-1))
temporal_positive_frequencies = N_ELEMENTS(WHERE(temporal_frequencies_orig ge 0.))
IF N_ELEMENTS(temporal_frequencies_orig) MOD 2 EQ 0 THEN temporal_frequencies = SHIFT(temporal_frequencies_orig, (temporal_positive_frequencies-2))
IF N_ELEMENTS(temporal_frequencies_orig) MOD 2 NE 0 THEN temporal_frequencies = SHIFT(temporal_frequencies_orig, (temporal_positive_frequencies-1))

; ALSO NEED TO ENSURE THE threedft ALIGNS WITH THE NEW FREQUENCY AXES DESCRIBED ABOVE
IF N_ELEMENTS(temporal_frequencies_orig) MOD 2 EQ 0 THEN BEGIN
    FOR x = 0, (xsize_cube - 1) DO BEGIN
        FOR y = 0, (ysize_cube - 1) DO threedft[x, y, *] = SHIFT(REFORM(threedft[x,y,*]), -1)
    ENDFOR
ENDIF
IF N_ELEMENTS(spatial_frequencies_orig) MOD 2 EQ 0 THEN BEGIN
    FOR z = 0, (zsize_cube - 1) DO threedft[*, *, z] = SHIFT(REFORM(threedft[*,*,z]), [-1, -1])
ENDIF

; CONVERT FREQUENCIES AND WAVENUMBERS OF INTEREST INTO (FFT) DATACUBE PIXELS
pixel_k1_positive = walsa_closest(k1, spatial_frequencies_orig)
pixel_k2_positive = walsa_closest(k2, spatial_frequencies_orig)
pixel_f1_positive = walsa_closest(f1/1000., temporal_frequencies)
pixel_f2_positive = walsa_closest(f2/1000., temporal_frequencies)
pixel_f1_negative = walsa_closest(-f1/1000., temporal_frequencies)
pixel_f2_negative = walsa_closest(-f2/1000., temporal_frequencies)


torus_depth  = FIX((pixel_k2_positive[0] - pixel_k1_positive[0])/2.) * 2.
torus_center = FIX(((pixel_k2_positive[0] - pixel_k1_positive[0])/2.) + pixel_k1_positive[0])
IF KEYWORD_SET(spatial_torus) AND NOT KEYWORD_SET(no_spatial_filt) THEN BEGIN
    ; CREATE A FILTER RING PRESERVING EQUAL WAVENUMBERS FOR BOTH kx AND ky
    ; DO THIS AS A TORUS TO PRESERVE AN INTEGRATED GAUSSIAN SHAPE ACROSS THE WIDTH OF THE ANNULUS, THEN INTEGRATE ALONG 'z'
    spatial_torus = FLTARR(xsize_cube, ysize_cube, torus_depth)
    FOR i = 0, (FIX(torus_depth/2.)) DO BEGIN
        spatial_ring = (walsa_radial_distances([1,xsize_cube,ysize_cube]) LE (torus_center-i)) - $
                       (walsa_radial_distances([1,xsize_cube,ysize_cube]) LE (torus_center+i+1))
        spatial_ring[WHERE(spatial_ring gt 0.)] = 1.
        spatial_ring[WHERE(spatial_ring ne 1.)] = 0.
        spatial_torus[*,*,i] = spatial_ring
        spatial_torus[*,*,torus_depth-i-1] = spatial_ring
    ENDFOR
    ; INTEGRATE THROUGH THE TORUS TO FIND THE SPATIAL FILTER
    spatial_ring_filter      = TOTAL(spatial_torus,3,/nan) / FLOAT(torus_depth)
    spatial_ring_filter      = spatial_ring_filter / MAX(spatial_ring_filter,/nan) ; TO ENSURE THE PEAKS ARE AT 1.0
    ;IF N_ELEMENTS(spatial_frequencies_orig) MOD 2 EQ 0 THEN spatial_ring_filter = SHIFT(spatial_ring_filter, [-1,-1])
ENDIF

IF NOT KEYWORD_SET(spatial_torus) AND NOT KEYWORD_SET(no_spatial_filt) THEN BEGIN 
    spatial_ring_filter      = (walsa_radial_distances([1,xsize_cube,ysize_cube]) LE (torus_center-(FIX(torus_depth/2.)))) - $
                               (walsa_radial_distances([1,xsize_cube,ysize_cube]) LE (torus_center+(FIX(torus_depth/2.))+1))
    spatial_ring_filter      = spatial_ring_filter / MAX(spatial_ring_filter,/nan) ; TO ENSURE THE PEAKS ARE AT 1.0
    spatial_ring_filter[WHERE(spatial_ring_filter NE 1.)] = 0.
    ;IF N_ELEMENTS(spatial_frequencies_orig) MOD 2 EQ 0 THEN spatial_ring_filter = SHIFT(spatial_ring_filter, [-1,-1])
ENDIF

IF KEYWORD_SET(no_spatial_filt) THEN BEGIN
    spatial_ring_filter      = FLTARR(xsize_cube, ysize_cube)
    spatial_ring_filter[*]   = 1.
ENDIF

IF NOT KEYWORD_SET(no_temporal_filt) AND KEYWORD_SET(temporal_torus) THEN BEGIN
    ; CREATE A GAUSSIAN TEMPORAL FILTER TO PREVENT ALIASING
    temporal_filter    = FLTARR(zsize_cube)
    temporal_filter[*] = 0.
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 25) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(3, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 25) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 30) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(4, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 30) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 40) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(5, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 40) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 45) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(6, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 45) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 50) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(7, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 50) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 55) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(8, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 55) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 60) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(9, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 60) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 65) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(10, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 65) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 70) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(11, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 70) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 80) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(12, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 80) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 90) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(13, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 90) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 100) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(14, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 100) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 110) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(15, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 110) AND ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) LT 130) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(16, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    IF ((pixel_f2_positive[0]-pixel_f1_positive[0]+1) GE 130) THEN temporal_Gaussian = GAUSSIAN_FUNCTION(17, WIDTH=(pixel_f2_positive[0]-pixel_f1_positive[0]+1), MAXIMUM=1, /double)
    temporal_filter[pixel_f1_positive(0):pixel_f2_positive(0)] = temporal_Gaussian
    temporal_filter[pixel_f2_negative(0):pixel_f1_negative(0)] = temporal_Gaussian
    temporal_filter = temporal_filter / MAX(temporal_filter,/nan) ; TO ENSURE THE PEAKS ARE AT 1.0
ENDIF

IF NOT KEYWORD_SET(no_temporal_filt) AND NOT KEYWORD_SET(temporal_torus) THEN BEGIN
    temporal_filter    = FLTARR(zsize_cube)
    temporal_filter[*] = 0.
    temporal_filter[pixel_f1_positive(0):pixel_f2_positive(0)] = 1.0
    temporal_filter[pixel_f2_negative(0):pixel_f1_negative(0)] = 1.0
ENDIF


IF KEYWORD_SET(no_temporal_filt) THEN BEGIN
    temporal_filter    = FLTARR(zsize_cube)
    temporal_filter[*] = 1.
ENDIF


; MAKE SOME FIGURES FOR PLOTTING - MAKES THINGS AESTHETICALLY PLEASING!
torus_map = MAKE_MAP(spatial_ring_filter, dx=spatial_frequencies[1]-spatial_frequencies[0], dy=spatial_frequencies[1]-spatial_frequencies[0], xc=0, yc=0, time='', units='arcsecs')
spatial_fft = TOTAL(threedft, 3,/nan)
spatial_fft_map = MAKE_MAP(ALOG10(spatial_fft), dx=spatial_frequencies[1]-spatial_frequencies[0], dy=spatial_frequencies[1]-spatial_frequencies[0], xc=0, yc=0, time='', units='arcsecs')
spatial_fft_filtered = spatial_fft * spatial_ring_filter
spatial_fft_filtered_map = MAKE_MAP(ALOG10(spatial_fft_filtered>1e-15), dx=spatial_frequencies[1]-spatial_frequencies[0], dy=spatial_frequencies[1]-spatial_frequencies[0], xc=0, yc=0, time='', units='arcsecs')
temporal_fft = TOTAL(TOTAL(threedft, 2,/nan), 1)
IF (xscreensize ge 1000) AND (yscreensize ge 1000) THEN WINDOW, 1, xsize=1500, ysize=1000, title='QUEEFF: FFT filter specs'
IF (xscreensize lt 1000) OR  (yscreensize lt 1000) THEN WINDOW, 1, xsize=smallest_screensize, ysize=FIX(smallest_screensize*0.8), title='QUEEFF: FFT filter specs'
x1 = 0.07
x2 = 0.33
x3 = 0.40
x4 = 0.66
x5 = 0.72
x6 = 0.98
y1 = 0.07
y2 = 0.47
y3 = 0.56
y4 = 0.96
LOADCT, 5, /silent
IF (xscreensize ge 1000) AND (yscreensize ge 1000) THEN BEGIN 
    plot_map, spatial_fft_map, charsize=2, xticklen=-.025, yticklen=-.025, xtitle='Wavenumber (k!Dx!N ; arcsec!U-1!N)', ytitle='Wavenumber (k!Dy!N ; arcsec!U-1!N)', title='Spatial FFT', dmin=MIN(spatial_fft_map.data,/nan)+1., dmax=MAX(spatial_fft_map.data,/nan)-1., position=[x1, y3, x2, y4]
    PLOTS, [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], [0, 0], line=2, thick=2, color=0
    PLOTS, [0, 0], [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], line=2, thick=2, color=0
    LOADCT,0,/silent
    plot_map, torus_map, charsize=2, xticklen=-.025, yticklen=-.025, xtitle='Wavenumber (k!Dx!N ; arcsec!U-1!N)', ytitle='', title='Spatial FFT filter', dmin=0, dmax=1, position=[x3, y3, x4, y4], /noerase
    PLOTS, [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], [0, 0], line=2, thick=2, color=255
    PLOTS, [0, 0], [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], line=2, thick=2, color=255
    LOADCT,5,/silent
    plot_map, spatial_fft_filtered_map, charsize=2, xticklen=-.025, yticklen=-.025, xtitle='Wavenumber (k!Dx!N ; arcsec!U-1!N)', ytitle='', title='Filtered spatial FFT', dmin=MIN(spatial_fft_map.data,/nan)+1., dmax=MAX(spatial_fft_map.data,/nan)-1., position=[x5, y3, x6, y4], /noerase
    PLOTS, [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], [0, 0], line=2, thick=2, color=255
    PLOTS, [0, 0], [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], line=2, thick=2, color=255
    PLOT, temporal_frequencies*1000., ABS(temporal_fft), /ylog, xst=1, xticklen=-.026, yticklen=-.011, charsize=2, xtitle='Frequency (mHz)', ytitle='Power (arb. units)', position=[x1+0.05, y1, x6, y2], /noerase
    temporal_fft_plot_ymin = 10^MIN(!y.crange,/nan)
    temporal_fft_plot_ymax = 10^MAX(!y.crange,/nan)
    PLOTS, [0, 0], [temporal_fft_plot_ymin, temporal_fft_plot_ymax], line=2, thick=2, color=0
    LOADCT,39,/silent
    OPLOT, temporal_frequencies*1000., (temporal_filter)>temporal_fft_plot_ymin, line=2, color=55, thick=2
    OPLOT, temporal_frequencies*1000., (ABS(temporal_fft * temporal_filter))>temporal_fft_plot_ymin, line=0, color=254, thick=2
    LOADCT,5,/silent
    WAIT, 0.5
ENDIF
IF (xscreensize lt 1000) OR  (yscreensize lt 1000) THEN BEGIN 
    plot_map, spatial_fft_map, charsize=1, xticklen=-.025, yticklen=-.025, xtitle='Wavenumber (k!Dx!N ; arcsec!U-1!N)', ytitle='Wavenumber (k!Dy!N ; arcsec!U-1!N)', title='Spatial FFT', dmin=MIN(spatial_fft_map.data,/nan)+1., dmax=MAX(spatial_fft_map.data,/nan)-1., position=[x1, y3, x2, y4]
    PLOTS, [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], [0, 0], line=2, thick=2, color=0
    PLOTS, [0, 0], [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], line=2, thick=2, color=0
    LOADCT,0,/silent
    plot_map, torus_map, charsize=1, xticklen=-.025, yticklen=-.025, xtitle='Wavenumber (k!Dx!N ; arcsec!U-1!N)', ytitle='', title='Spatial FFT filter', dmin=0, dmax=1, position=[x3, y3, x4, y4], /noerase
    PLOTS, [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], [0, 0], line=2, thick=2, color=255
    PLOTS, [0, 0], [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], line=2, thick=2, color=255
    LOADCT,5,/silent
    plot_map, spatial_fft_filtered_map, charsize=1, xticklen=-.025, yticklen=-.025, xtitle='Wavenumber (k!Dx!N ; arcsec!U-1!N)', ytitle='', title='Filtered spatial FFT', dmin=MIN(spatial_fft_map.data,/nan)+1., dmax=MAX(spatial_fft_map.data,/nan)-1., position=[x5, y3, x6, y4], /noerase
    PLOTS, [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], [0, 0], line=2, thick=2, color=255
    PLOTS, [0, 0], [MIN(spatial_frequencies,/nan), MAX(spatial_frequencies,/nan)], line=2, thick=2, color=255
    PLOT, temporal_frequencies*1000., ABS(temporal_fft), /ylog, xst=1, charsize=1, xticklen=-.026, yticklen=-.011, xtitle='Frequency (mHz)', ytitle='Power (arb. units)', position=[x1+0.05, y1, x6, y2], /noerase
    temporal_fft_plot_ymin = 10^MIN(!y.crange,/nan)
    temporal_fft_plot_ymax = 10^MAX(!y.crange,/nan)
    PLOTS, [0, 0], [temporal_fft_plot_ymin, temporal_fft_plot_ymax], line=2, thick=2, color=0
    LOADCT,39,/silent
    OPLOT, temporal_frequencies*1000., (temporal_filter)>temporal_fft_plot_ymin, line=2, color=55, thick=2
    OPLOT, temporal_frequencies*1000., (ABS(temporal_fft * temporal_filter))>temporal_fft_plot_ymin, line=0, color=254, thick=2
    LOADCT,5,/silent
    WAIT, 0.5
ENDIF

; APPLY THE GAUSSIAN FILTERS TO THE DATA TO PREVENT ALIASING
FOR i = 0, (zsize_cube-1) DO threedft[*,*,i] = REFORM(threedft[*,*,i]) * spatial_ring_filter
FOR x = 0, (xsize_cube-1) DO BEGIN 
    FOR y = 0, (ysize_cube-1) DO BEGIN 
        threedft[x,y,*] = REFORM(threedft[x,y,*]) * temporal_filter
    ENDFOR
ENDFOR

; ALSO NEED TO ENSURE THE threedft ALIGNS WITH THE OLD FREQUENCY AXES USED BY THE /center CALL
IF N_ELEMENTS(temporal_frequencies_orig) MOD 2 EQ 0 THEN BEGIN
    FOR x = 0, (xsize_cube - 1) DO BEGIN
        FOR y = 0, (ysize_cube - 1) DO threedft[x, y, *] = SHIFT(REFORM(threedft[x,y,*]), 1)
    ENDFOR
ENDIF
IF N_ELEMENTS(spatial_frequencies_orig) MOD 2 EQ 0 THEN BEGIN
    FOR z = 0, (zsize_cube - 1) DO threedft[*, *, z] = SHIFT(REFORM(threedft[*,*,z]), [1, 1])
ENDIF

new_cube = REAL_PART(FFT(threedft, 1, /double, /center))

LOADCT,0, /silent

filtered_cube = new_cube

PRINT
if (mode eq 0) then print,' mode = 0: log(power)'
if (mode eq 1) then print,' mode = 1: linear power' 
if (mode eq 2) then print,' mode = 2: sqrt(power)'

!P.Multi = 0
Cleanplot, /Silent

print, ''
print, 'COMPLETED!'
print,''

END

This code uses the following routine (originanly from Rob Rutten) to compute the k-ω power.

WaLSA_plotkopower_funct.pro

; -----------------------------------------------------------------------------------------------------
; WaLSAtools: Wave analysis tools
; Copyright (C) 2025 WaLSA Team - Shahin Jafarzadeh et al.
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
; 
; http://www.apache.org/licenses/LICENSE-2.0
; 
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
; 
; Note: If you use WaLSAtools for research, please consider citing:
; Jafarzadeh, S., Jess, D. B., Stangalini, M. et al. 2025, Nature Reviews Methods Primers, in press.
; -----------------------------------------------------------------------------------------------------


; NAME: plotkopower --> walsa_plotkopower_funct
;
; PURPOSE: plot "k-omega" diagram = Fourier power at temporal frequency
;   f against horizontal spatial wavenumber k_h
;
; CALLING SEQUENCE:
;   plotkopower, cube, arcsecpx, cadence, plotfilename, $
;     apod=apod, $
;     kmax=kmax, fmax=fmax, minpower=minpower, maxpower=maxpower, $
;     contours=contours, lamb=lamb, fundamental=fundamental

;
; INPUTS:
;   cube: (x,y,t) data cube, any type
;   arcsecpx = angular image scale in arcsec/px
;   cadence: [regular] sampling cadence in sec
;   plotfilename: postscript output file name
;     optional keywords: 
;   apod: fractional extent of apodization edges; default 0.1
;   kmax: maximum k_h axis as fraction of Nyquist value, default 0.2
;   fmax: maximum f axis as fraction of Nyquist value, default 0.5
;   minpower: minimum of power plot range, default maxpower-5
;   maxpower: maximum of power plot range, default alog10(max(power))
;   contours: set true to plot contours
;   lamb: value > 0 overplots Lamb line omeha = c_s kn at c_s = lamb km/s 
;   fundamental: set true to overplot fundamental mode omega=sqrt(g k_h)
;
; MODIFICATION HISTORY: 
;  Sep 2010: Rob Rutten (RR) assembly of Alfred de Wijn's routines.
;  Okt 2010: RR optional overplots Lamb line and fundamental mode
;  SJ modified: based on latest version of this code (Mar 2 2012): Mar 2011: RR -1 > +1 line 162 from Alfred de Wijn
;-

;-----------------------------------------------------------------------
function avgstd, array, stdev=stdev
    ; get array average + standard deviation
    ;RR Aug 23 2010 found in Mabula Haverkamp's IDL, later also AdW
    ;RR not the same as avg.pro in ssw
    ;RR not the same as avg.pro in Pit Suetterlin DOT software
    ;RR so renamed to avgstd
  avrg = total(array,/nan)/n_elements(array)
  stdev = sqrt(total((array-avrg)^2,/nan)/n_elements(array))
  return, avrg
end

;----------------------------------------------------------------------------
function linear, x, p            ;RR used in temporal detrending 
    ymod = p[0] + x * p[1]
    return, ymod
end

;----------------------------------------------------------------------------
function gradient, x, y, p       ;RR used in spatail detrending
    zmod = p[0] + x * p[1] + y * p[2]
    return, zmod
end

;----------------------------------------------------------------------------
function apod3dcube,cube,apod
    ; apodizes cube in all three coordinates, with detrending 

  ; get cube dimensions
  sizecube=size(cube)
  nx=sizecube[1]
  ny=sizecube[2]
  nt=sizecube[3]
  apocube = fltarr(nx,ny,nt)

  ; define temporal apodization
  apodt = fltarr(nt)+1
  if (apod ne 0) then begin
    apodrimt = nt*apod
    apodt[0] = (sin(!pi/2.*findgen(apodrimt)/apodrimt))^2
    apodt = apodt*shift(rotate(apodt,2),1)   ;RR had ik nooit verzonnen
  endif 

  ; temporal detrending, not per px, only mean-image trend 
  ttrend = fltarr(nt)
  tf = findgen(nt) + 1
  for it=0, nt-1 do begin
    img = cube[*,*,it]
    ttrend[it] = avgstd(img)
  endfor
  fitp = mpfitfun('linear', tf, ttrend, fltarr(nt)+1, [1000.,0.],/quiet)
  fit = fitp[0] + tf * fitp[1]

  ; temporal apodization per (x,y) column
  ;RR do not reinsert trend to keep [0,0] Fourier pixel from dominating
  for it=0, nt-1 do begin
    img = cube[*,*,it]
    apocube[*,*,it] = (img-fit[it])*apodt[it]  ;RR  + ttrend[it]
  endfor

  ; define spatial apodization
  apodx = fltarr(nx)+1
  apody = fltarr(ny)+1
  if (apod ne 0) then begin
    apodrimx=apod*nx
    apodrimy=apod*ny
    apodx[0] = (sin(!pi/2.*findgen(apodrimx)/apodrimx))^2
    apody[0] = (sin(!pi/2.*findgen(apodrimy)/apodrimy))^2
    apodx = apodx*shift(rotate(apodx,2),1)
    apody = apody*shift(rotate(apody,2),1)
    apodxy = apodx # apody
 endif

 if (apod eq 0) then begin
    apodxy = apodx # apody
 endif

  ; spatial gradient removal + apodizing per image
  xf = fltarr(nx,ny)+1.
  yf = xf
  for it=0, nt-1 do begin
    img = apocube[*,*,it]
    avg = avgstd(img)
    ;RR mpfit2dfun = ssw/gen/idl/fitting/mpfit/mpfit2dfun.pro
    fitp = mpfit2dfun('gradient',xf,yf,img,fltarr(nx,ny)+1,[1000.,0.,0.],/quiet)
    fit = fitp[0]+xf*fitp[1]+yf*fitp[2]
    apocube[*,*,it] = (img-fit)*apodxy + avg
  endfor

  ; done
  return,apocube
end

;---------------------------------------------------------------------------
function ko_dist, sx, sy, double=double
    ; set up Pythagoras distance array from origin 
    ;RR from Alfred de Wijn email Aug 30 2010 
  dx = rebin(dindgen(sx/2+1)/(sx/2),sx/2+1,sy/2+1)
  dy = rotate(rebin(dindgen(sy/2+1)/(sy/2),sy/2+1,sx/2+1),1)
  dxy = sqrt(dx^2+dy^2)*(min([sx,sy],/nan)/2+1.)
  afstanden = dblarr(sx,sy)
  afstanden[0,0] = dxy
  ; get other quadrants
  afstanden[sx/2,0] = rotate(dxy[1:*,*],5)       ;RR 5 = 90 deg
  afstanden[0,sy/2] = rotate(dxy[*,1:*],7)       ;RR 7 = 270 deg
  afstanden[sx/2,sy/2] = rotate(dxy[1:*,1:*],2)  ;RR 2 = 180 deg
  if not keyword_set(double) then afstanden = fix(round(afstanden))
  return, afstanden
end

;---------------------------------------------------------------------------
function averpower,cube
  ; compute 2D (k_h,f) power array by circular averaging over k_x, k_y

  ; get cube dimensions
  sizecube=size(cube)
  nx=sizecube[1]
  ny=sizecube[2]
  nt=sizecube[3]

  ; forward fft and throw away half of it
  ; perform fft in time direction first
  fftcube = (fft(fft((fft(cube,-1, dimension=3))[*,*,0:nt/2],-1,dimension=1),dimension=2))

  ; set up distances 
  fftfmt = size(fftcube)
  afstanden = ko_dist(nx,ny)   ;RR integer-rounded Pythagoras array
  ; maxdist = min([nx,ny])/2-1   ;RR largest quarter circle
  maxdist = min([nx,ny],/nan)/2+1   ;RR largest quarter circle --> from RR on Mar 2011!

  ; get average power over all k_h distances, building power(k_h,f)
  avpow = fltarr(maxdist+1,nt/2+1)
  for i=0, maxdist do begin
    waar = where(afstanden eq i)
    for j=0, nt/2 do begin
      w1 = (fftcube[*,*,j])[waar]
      avpow[i,j] = total(w1*conj(w1),/nan)/n_elements(waar)
    endfor
;;  writeu,-1,string(format='(%"\rcomputing avpow... ",i6,"/",i6)',i,maxdist)
  endfor

  ; done
  return, avpow
end

;---------------------------------------------------------------------------
pro koplotpow,avpow,arcsecpx,cadence,kmax,fmax,minpower,maxpower
    ; plotting program
  sizepow = size(avpow)

  ; select extent = fractions of Nyquist values
  plotrange = [fix(sizepow[1]*kmax),fix(sizepow[2]*fmax)]
  plotpow = avpow[0:plotrange[0]-1,0:plotrange[1]-1]

  ;RR 5x5 resizing, I guess for better tick positioning and contours
  xas = 2.*!pi/(arcsecpx*2)*findgen(plotrange[0])/(sizepow[1]-1)
  rexas = 2.*!pi/(arcsecpx*2)*findgen(plotrange[0]*5)/(sizepow[1]*5-1)
  yas = 1./(cadence*2)*findgen(plotrange[1])/(sizepow[2]-1)*1e3
  reyas = 1./(cadence*2)*findgen(plotrange[1]*5)/(sizepow[2]*5-1)*1e3
  plotpowrebin = convol(rebin(plotpow,plotrange[0]*5,plotrange[1]*5),fltarr(6,6)+1/(6.*6.),/edge_truncate)
;  plotpowrebin=plotpow
  xrange = [min(xas,/nan),max(xas,/nan)]
  yrange = [min(yas,/nan),max(yas,/nan)]


tvframe,alog10(plotpowrebin) > minpower < maxpower,/ba,/sa,/as,xrange=xrange,yrange=yrange,xtitle='k_h [arcsec!U-1!N]',ytitle='f [mHz]'



  ; plot wavelength axis along top
  if (xrange[1] lt 10) then begin ;RR I wonder how to automate this  
    wavtickspos = [10, 5, 3, 2, 1.5, 1]
    wavticksn = ['10','5','3','2','1.5','1']
  endif else if (xrange[1] lt 20) then begin
    wavtickspos = [10, 5, 3, 2, 1.5, 1, 0.5]
    wavticksn = ['10','5','3','2','1.5','1','0.5']
  endif else if (xrange[1] lt 50) then begin
    wavtickspos = [5.0, 2.0, 1.0, 0.5, 0.2]
    wavticksn = ['5','2','1','0.5','0.2']
  endif else begin
    wavtickspos = [5.0, 2.0, 1.0, 0.5, 0.2, 0.1, 0.05]
    wavticksn = ['5','2','1','0.5','0.2','0.1','0.05']
  endelse
  wavticks=n_elements(wavtickspos)-1
  wavticksv = 2.*!pi/wavtickspos   ;RR wavelength from circle frequency

  axis, /xaxis, xticks=wavticks, xtickv=wavticksv, xtickname=wavticksn,ticklen=-0.015/0.53, xminor=1, xtitle='wavelength [arcsec]'

  ; plot period axis along righthand side 
  if (yrange[1] lt 10) then begin ;RR I wonder how to automate this 
    pertickspos = [20, 10, 5, 3, 2, 1]
    perticksn = ['20','10','5','3','2','1']
  endif else if (yrange[1] lt 20) then begin
    pertickspos = [10, 5, 3, 2, 1.5, 1, 0.5]
    perticksn = ['10', '5','3','2','1.5','1','0.5']
  endif else if (yrange[1] lt 50) then begin
    pertickspos = [10, 5, 2, 1, 0.5, 0.2, 0.1]
    perticksn = ['10','5','2','1','0.5','0.2','0.1']
  endif else if (yrange[1] lt 100) then begin
    pertickspos = [2, 1, 0.5, 0.2, 0.1]
    perticksn = ['2','1','0.5','0.2','0.1']
  endif else begin
    pertickspos = [0.5, 0.2, 0.1, 0.05, 0.02]
    perticksn = ['0.5','0.2','0.1','0.05','0.02']
  endelse
  perticks=n_elements(pertickspos)-1
  perticksv = 1e3/pertickspos/60.  ;RR period, from mHz to min
  axis, /yaxis, yticks=perticks, ytickv=perticksv, ytickname=perticksn,ticklen=-0.015/0.7, yminor=1, ytitle='period [min]'



end

; --------------------------- main part ------------------------------

FUNCTION walsa_plotkopower_funct,cube,sp_out,arcsecpx,cadence,apod=apod,kmax=kmax,fmax=fmax,minpower=minpower,maxpower=maxpower
    ; wrapper calling the above subroutines

if (n_elements(apod) ne 0) then apod=apod else apod=0.1
if (n_elements(kmax) ne 0) then kmax=kmax else kmax=1.0
if (n_elements(fmax) ne 0) then fmax=fmax else fmax=1.0

if (kmax gt 1) then kmax=1
if (fmax gt 1) then fmax=1



; apodize the cube
apocube=apod3dcube(cube,apod)

; compute radially-averaged power
avpow=averpower(apocube)

sp_out=avpow


; set min, max cutoffs
maxavpow=alog10(max(avpow,/nan))
minavpow=alog10(min(avpow,/nan))
print, ' max log(power) = ', maxavpow, ' min log(power) = ',minavpow
if (n_elements(maxpower) ne 0) then maxpower=maxpower else maxpower=maxavpow
if (n_elements(minpower) ne 0) then minpower=minpower else minpower=maxpower-5

;print,kmax,fmax

; plot the diagram
;koplotpow,avpow,arcsecpx,cadence,kmax,fmax,minpower,maxpower

; done

RETURN,sp_out


end

B-ω Analysis

WaLSA_bomega

This routine computes and plots B-ω diagram, based on the approach introduced in this scientific article.

WaLSA_bomega.pro

; -----------------------------------------------------------------------------------------------------
; WaLSAtools: Wave analysis tools
; Copyright (C) 2025 WaLSA Team - Shahin Jafarzadeh et al.
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
; 
; http://www.apache.org/licenses/LICENSE-2.0
; 
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
; 
; Note: If you use WaLSAtools for research, please consider citing:
; Jafarzadeh, S., Jess, D. B., Stangalini, M. et al. 2025, Nature Reviews Methods Primers, in press.
; -----------------------------------------------------------------------------------------------------


;+
; NAME: WaLSA_Bomega
;       part of -- WaLSAtools --
;
; PURPOSE:   
;   Compute (and plot) B-ω diagram: average FFT power spectra 
;            for various magnetic-field bins within the FoV
;            Note: (1) The entire cube is detrended (linearly) 
;            and apodized (in both spatial and temporal domains)
;            prior to the FFT analysis
;            (2) Plotted B-ω is smoothed for better visibility
;
; CALLING SEQUENCE:
;   walsa_bomega, datacube, Bmap, time, power=power, frequencies=frequencies, Barray=Barray
;
; + INPUTS:
;   datacube:   datacube in the form of [x, y, t]
;   bmap:       magnetic-field map (in G), in the form of [x, y]
;               -- should spatially be the same size as the datacube
;   time:       time of the observations in sec
;
; + OPTIONAL KEYWORDS:
;   binsize:         magnetic-field bin size (default: 50 G)
;   silent:          if set, the B-ω map is not displayed.
;   clt:             color table number (idl ctload)
;   koclt:           custom color tables for k-ω diagram (currently available: 1 and 2)
;   threemin:        if set, a horizontal line marks the three-minute periodicity
;   fivemin:         if set, a horizontal line marks the five-minute periodicity
;   xlog:            if set, x-axis (magnetic field) is plotted in logarithmic scale
;   ylog:            if set, y-axis (frequency) is plotted in logarithmic scale
;   xrange:          x-axis (wavenumber) range
;   yrange:          y-axis (frequency) range
;   noy2:            if set, 2nd y-axis (period, in sec) is not plotted
;                    (p = 1000/frequency)
;   smooth:          if set, power is smoothed
;   normalizedbins   if set, power at each bin is normalized to its maximum value
;                    (this facilitates visibility of relatively small power)
;   xtickinterval    x-asis (i.e., magnetic fields) tick intervals in G (default: 400 G)   
;   epsfilename:     if provided (as a string), an eps file of the k-ω diagram is made
;   mode:            outputted power mode: 0 = log(power) (default), 1 = linear power, 2 = sqrt(power) = amplitude
; ---- detrending, and apodization parameters----
;   apod:           extent of apodization edges (of a Tukey window); default 0.1
;   nodetrendapod:  if set, neither detrending nor apodization is performed!
;   pxdetrend:      subtract linear trend with time per pixel. options: 1=simple, 2=advanced; default: 2
;   polyfit:        the degree of polynomial fit to the data to detrend it.
;                   if set, instead of linear fit this polynomial fit is performed.
;   meantemporal:   if set, only a very simple temporal detrending is performed by subtracting the mean signal from the signal.
;                   i.e., the fitting procedure (linear or higher polynomial degrees) is omitted.
;   meandetrend:    if set, subtract linear trend with time for the image means (i.e., spatial detrending)
;   recon:          optional keyword that will Fourier reconstruct the input timeseries.
;                   note: this does not preserve the amplitudes and is only useful when attempting 
;                   to examine frequencies that are far away from the 'untrustworthy' low frequencies.
;
; + OUTPUTS:
;   power:          B-ω map, a stack of average power spectra (in magnetic-field bins)
;                   along y axis -> The x and y axes are B (in G) and 
;                   frequency (in mHz); in dn^2/mhz, i.e., normalized to frequency resolution (see mode for the scale)
;   frequencies:    frequencies of the power spectra 
;                   (i.e., values of the y axis of the B-ω map)
;   barray:         magnetic-field values of the middle of the bins
;                   (i.e., values of the x axis of the B-ω map)
;
;
; + CREDITS:
;   Author: Shahin Jafarzadeh, March 2021. 
;   Note: if YOU USE THIS CODE, then PLEASE CITE THE ORIGINAL PUBLICATION WHERE IT WAS INTRODUCED:
;          Stangalini et al. 2021, A&A, in press (https://ui.adsabs.harvard.edu/abs/2021arXiv210311639S)
;-

pro walsa_bomega, datacube, Bmap, cadence=cadence, time=time, binsize=binsize, power=power, frequencies=frequencies, Barray=Barray, $
                  silent=silent, clt=clt, koclt=koclt, threemin=threemin, fivemin=fivemin, xlog=xlog, ylog=ylog, $ ; plotting keywords
                  xrange=xrange, yrange=yrange, epsfilename=epsfilename, noy2=noy2, smooth=smooth, normalizedbins=normalizedbins, $
                  xtickinterval=xtickinterval, mode=mode

if n_elements(cadence) eq 0 then cadence=walsa_mode(walsa_diff(time))

nx = N_ELEMENTS(datacube[*,0,0])
ny = N_ELEMENTS(datacube[0,*,0])
nt = N_ELEMENTS(datacube[0,0,*])

temporal_Nyquist = 1. / (cadence * 2.)

print,' '
print,'The input datacube is of size: ['+strtrim(nx,2)+', '+strtrim(ny,2)+', '+strtrim(nt,2)+']'
print,' '
print,'Temporally, the important values are:'
print,'    2-element duration (Nyquist period) = '+strtrim((cadence * 2.),2)+' seconds'
print,'    Time series duration = '+strtrim(cadence*nt,2)+' seconds'
print,'    Nyquist frequency = '+strtrim(temporal_Nyquist*1000.,2)+' mHz'
print, ' '

nxb = N_ELEMENTS(Bmap[*,0,0])
nyb = N_ELEMENTS(Bmap[0,*,0])

dimensions = GET_SCREEN_SIZE(RESOLUTION=resolution)
xscreensize = dimensions[0]
yscreensize = dimensions[1]
IF (xscreensize le yscreensize) THEN smallest_screensize = xscreensize
IF (yscreensize le xscreensize) THEN smallest_screensize = yscreensize

if nx gt nxb OR ny gt nyb then begin
    print, ' [!] The datacube and B-map must have the same [x,y] size.'
    print, ' '
    stop
endif

if n_elements(binsize) eq 0 then binsize = 50. ; in G
if n_elements(silent) eq 0 then silent = 0
if n_elements(noy2) eq 0 then noy2 = 0
if n_elements(normalizedbins) eq 0 then normalizedbins = 0 else normalizedbins = 1 
if n_elements(epsfilename) eq 0 then eps = 0 else eps = 1 
if n_elements(xtickinterval) eq 0 then xtickinterval = 400 ; in G
if n_elements(mode) eq 0 then mode=0
if n_elements(xlog) eq 0 then xlog = 0
if n_elements(ylog) eq 0 then ylog = 0
if n_elements(nodetrendapod) eq 0 then nodetrendapod = 0 else nodetrendapod = 1

Bmap = ABS(Bmap)

brange = minmax(Bmap, /nan)
nbin = floor((brange[1]-brange[0])/binsize)

; detrend and apodize the cube
if nodetrendapod eq 0 then begin
    print, ' '
    print, ' -- Detrend and apodize the cube .....'
    datacube=walsa_detrend_apod(datacube,apod,meandetrend,pxdetrend,polyfit=polyfit,meantemporal=meantemporal,recon=recon,cadence=cadence) 
endif

frequencies = 1./(cadence*2)*findgen(nt/2+1)/(nt/2)
nff=n_elements(frequencies)
frequencies = frequencies[1:nff-1]
frequencies = frequencies*1000. ; in mHz

numf = n_elements(frequencies)

Barray = fltarr(nbin)
bopower = fltarr(nbin,numf)
fac = 0

PRINT
for i=0L, nbin-1 do begin
    ii = where(Bmap le brange[1]-(fac*float(binsize)) AND Bmap gt brange[1]-((fac+1)*float(binsize)),num)

    Barray[i] = (brange[1]-(fac*float(binsize)))-(float(binsize)/2.)

    if num gt 0 then begin
        coords = array_indices(Bmap,ii)
        xx = reform(coords[0,*])  &  yy = reform(coords[1,*])
        nxy = n_elements(xx)

        poweravg = fltarr(numf)
        for ixy=0L, nxy-1 do begin
            poweri = (2.*(ABS((fft(reform(datacube[xx[ixy],yy[ixy],*]),-1))[0:nt/2.])^2))/frequencies[0] ; in DN^2/mHz
            poweravg = poweravg + poweri[1:nff-1]
        endfor
        poweravg = poweravg/float(nxy)

        if normalizedbins eq 1 then poweravg = 100.*poweravg/max(poweravg)

        bopower[i,*] = poweravg
    endif

    print,string(13b)+'.... % finished (through bins, from larger B): ',float(i)*100./(nbin-1),format='(a,f4.0,$)'

    if brange[1]-((fac+1)*float(binsize)) le 0 then break
    fac = fac+1
endfor
PRINT

bopower = reverse(bopower,1)
Barray = reverse(Barray)
nb = n_elements(Barray)

Barray[0]=brange[0]
Barray[nb-1]=brange[1]

ppp = bopower

Gaussian_kernel = GAUSSIAN_FUNCTION([0.65,0.65], WIDTH=3, MAXIMUM=1, /double)
Gaussian_kernel_norm = TOTAL(Gaussian_kernel,/nan)
bopower = CONVOL(bopower,  Gaussian_kernel, Gaussian_kernel_norm, /edge_truncate)

if mode eq 0 then bopower = ALOG10(bopower)
if mode eq 2 then bopower = SQRT(bopower)

nlevels = 256
step = (Max(bopower) - Min(bopower)) / nlevels
vals = Indgen(nlevels) * step + Min(bopower)

bopower = (bopower)[0:*,0:*]>MIN((bopower)[0:*,0:*],/nan)<MAX((bopower)[0:*,0:*],/nan)

if silent eq 0 then begin

    LOADCT, 0, /silent
    !p.background = 255.
    !p.color = 0.
    x1 = 0.12 
    x2 = 0.86 
    y1 = 0.10
    y2 = 0.80

    if n_elements(clt) eq 0 then clt = 13 else clt=clt 
    ctload, clt, /silent 
    if n_elements(koclt) ne 0 then walsa_powercolor, koclt

    !p.background = 255.
    !p.color = 0.

    positioncb=[x1,y2+0.05,x2,y2+0.07]

    if EPS eq 1 then begin
        walsa_eps, size=[20,22]
        !p.font=0
        device,set_font='Times-Roman'
        !p.charsize=1.3
        !x.thick=4.
        !y.thick=4.
        !x.ticklen=-0.025
        !y.ticklen=-0.025
    endif else begin
        if (xscreensize ge 1000) AND (yscreensize ge 1000) then begin 
            WINdoW, 0, xsize=1000, ysize=1000, title='B-omega diagram'
            !p.charsize=1.7
            !p.charthick=1
            !x.thick=2
            !y.thick=2
            !x.ticklen=-0.025
            !y.ticklen=-0.025
        endif
        if (xscreensize lt 1000) OR  (yscreensize lt 1000) then begin 
            WINdoW, 0, xsize=FIX(smallest_screensize*0.9), ysize=FIX(smallest_screensize*0.9), title='B-omega diagram'
            !p.charsize=1
            !p.charthick=1
            !x.thick=2
            !y.thick=2
            !x.ticklen=-0.025
            !y.ticklen=-0.025       
        endif
    endelse

    walsa_pg_plotimage_komega, bopower, Barray, frequencies, $
        ytitle='Frequency (mHz)', xtitle='B (G)', xst=1, yst=8, xlog=xlog, ylog=ylog, xrange=xrange, yrange=yrange, eps=eps, $
        position=[x1, y1, x2, y2], noy2=noy2, nox2=1, smooth=smooth, threemin=threemin, fivemin=fivemin, xtickinter=xtickinterval

    tickmarknames = STRARR(4)
    tickmarknames[0] = STRING(MIN(bopower[1:*,1:*],/nan), FORMAT='(F5.1)')
    tickmarknames[1] = STRING(((MAX(bopower[1:*,1:*],/nan)-MIN(bopower[1:*,1:*],/nan)) * 0.33) + MIN(bopower[1:*,1:*],/nan), FORMAT='(F5.1)')
    tickmarknames[2] = STRING(((MAX(bopower[1:*,1:*],/nan)-MIN(bopower[1:*,1:*],/nan)) * 0.67) + MIN(bopower[1:*,1:*],/nan), FORMAT='(F4.1)')
    tickmarknames[3] = STRING(MAX(bopower[1:*,1:*],/nan), FORMAT='(F4.1)')

    if normalizedbins eq 1 then cbtitle='Log!d10!n(Normalized Oscillation Power)' else cbtitle='Log!d10!n(Oscillation Power)'

    cgcolorbar, bottom=0, ncolors=255,  minrange=MIN(bopower[1:*,1:*],/nan), maxrange=MAX(bopower[1:*,1:*],/nan), /top, $
        ticknames=tickmarknames, yticklen=0.00001, position=positioncb, title=cbtitle

    if EPS eq 1 then walsa_endeps, filename=epsfilename, /noboundingbox
endif

power = bopower

PRINT
if (mode eq 0) then print,' mode = 0: log(power)'
if (mode eq 1) then print,' mode = 1: linear power' 
if (mode eq 2) then print,' mode = 2: sqrt(power)'

!P.Multi = 0
Cleanplot, /Silent

PRINT
PRINT, 'COMPLETED!'
PRINT

end

Detrending and Apodisation

WaLSA_detrend_apod

All signals are detrended (linearly, or using higher-order polynomial fits if desired) and apodised (using a Tukey window, i.e., tapered cosine) prior to all spectral analyses (unless otherwise it is omitted). Here is the code used for detrending and apodising the signals. The spatial apodisation for the k-ω diagram, are performed inside the WaLSA_plotkopower_funct.pro.

WaLSA_detrend_apod.pro

; -----------------------------------------------------------------------------------------------------
; WaLSAtools: Wave analysis tools
; Copyright (C) 2025 WaLSA Team - Shahin Jafarzadeh et al.
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
; 
; http://www.apache.org/licenses/LICENSE-2.0
; 
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
; 
; Note: If you use WaLSAtools for research, please consider citing:
; Jafarzadeh, S., Jess, D. B., Stangalini, M. et al. 2025, Nature Reviews Methods Primers, in press.
; -----------------------------------------------------------------------------------------------------


;+
; NAME: WaLSA_detrend_apod
;       part of -- WaLSAtools --
;
; PURPOSE:
;   Detrend and apodise time series to account for their nonlinear and nonstationary nature
;   (apodisation is also known as 'windowing' in some other communities)
;
; DESCRIPTION:
;   Each 1D time series is detrended by subtracting the signal from its fitted linear or higher 
;   polynomial degrees function (or only from the mean signal), and finally apodised by the 
;   Tukey (tapered cosine) window to bring the edges to zero. If set, the edge effects are 
;   more conservatively examined by means of a wavelet-based approach (walsa_wave_recon.pro)
;
; CALLING SEQUENCE:
;   corrected_cube = walsa_detrend_apod(cube,apod,meandetrend,pxdetrend)
;
; + INPUTS:
;   data:           1D time series, or (x,y,t) datacube, any type
;                   (an ordered sequence of data points, typically some variable quantity 
;                   measured at successive times.)
;   apod:           extent of apodization edges of the Tukey (tapered cosine) window
;                   default: 0.1
;   pxdetrend:      subtract linear trend with time per pixel. options: 1=simple, 2=advanced
;                   default: 2
;
; + OPTIONAL KEYWORDS:
;   polyfit:            the degree of polynomial fit to the data to detrend it.
;                       if set, instead of linear fit this polynomial fit is performed.
;   meantemporal:       if set, only a very simple temporal detrending is performed by subtracting 
;                       the mean signal from the signal.
;                       i.e., the fitting procedure (linear or higher polynomial degrees) is omitted.
;   meandetrend:        if set, subtract linear trend with time for the image means 
;                       (i.e., spatial detrending)
;   recon:              optional keyword that will Fourier reconstruct the input timeseries.
;                       note: this does not preserve the amplitudes and is only useful when attempting 
;                       to examine frequencies that are far away from the 'untrustworthy' low frequencies.
;   cadence:            cadence of the observations. it is required when recon is set.
;   resample_original   if recon is set, then this keyword allow setting a range (i.e., min_resample and max_resample)
;                       to which the unpreserved amplitudes are approximated.
;   min_resample        minimum value for resample_original. Default: min of each 1D array (time series) in data.
;   max_resample        maximum value for resample_original. Default: max of each 1D array (time series) in data.
; 
; + OUTPUTS:
;   corrected_cube:     The detrended and apodised cube
;
; MODIFICATION HISTORY
;
;  2010 plotpowermap: Rob Rutten, assembly of Alfred de Wijn's routines 
;                     (https://webspace.science.uu.nl/~rutte101/rridl/cubelib/plotpowermap.pro)
;  2018-2021 modified/extended by Shahin Jafarzadeh & David B. Jess
;-

FUNCTION linear, x, p   
  ; used by mpfitfun
  ymod=p[0] + x * p[1] 
  return, ymod
END

FUNCTION walsa_detrend_apod,cube,apod,meandetrend,pxdetrend,polyfit=polyfit,meantemporal=meantemporal,$
                            recon=recon,cadence=cadence,resample_original=resample_original,min_resample=min_resample,$
                            max_resample=max_resample,silent=silent, dj=dj, lo_cutoff=lo_cutoff, hi_cutoff=hi_cutoff, upper=upper

  if (n_elements(apod) ne 0) then apod=apod else apod=0.1
  if (n_elements(polyfit) eq 0) then apolyfit=0 else apolyfit=1
  if not keyword_set(meandetrend) then meandetrend=0
  if not keyword_set(silent) then silent=0
  if (n_elements(pxdetrend) ne 0) then pxdetrend=pxdetrend else pxdetrend=2

  if silent eq 0 then begin
      print, ' '
      print, ' -- Detrend and apodize the cube .....'
      print, ' '
  endif

  sizecube = size(cube)
  if sizecube[0] ne 3 then begin
      if sizecube[0] eq 1 then begin
          blablacube = fltarr(1,1,sizecube[1])
          blablacube[0,0,*] = cube
          cube = blablacube
      endif else begin
          print, ' '
          print, ' [!] The datacube must have either 1 or 3 dimension(s).'
          print, ' '
          stop
      endelse
  endif

  sizecube=size(cube)
  nx=sizecube[1]
  ny=sizecube[2]
  nt=sizecube[3]
  apocube=cube
  tf=findgen(nt) + 1
  col=fltarr(nt)
  apodt = fltarr(nt)+1
  if (apod ne 0) then begin ; Tukey window
      apodrim = apod*nt
      apodt[0] = (sin(!pi/2.*findgen(apodrim)/apodrim))^2
      apodt = apodt*shift(rotate(apodt,2),1)
  endif 

  ; meandetrend: get spatially-averaged trend 
  fit=0  
  if (meandetrend) then begin
      avgseq=fltarr(nt)
      for it=0,nt-1 do avgseq[it]=total(cube[*,*,it])
      avgseq=avgseq/(double(nx)*double(ny)) 
      meanfitp = mpfitfun('linear',tf,avgseq,fltarr(nt)+1,[1000.,0.],/quiet)
      meanfit=meanfitp[0]+tf*meanfitp[1]-total(avgseq)/double(nt)
  endif 

  ; apodize per [x,y] temporal column = time sequence per pixel
  for ix=long(0),long(nx)-1 do begin  
      for iy=long(0),long(ny)-1 do begin
          col=cube[ix,iy,*]
          IF KEYWORD_SET(recon) THEN col = walsa_wave_recon(reform(col), cadence, dj=dj, lo_cutoff=lo_cutoff, hi_cutoff=hi_cutoff, upper=upper)
          meancol=walsa_avgstd(col)
          if (meandetrend) then col=col-meanfit
          if n_elements(meantemporal) eq 0 then begin 
              if apolyfit eq 0 then begin
                  if (pxdetrend eq 1) then begin
                      pxfitp=poly_fit(tf,col,1) 
                      col=col-pxfitp[0]-tf*pxfitp[1]+meancol  
                  endif
                  if (pxdetrend eq 2) then begin
                      pxfitp = mpfitfun('linear',tf,col,fltarr(nt)+1,[meancol, 0.],/quiet)
                      col=col-pxfitp[0]-tf*pxfitp[1]+meancol  
                  endif
              endif else begin
                  lc_fit = GOODPOLY(FINDGEN(nt), col, polyfit, 3, lc_yfit, lc_newx, lc_newy)
                  col=col-lc_yfit
              endelse
          endif
          ocol=(col-meancol)*apodt+meancol
          if not KEYWORD_SET(min_resample) then min_resample = min(cube[ix,iy,*])
          if not KEYWORD_SET(max_resample) then max_resample = max(cube[ix,iy,*])
          IF KEYWORD_SET(recon) THEN if KEYWORD_SET(resample_original) then ocol = scale_vector(ocol,min_resample,max_resample)
          apocube[ix,iy,*] = ocol
      endfor
      if silent eq 0 then if long(nx) gt 1 then if (pxdetrend ne 0) then $ 
         writeu,-1,string(format='(%"\r == detrend next row... ",i5,"/",i5)',ix+1,nx)
  endfor
  if silent eq 0 then if long(nx) gt 1 then if (pxdetrend ne 0) then PRINT

  return, reform(apocube)
END

Wavelet Analysis

WaLSA_wavelet

A modified/extended variant of wavelet.pro (of Torrence & Compo) to compute wavelet power spectrum and its related parameters.

WaLSA_wavelet.pro

; -----------------------------------------------------------------------------------------------------
; WaLSAtools: Wave analysis tools
; Copyright (C) 2025 WaLSA Team - Shahin Jafarzadeh et al.
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
; 
; http://www.apache.org/licenses/LICENSE-2.0
; 
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
; 
; Note: If you use WaLSAtools for research, please consider citing:
; Jafarzadeh, S., Jess, D. B., Stangalini, M. et al. 2025, Nature Reviews Methods Primers, in press.
; -----------------------------------------------------------------------------------------------------


;+
; NAME: WaLSA_wavelet
;       part of -- WaLSAtools --
;       modification of WAVELET.pro
;
; PURPOSE:   
;   Compute the WAVELET transform of a 1D time series.
;           Based on the original code from Torrenc & Compo
;           (see Torrence, C. and G. P. Compo, 1998: A Practical Guide to
;           Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.)
;
; CALLING SEQUENCE:
;   wave = WaLSA_wavelet(Y,DT)
;
; INPUTS:
;
;    Y = the time series of length N.
;
;    DT = amount of time between each Y value, i.e. the sampling time.
;
; OUTPUTS:
;
;    WAVE is the WAVELET transform of Y. This is a complex array
;    of dimensions (N,J+1). FLOAT(WAVE) gives the WAVELET amplitude,
;    ATAN(IMAGINARY(WAVE),FLOAT(WAVE)) gives the WAVELET phase.
;    The WAVELET power spectrum is ABS(WAVE)^2.
;
;
; OPTIONAL KEYWORD INPUTS:
;
;    S0 = the smallest scale of the wavelet.  Default is 2*DT.
;
;    DJ = the spacing between discrete scales. Default is 0.125.
;         A smaller # will give better scale resolution, but be slower to plot.
;
;    J = the # of scales minus one. Scales range from S0 up to S0*2^(J*DJ),
;        to give a total of (J+1) scales. Default is J = (LOG2(N DT/S0))/DJ.
;
;    MOTHER = A string giving the mother wavelet to use.
;            Currently, 'Morlet','Paul','DOG' (derivative of Gaussian)
;            are available. Default is 'Morlet'.
;
;    PARAM = optional mother wavelet parameter.
;            For 'Morlet' this is k0 (wavenumber), default is 6.
;            For 'Paul' this is m (order), default is 4.
;            For 'DOG' this is m (m-th derivative), default is 2.
;
;    PAD = if set, then pad the time series with enough zeroes to get
;         N up to the next higher power of 2. This prevents wraparound
;         from the end of the time series to the beginning, and also
;         speeds up the FFT's used to do the wavelet transform.
;         This will not eliminate all edge effects (see COI below).
;
;    LAG1 = LAG 1 Autocorrelation, used for SIGNIF levels. Default is 0.0
;
;    SIGLVL = significance level to use. Default is 0.95
;
;    VERBOSE = if set, then print out info for each analyzed scale.
;
;    RECON = if set, then reconstruct the time series, and store in Y.
;            Note that this will destroy the original time series,
;            so be sure to input a dummy copy of Y.
;
;    FFT_THEOR = theoretical background spectrum as a function of
;                Fourier frequency. This will be smoothed by the
;                wavelet function and returned as a function of PERIOD.
;
;   colornoise = if set, noise background is based on Auchère et al. 2017, ApJ, 838, 166 / 2016, ApJ, 825, 110
;
; OPTIONAL KEYWORD OUTPUTS:
;
;    PERIOD = the vector of "Fourier" periods (in time units) that corresponds
;           to the SCALEs.
;
;    POWER = Wavelet power spectrum
;
;    SCALE = the vector of scale indices, given by S0*2^(j*DJ), j=0...J
;            where J+1 is the total # of scales.
;
;    COI = if specified, then return the Cone-of-Influence, which is a vector
;        of N points that contains the maximum period of useful information
;        at that particular time.
;        Periods greater than this are subject to edge effects.
;        This can be used to plot COI lines on a contour plot by doing:
;            IDL>  CONTOUR,wavelet,time,period
;            IDL>  PLOTS,time,coi,NOCLIP=0
;
;    YPAD = returns the padded time series that was actually used in the
;         wavelet transform.
;
;    DAUGHTER = if initially set to 1, then return the daughter wavelets.
;         This is a complex array of the same size as WAVELET. At each scale
;         the daughter wavelet is located in the center of the array.
;
;    SIGNIF = output significance levels as a function of PERIOD
;
;    FFT_THEOR = output theoretical background spectrum (smoothed by the
;                wavelet function), as a function of PERIOD.
;
;    Plot = if set, the wavelet power spectrum is plotted.
;           
;    colorct = the IDL color table number. Default: 20
;           
;    w = window number (for IDL). Default: 6
;           
; ---- detrending, and apodization parameters----
;   apod:           extent of apodization edges (of a Tukey window); default 0.1
;   nodetrendapod:  if set, neither detrending nor apodization is performed!
;   pxdetrend:      subtract linear trend with time per pixel. options: 1=simple, 2=advanced; default: 2
;   polyfit:        the degree of polynomial fit to the data to detrend it.
;                   if set, instead of linear fit this polynomial fit is performed.
;   meantemporal:   if set, only a very simple temporal detrending is performed by subtracting the mean signal from the signal.
;                   i.e., the fitting procedure (linear or higher polynomial degrees) is omitted.
;   meandetrend:    if set, subtract linear trend with time for the image means (i.e., spatial detrending)
;
; [ Defunct INPUTS:
; [   OCT = the # of octaves to analyze over.           ]
; [         Largest scale will be S0*2^OCT.             ]
; [         Default is (LOG2(N) - 1).                   ]
; [   VOICE = # of voices in each octave. Default is 8. ]
; [          Higher # gives better scale resolution,    ]
; [          but is slower to plot.                     ]
; ]
;
;;----------------------------------------------------------------------------
;
; EXAMPLE:
;
;    IDL> ntime = 256
;    IDL> y = RANDOMN(s,ntime)       ;*** create a random time series
;    IDL> dt = 0.25
;    IDL> time = FINDGEN(ntime)*dt   ;*** create the time index
;    IDL> 
;    IDL> wave = WaLSA_wavelet(y,dt,PERIOD=period,PAD=1,COI=coi,MOTHER='Morlet',/RECON,dj=0.025,scale=scale,SIGNIF=SIGNIF,SIGLVL=0.99,/apod,/plot)
;
;;----------------------------------------------------------------------------
; This routine is originally based on WAVELET.pro
; Copyright (C) 1995-2004, Christopher Torrence and Gilbert P. Compo
;
; This software may be used, copied, or redistributed as long as it is not
; sold and this copyright notice is reproduced on each copy made.
; This routine is provided as is without any express or implied warranties
; whatsoever.
;
; Notice: Please acknowledge the use of the above software in any publications:
;    ``Wavelet software was provided by C. Torrence and G. Compo,
;      and is available at URL: http://paos.colorado.edu/research/wavelets/''.
;
; Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
;            Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.
;
; Please send a copy of such publications to either C. Torrence or G. Compo:
;  Dr. Christopher Torrence               Dr. Gilbert P. Compo
;  Research Systems, Inc.                 Climate Diagnostics Center
;  4990 Pearl East Circle                 325 Broadway R/CDC1
;  Boulder, CO 80301, USA                 Boulder, CO 80305-3328, USA
;  E-mail: chris[AT]rsinc[DOT]com         E-mail: compo[AT]colorado[DOT]edu
;;----------------------------------------------------------------------------
; Modified/extended by Shahin Jafarzadeh 2016-2021
;-

FUNCTION morlet, $ ;*********************************************** MORLET
    k0,scale,k,period,coi,dofmin,Cdelta,psi0

    IF (k0 EQ -1) THEN k0 = 6d
    n = N_ELEMENTS(k)
    expnt = -(scale*k - k0)^2/2d*(k GT 0.)
    dt = 2*!PI/(n*k(1))
    norm = SQRT(2*!PI*scale/dt)*(!PI^(-0.25))   ; total energy=N   [Eqn(7)]
    morlet = norm*EXP(expnt > (-100d))
    morlet = morlet*(expnt GT -100)  ; avoid underflow errors
    morlet = morlet*(k GT 0.)  ; Heaviside step function (Morlet is complex)
    fourier_factor = (4*!PI)/(k0 + SQRT(2+k0^2)) ; Scale-->Fourier [Sec.3h]
    period = scale*fourier_factor
    coi = fourier_factor/SQRT(2)   ; Cone-of-influence [Sec.3g]
    dofmin = 2   ; Degrees of freedom with no smoothing
    Cdelta = -1
    IF (k0 EQ 6) THEN Cdelta = 0.776 ; reconstruction factor
    psi0 = !PI^(-0.25)
;   PRINT,scale,n,SQRT(TOTAL(ABS(morlet)^2,/DOUBLE))
    RETURN,morlet
END

FUNCTION paul, $ ;************************************************** PAUL
    m,scale,k,period,coi,dofmin,Cdelta,psi0

    IF (m EQ -1) THEN m = 4d
    n = N_ELEMENTS(k)
    expnt = -(scale*k)*(k GT 0.)
    dt = 2d*!PI/(n*k(1))
    norm = SQRT(2*!PI*scale/dt)*(2^m/SQRT(m*FACTORIAL(2*m-1)))
    paul = norm*((scale*k)^m)*EXP(expnt > (-100d))*(expnt GT -100)
    paul = paul*(k GT 0.)
    fourier_factor = 4*!PI/(2*m+1)
    period = scale*fourier_factor
    coi = fourier_factor*SQRT(2)
    dofmin = 2   ; Degrees of freedom with no smoothing
    Cdelta = -1
    IF (m EQ 4) THEN Cdelta = 1.132 ; reconstruction factor
    psi0 = 2.^m*FACTORIAL(m)/SQRT(!PI*FACTORIAL(2*m))
;   PRINT,scale,n,norm,SQRT(TOTAL(paul^2,/DOUBLE))*SQRT(n)
    RETURN,paul
END

FUNCTION dog, $ ;*************************************************** DOG
    m,scale,k,period,coi,dofmin,Cdelta,psi0

    IF (m EQ -1) THEN m = 2
    n = N_ELEMENTS(k)
    expnt = -(scale*k)^2/2d
    dt = 2d*!PI/(n*k(1))
    norm = SQRT(2*!PI*scale/dt)*SQRT(1d/GAMMA(m+0.5))
    I = DCOMPLEX(0,1)
    gauss = -norm*(I^m)*(scale*k)^m*EXP(expnt > (-100d))*(expnt GT -100)
    fourier_factor = 2*!PI*SQRT(2./(2*m+1))
    period = scale*fourier_factor
    coi = fourier_factor/SQRT(2)
    dofmin = 1   ; Degrees of freedom with no smoothing
    Cdelta = -1
    psi0 = -1
    IF (m EQ 2) THEN BEGIN
        Cdelta = 3.541 ; reconstruction factor
        psi0 = 0.867325
    ENDIF
    IF (m EQ 6) THEN BEGIN
        Cdelta = 1.966 ; reconstruction factor
        psi0 = 0.88406
    ENDIF
;   PRINT,scale,n,norm,SQRT(TOTAL(ABS(gauss)^2,/DOUBLE))*SQRT(n)
    RETURN,gauss
END

;****************************************************************** WAVELET
FUNCTION walsa_wavelet,y1,dt, $   ;*** required inputs
    S0=s0,DJ=dj,J=j, $   ;*** optional inputs
    PAD=pad,MOTHER=mother,PARAM=param, $
    VERBOSE=verbose,NO_WAVE=no_wave,RECON=recon, $
    LAG1=lag1,SIGLVL=siglvl,DOF=dof,GLOBAL=global, $   ;*** optional inputs
    SCALE=scale,PERIOD=period,YPAD=ypad, $  ;*** optional outputs
    DAUGHTER=daughter,COI=coi, removespace=removespace, koclt=koclt, $
    SIGNIF=signif,FFT_THEOR=fft_theor, $
    OCT=oct,VOICE=voice,log=log,silent=silent, normal=normal, $
    plot=plot,colorct=colorct,w=w, apod=apod, nodetrendapod=nodetrendapod, $
    pxdetrend=pxdetrend, meandetrend=meandetrend,power=power,$
    polyfit=polyfit,meantemporal=meantemporal,colornoise=colornoise,$
    clt=clt,epsfilename=epsfilename

    ON_ERROR,2
    r = CHECK_MATH(0,1)
    n = N_ELEMENTS(y1)
    n1 = n
    base2 = FIX(ALOG(n)/ALOG(2) + 0.4999)   ; power of 2 nearest to N

    ;....check keywords & optional inputs
    if n_elements(log) eq 0 THEN log = 0
    if n_elements(pad) eq 0 THEN pad = 1
    IF (N_ELEMENTS(s0) LT 1) THEN s0 = 2.0*dt
    IF (N_ELEMENTS(voice) EQ 1) THEN dj = 1./voice
    IF (N_ELEMENTS(dj) LT 1) THEN dj = 0.025
    if n_elements(colornoise) eq 0 then colornoise=0
    IF (N_ELEMENTS(oct) EQ 1) THEN J = FLOAT(oct)/dj
    IF (N_ELEMENTS(J) LT 1) THEN J=FIX((ALOG(FLOAT(n)*dt/s0)/ALOG(2))/dj)  ;[Eqn(10)]
    IF (N_ELEMENTS(mother) LT 1) THEN mother = 'MORLET'
    IF (N_ELEMENTS(param) LT 1) THEN param = -1
    IF (N_ELEMENTS(siglvl) LT 1) THEN siglvl = 0.95
    IF (N_ELEMENTS(lag1) LT 1) THEN lag1 = 0.0
    if n_elements(plot) eq 0 then plot = 0
    if n_elements(nodetrendapod) eq 0 then nodetrendapod=0
    lag1 = lag1(0)
    verbose = KEYWORD_SET(verbose)
    do_daughter = KEYWORD_SET(daughter)
    do_wave = NOT KEYWORD_SET(no_wave)
    recon = KEYWORD_SET(recon)

    if colornoise ne 0 then begin ; Auchère et al. 2017, ApJ, 838, 166 / 2016, ApJ, 825, 110
        nt = n_elements(y1)
        J = FIX(alog(nt/2.0)/alog(2)/dj) 
        s0 = 2*dt
    endif

    IF KEYWORD_SET(global) THEN MESSAGE, $
        'Please use WAVE_SIGNIF for global significance tests'

    ; detrend and apodize the cube
    if nodetrendapod eq 0 then $
        y1=walsa_detrend_apod(y1,apod,meandetrend,pxdetrend,polyfit=polyfit,meantemporal=meantemporal,cadence=dt,silent=silent)

;....construct time series to analyze, pad if necessary
    ypad = y1 - TOTAL(y1)/n    ; remove mean
    IF KEYWORD_SET(pad) THEN BEGIN   ; pad with extra zeroes, up to power of 2
        ypad = [ypad,FLTARR(2L^(base2 + 1) - n)]
        n = N_ELEMENTS(ypad)
    ENDIF

;....construct SCALE array & empty PERIOD & WAVE arrays
    na = J + 1                  ; # of scales
    scale = DINDGEN(na)*dj      ; array of j-values
    scale = 2d0^(scale)*s0      ; array of scales  2^j   [Eqn(9)]
    period = FLTARR(na,/NOZERO) ; empty period array (filled in below)
    wave = COMPLEXARR(n,na,/NOZERO)  ; empty wavelet array
    IF (do_daughter) THEN daughter = wave   ; empty daughter array

;....construct wavenumber array used in transform [Eqn(5)]
    k = (DINDGEN(n/2) + 1)*(2*!PI)/(DOUBLE(n)*dt)
    k = [0d,k,-REVERSE(k(0:(n-1)/2 - 1))]

;....compute FFT of the (padded) time series
    yfft = FFT(ypad,-1,/DOUBLE)  ; [Eqn(3)]

    IF (verbose) THEN BEGIN  ;verbose
        PRINT
        PRINT,mother
        PRINT,'#points=',n1,'   s0=',s0,'   dj=',dj,'   J=',FIX(J)
        IF (n1 NE n) THEN PRINT,'(padded with ',n-n1,' zeroes)'
        PRINT,['j','scale','period','variance','mathflag'], $
            FORMAT='(/,A3,3A11,A10)'
    ENDIF  ;verbose
    IF (N_ELEMENTS(fft_theor) EQ n) THEN fft_theor_k = fft_theor ELSE $
        fft_theor_k = (1-lag1^2)/(1-2*lag1*COS(k*dt)+lag1^2)  ; [Eqn(16)]
    fft_theor = FLTARR(na)

;....loop thru each SCALE
    FOR a1=0,na-1 DO BEGIN  ;scale
            psi_fft = CALL_FUNCTION(mother, param, scale(a1), k, period1, coi, dofmin, Cdelta, psi0)
        IF (do_wave) THEN $
            wave(*,a1) = FFT(yfft*psi_fft,1,/DOUBLE)  ;wavelet transform[Eqn(4)]
        period(a1) = period1   ; save period
        fft_theor(a1) = TOTAL((ABS(psi_fft)^2)*fft_theor_k)/n
        IF (do_daughter) THEN $
            daughter(*,a1) = FFT(psi_fft,1,/DOUBLE)   ; save daughter
        IF (verbose) THEN PRINT,a1,scale(a1),period(a1), $
                TOTAL(ABS(wave(*,a1))^2),CHECK_MATH(0), $
                FORMAT='(I3,3F11.3,I6)'
    ENDFOR  ;scale

    coi = coi*[FINDGEN((n1+1)/2),REVERSE(FINDGEN(n1/2))]*dt   ; COI [Sec.3g]

    IF (do_daughter) THEN $   ; shift so DAUGHTERs are in middle of array
        daughter = [daughter(n-n1/2:*,*),daughter(0:n1/2-1,*)]

;....significance levels [Sec.4]
    sdev = (MOMENT(y1))(1)
    fft_theor = sdev*fft_theor  ; include time-series variance
    dof = dofmin
    signif = fft_theor*CHISQR_CVF(1. - siglvl,dof)/dof   ; [Eqn(18)]

    IF (recon) THEN BEGIN  ; Reconstruction [Eqn(11)]
        IF (Cdelta EQ -1) THEN BEGIN
            y1 = -1
            MESSAGE,/INFO, $
                'Cdelta undefined, cannot reconstruct with this wavelet'
        ENDIF ELSE BEGIN
            y1=dj*SQRT(dt)/(Cdelta*psi0)*(FLOAT(wave) # (1./SQRT(scale)))
            y1 = y1[0:n1-1]
        ENDELSE
    ENDIF

    time = findgen(n1)*dt
    amplitude = wave(0:n1-1,*) ; get rid of padding before returning amplitudes
    power = ABS(amplitude)^2

    if plot ne 0 then begin
        powplot = power
        perplot = period
        timplot = time
        coiplot = coi
        sigplot = signif
        walsa_plot_wavelet_spectrum, powplot, perplot, timplot, coiplot, sigplot, clt=clt, w=w, log=log, removespace=removespace, $
                                     koclt=koclt, normal=normal, epsfilename=epsfilename
    endif

    RETURN,amplitude 

END

This code also uses the following routine to plot the wavelet power spectrum (along with confidence levels and cone-of-influence regions).

WaLSA_plot_wavelet_spectrum.pro

; -----------------------------------------------------------------------------------------------------
; WaLSAtools: Wave analysis tools
; Copyright (C) 2025 WaLSA Team - Shahin Jafarzadeh et al.
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
; 
; http://www.apache.org/licenses/LICENSE-2.0
; 
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
; 
; Note: If you use WaLSAtools for research, please consider citing:
; Jafarzadeh, S., Jess, D. B., Stangalini, M. et al. 2025, Nature Reviews Methods Primers, in press.
; -----------------------------------------------------------------------------------------------------


pro walsa_plot_wavelet_spectrum, power, period, time, coi, significance, clt=clt, w=w, log=log, removespace=removespace, $
                                 koclt=koclt, normal=normal, epsfilename=epsfilename, maxperiod=maxperiod, mHz=mHz

if n_elements(log) eq 0 then log=1
if n_elements(mHz) eq 0 then mHz=1
if n_elements(removespace) eq 0 then removespace=0
if n_elements(normal) eq 0 then normal=0
if n_elements(epsfilename) eq 0 then eps = 0 else eps = 1 

nt = n_elements(reform(time)) & np = n_elements(reform(period))
significance = REBIN(TRANSPOSE(significance),nt,np)

fundf = 1000./(time[nt-1]) ; fundamental frequency (frequency resolution) in mHz
if n_elements(maxperiod) eq 0 then maxp = 1000./fundf else maxp = maxperiod ; longest period to be plotted
if n_elements(maxperiod) eq 0 then if removespace ne 0 then maxp = max(coi) ; remove areas below the COI
iit = closest_index(maxp,period)
period = period[0:iit]
significance = reform(significance[*,0:iit])
power = reform(power[*,0:iit])

power = reverse(power,2)
significance = reverse(significance,2)

sigi = power/significance

if n_elements(w) eq 0 then w=6

dimensions = GET_SCREEN_SIZE(RESOLUTION=resolution)
xscreensize = dimensions[0]
yscreensize = dimensions[1]
IF (xscreensize le yscreensize) THEN smallest_screensize = xscreensize
IF (yscreensize le xscreensize) THEN smallest_screensize = yscreensize

if EPS eq 1 then begin
    walsa_eps, size=[18,13]
    !p.font=0
    device,set_font='Times-Roman'
    charsize = 1.3
    !x.thick=4.
    !y.thick=4.
    !x.ticklen=-0.033
    !y.ticklen=-0.021
    barthick = 550
    distbar = 550
    coithick = 3.
    arrowsize = 20.
    arrowthick = 3.5
    c_thick = 3.
    h_thick = 1.4
    ; arrowheadsize = 10.
endif else begin
    if (xscreensize ge 1000) AND (yscreensize ge 1000) then begin 
        window, w, xs=900, ys = 650, title=strtrim(w,2)+': Wavelet Power Spectrum'
        charsize=2.0
        !x.thick=2.  
        !y.thick=2.  
        !x.ticklen=-0.033
        !y.ticklen=-0.021
        ; !X.MINOR=6
        distbar = 30
        barthick = 30
        coithick = 2
        c_thick = 2.
        h_thick = 1.
    endif
    if (xscreensize lt 1000) OR  (yscreensize lt 1000) then begin 
        window, w, xs=FIX(smallest_screensize*0.9), ys=FIX(smallest_screensize*0.9), title=strtrim(w,2)+': Wavelet Power Spectrum'
        charsize=1.7
        !x.thick=2
        !y.thick=2
        !x.ticklen=-0.033
        !y.ticklen=-0.021
        distbar = 25
        barthick = 25
        coithick = 2
        c_thick = 2.
        h_thick = 1.
    endif
endelse

colset
device,decomposed=0

xtitle = 'Time (s)'
ztitle = 'Power!C'
if log ne 0 then ztitle = 'Log!d10!n(Power)!C'
if normal ne 0 then begin
    ztitle = 'Normalised Power!C'
    if log ne 0 then ztitle = 'Log!d10!n(Normalised Power)!C'
endif

ii = where(power lt 0., cii)
if cii gt 0 then power(ii) = 0.
if normal ne 0 then power = 100.*power/max(power)

xrg = [min(time),max(time)]
yrg = [max(period),min(period)]

if n_elements(clt) eq 0 then clt=20
loadct, clt
if n_elements(koclt) ne 0 then walsa_powercolor, koclt

if log ne 0 then power = alog10(power)

walsa_image_plot, power, xrange=xrg, yrange=yrg, $
          nobar=0, zrange=minmax(power,/nan), /ylog, $
          contour=0, /nocolor, charsize=charsize, $
          ztitle=ztitle, xtitle=xtitle,  $
          exact=1, aspect=0, cutaspect=0, ystyle=5, $
          barpos=1, zlen=-0.6, distbar=distbar, $ 
          barthick=barthick, position=[0.14,0.14,0.87,0.87]

cgAxis, YAxis=0, YRange=yrg, ystyle=1, /ylog, title='Period (s)', charsize=charsize
if mHz then cgAxis, YAxis=1, YRange=[1000./yrg[0],1000./yrg[1]], ystyle=1, /ylog, title='Frequency (mHz)', charsize=charsize $ else
    cgAxis, YAxis=1, YRange=[1./yrg[0],1./yrg[1]], ystyle=1, /ylog, title='Frequency (Hz)', charsize=charsize

; plot the Cone-of-Influence
plots, time, coi, noclip=0, linestyle=0, thick=coithick, color=cgColor('Black')
; shade the area above the Cone-of-Influence, with hashed lines:
ncoi = n_elements(coi)
y = fltarr(ncoi)
for j=0, ncoi-1 do y(j) = maxp
walsa_curvefill, time, y, coi, color = cgColor('Black'), thick=h_thick, /LINE_FILL, ORIENTATION=45 
walsa_curvefill, time, y, coi, color = cgColor('Black'), thick=h_thick, /LINE_FILL, ORIENTATION=-45

; contours mark significance level
cgContour, sigi, /noerase, levels=1., XTICKFORMAT="(A1)", YTICKFORMAT="(A1)", $
    xthick=1.e-40, ythick=1.e-40, xticklen=1.e-40, yticklen=1.e-40, xticks=1.e-40, yticks=1.e-40, $
    c_colors=[cgColor('Navy')], label=0, $
    c_linestyle=0, c_thick=c_thick

if EPS eq 1 then walsa_endeps, filename=epsfilename, /pdf

end

Cross Correlations: 1D power spectra

WaLSA_cross_spectrum

Calculating cross-spectrum (also known as co-spectrum or cross-power), coherence, and phase relationships between two time series, where the 1D power spectra are obtained with FFT (Fast Fourier Transform), Lomb-Scargle, Wavelet (global, oglobal, and sensible), and HHT (Hilbert-Huang Transform), using the WaLSA_speclizer.pro.

WaLSA_cross_spectrum.pro

; -----------------------------------------------------------------------------------------------------
; WaLSAtools: Wave analysis tools
; Copyright (C) 2025 WaLSA Team - Shahin Jafarzadeh et al.
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
; 
; http://www.apache.org/licenses/LICENSE-2.0
; 
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
; 
; Note: If you use WaLSAtools for research, please consider citing:
; Jafarzadeh, S., Jess, D. B., Stangalini, M. et al. 2025, Nature Reviews Methods Primers, in press.
; -----------------------------------------------------------------------------------------------------


;+
; NAME: WaLSA_cross_spectrum
;       part of -- WaLSAtools --
;
; PURPOSE:
; 
;   Calculate cross-spectral relationships of two time series whose amplitudes and powers
;             are computed using the WaLSA_speclizer routine.
;             The cross-spectrum is complex valued, thus its magnitude is returned as the
;             co-spectrum. The phase lags between the two time series are are estimated 
;             from the imaginary and real arguments of the complex cross spectrum.
;             The coherence is calculated from the normalized square of the amplitude 
;             of the complex cross-spectrum
;
; CALLING SEQUENCE:
;   walsa_cross_spectrum(data1=data1,data2=data2,time=time,phase_angle=phase_angle,coherence=coherence,cospectrum=cospectrum,frequencies=frequencies)
;
; + INPUTS:
;           
;   data1:    first (1D) time series
;   data2:    second (1D) time series, co-aligned with data1
;   time:     observing times in seconds (1D array)
;
; + OPTIONAL KEYWORDS:
; ----type of analysis----
;   fft:            if set, Fast Fourier Transform (FFT) power spectrum is computed: for regular (evenly sampled) time series.
;   lombscargle:    if set, Lomb-Scargle power spectrum is computed: for irregular (unevenly sampled) time series.
;   hht:            if set, a power spectrum associated to EMD + Hilbert Transform is computed: for regular (evenly sampled) time series.
; ----padding, detrending, and apodization parameters----
;   padding:        oversampling factor: zero padding (increasing timespan) to increase frequency resolution (NOTE: doesn't add information)
;   apod:           extent of apodization edges (of a Tukey window); default 0.1
;   nodetrendapod:  if set, neither detrending nor apodization is performed!
;   pxdetrend:      subtract linear trend with time per pixel. options: 1=simple, 2=advanced; default: 2
;   polyfit:        the degree of polynomial fit to the data to detrend it.
;                   if set, instead of linear fit this polynomial fit is performed.
;   meantemporal:   if set, only a very simple temporal detrending is performed by subtracting the mean signal from the signal.
;                   i.e., the fitting procedure (linear or higher polynomial degrees) is omitted.
;   meandetrend:    if set, subtract linear trend with time for the image means (i.e., spatial detrending)
;   recon:          optional keyword that will Fourier reconstruct the input timeseries.
;                   note: this does not preserve the amplitudes and is only useful when attempting 
;                   to examine frequencies that are far away from the 'untrustworthy' low frequencies.
; ----significance-level parameters----
;   siglevel:       significance level (default: 0.05 = 5% significance level = 95% confidence level)
;   nperm:          number of random permutations for the significance test -- the larger the better (default: 1000)
;   nosignificance: if set, no significance level is calculated.
; ----HHT parameters/options----
;   stdlimit:       standard deviation to be achieved before accepting an IMF (recommended value between 0.2 and 0.3; perhaps even smaller); default: 0.2
;   nfilter:        Hanning window width for two dimensional smoothing of the Hilbert spectrum. default: 3 
;                   (an odd integer, prefrabely equal to or larger than 3; equal to 0 to avoid the windowing)
;
;   n_segments:     number of euqal segments (to which both datasets are broken prior to the analyses; default: 1)
;                   Each of these segments is considered an independent realisation of the underlying process. 
;                   The cross spectrum for each segement are averaged together to provide phase and coherence estimates at each frequency.
;
; + OUTPUTS:
;
;   cospectrum:     co-spectrum, i.e., magnitude of the complex cross spectrum
;   frequencies:    an array of frequencies. unit: mHz
;   phase_angle:    phase angles
;   coherence:      coherence of two series
;   signif_cross:   significance levels for the cospectrum (1D array)
;   signif_coh:     significance levels for the coherence (1D array)
;   d1_power:       power spectrum of data1
;   d2_power:       power spectrum of data2
;
;  Shahin Jafarzadeh & David B. Jess | WaLSA Team
;  + some routines/recipe from CROSS_SPECTRUM.pro of Simon Vaughan
;-

function walsa_getcross_spectrum, data1, data2, time, phase_angle=phase_angle, coherence=coherence, frequencies=frequencies, $
                                     fft=fft, lombscargle=lombscargle, hht=hht,$ ; type of analysis
                                     padding=padding, apod=apod, nodetrendapod=nodetrendapod, pxdetrend=pxdetrend, meandetrend=meandetrend,$ ; padding and apodization parameters
                                     polyfit=polyfit, meantemporal=meantemporal, recon=recon,$
                                     stdlimit=stdlimit, nfilter=nfilter, n_segments=n_segments, d1_power=d1_power, d2_power=d2_power

  cadence=walsa_mode(walsa_diff(time))
  nt = n_elements(data1)
  ; number of segments
  if (n_elements(n_segments) eq 0) then n_segments=1 ; number of segments to break the time series into.
  mn = nt/n_segments
  n_cut = mn*n_segments
  x_1 = reform(data1[0:n_cut-1],mn,n_segments)
  x_2 = reform(data2[0:n_cut-1],mn,n_segments)
  xtime = reform(time[0:n_cut-1],mn,n_segments)


  frequencies = 1./(cadence*2)*findgen(mn/2+1)/(mn/2)
  nff=n_elements(frequencies)
  frequencies = frequencies[1:nff-1]
  frequencies = frequencies*1000. ; in mHz
  nf=n_elements(frequencies)

  amplitude1 = complexarr(nf,n_segments)
  amplitude2 = complexarr(nf,n_segments)

  for iseg=0L, n_segments-1 do begin

      power1s = walsa_speclizer(reform(x_1[*,iseg]),reform(xtime[*,iseg]),$ ; main inputs
                              frequencies=frequencies,amplitude=amplitude1s,$ ; main (additional) outputs
                              fft=fft, lombscargle=lombscargle, hht=hht,$ ; type of analysis
                              padding=padding, apod=apod, nodetrendapod=nodetrendapod, pxdetrend=pxdetrend, meandetrend=meandetrend,$ ; padding and apodization parameters
                              polyfit=polyfit,meantemporal=meantemporal,recon=recon,$
                              nosignificance=1,$ ; significance-level parameters
                              stdlimit=stdlimit, nfilter=nfilter, $ ; HHT parameters/options
                              mode=1, /silent) ; power calibration

      power2s = walsa_speclizer(reform(x_2[*,iseg]),reform(xtime[*,iseg]),$ ; main inputs
                              frequencies=frequencies,amplitude=amplitude2s,$ ; main (additional) outputs
                              fft=fft, lombscargle=lombscargle, hht=hht,$ ; type of analysis
                              padding=padding, apod=apod, nodetrendapod=nodetrendapod, pxdetrend=pxdetrend, meandetrend=meandetrend,$ ; padding and apodization parameters
                              polyfit=polyfit,meantemporal=meantemporal,recon=recon,$
                              nosignificance=1,$ ; significance-level parameters
                              stdlimit=stdlimit, nfilter=nfilter, $ ; HHT parameters/options
                              mode=1, /silent) ; power calibration

      amplitude1[*,iseg] = amplitude1s
      amplitude2[*,iseg] = amplitude2s
  endfor
  amplitude1 = reform(amplitude1)
  amplitude2 = reform(amplitude2)

  power1 = abs(amplitude1)^2
  power2 = abs(amplitude2)^2

  ; complex cross-power spectrum
  cross_power = amplitude1 * CONJ(amplitude2)
  ; co-spectrum (real parts of cross-power spectrum)
  cospectrum = ABS(cross_power)
; ----------------------------------------------------------
; Average over the ensamble of time series segments and adjacent frequencies
; average the second-order quantities: C, P_1, P_2 over the ensemble of M segments
  if (n_segments gt 1) then begin
      binC = total(cross_power,2)/float(n_segments)
      binP_1 = total(power1,2)/float(n_segments)
      binP_2 = total(power2,2)/float(n_segments)
  endif else begin
      binC = cross_power
      binP_1 = power1
      binP_2 = power2
  endelse
; ----------------------------------------------------------
; calculate coherence
  coherence = abs(binC)^2 / (binP_1 * binP_2)
; calculate the phase lag (phase of complex cross spectrum)
  phase_angle = atan(binC,/phase)*(180./!pi)
; ----------------------------------------------------------
  cospectrum = abs(binC)

  d1_power=binP_1
  d2_power=binP_2

return, cospectrum
end
;==================================================== MAIN ROUTINE ====================================================
pro walsa_cross_spectrum, data1=data1, data2=data2, time=time, phase_angle=phase_angle, coherence=coherence, frequencies=frequencies, cospectrum=cospectrum, $
                             fft=fft, lombscargle=lombscargle, hht=hht, welch=welch,$ ; type of analysis
                             padding=padding, apod=apod, nodetrendapod=nodetrendapod, pxdetrend=pxdetrend, meandetrend=meandetrend,$ ; padding and apodization parameters
                             polyfit=polyfit,meantemporal=meantemporal,recon=recon,resample_original=resample_original,nperm=nperm,siglevel=siglevel,$
                             stdlimit=stdlimit, nfilter=nfilter, $ ; HHT parameters/options
                             nosignificance=nosignificance, signif_coh=signif_coh, signif_cross=signif_cross, n_segments=n_segments, d1_power=d1_power, d2_power=d2_power

 if n_elements(nosignificance) eq 0 then nosignificance = 0 
 if n_elements(nperm) eq 0 then nperm = 50 
 if n_elements(siglevel) eq 0 then siglevel=0.05 ; 5% significance level = 95% confidence level

 time1 = time & time2 = time

 if n_elements(silent) eq 0 then silent=0

 sizecube1 = size(reform(data1))
 sizecube2 = size(reform(data2))

 givewarning = 0
 if sizecube1[0] eq 1 and sizecube1[0] eq 1 then begin 
     if sizecube1[1] ne sizecube1[1] then givewarning = 1
     if n_elements(time) ne sizecube1[1] then givewarning = 1
     if n_elements(time) ne sizecube2[1] then givewarning = 1
 endif else givewarning = 1
 if givewarning eq 1 then begin
      print, ' '
      print, ' [!] data1, data2, and time must be one diemnsional and have identical lengths.'
      print, ' '
      stop
 endif

 if silent eq 0 then begin
     cadence=walsa_mode(walsa_diff(time))
     temporal_Nyquist = 1. / (cadence * 2.)
     print,' '
     print,'The input datacubes are of size: ['+ARR2STR(sizecube1[1],/trim)+']'
     print,'Temporally, the important values are:'
     print,'    2-element duration (Nyquist period) = '+ARR2STR((cadence * 2.),/trim)+' seconds'
     print,'    Time series duration = '+ARR2STR(cadence*sizecube1[1],/trim)+' seconds'
     print,'    Nyquist frequency = '+ARR2STR(temporal_Nyquist*1000.,/trim)+' mHz'
     print, ' '
 endif

 cospectrum = walsa_getcross_spectrum(data1, data2, time, phase_angle=phase_angle, coherence=coherence, frequencies=frequencies, $
                                     fft=fft, lombscargle=lombscargle, hht=hht,$ ; type of analysis
                                     padding=padding, apod=apod, nodetrendapod=nodetrendapod, pxdetrend=pxdetrend, meandetrend=meandetrend,$ ; padding and apodization parameters
                                     polyfit=polyfit, meantemporal=meantemporal, recon=recon,$
                                     stdlimit=stdlimit, nfilter=nfilter, n_segments=n_segments, d1_power=d1_power, d2_power=d2_power)

 d1_power = d1_power/frequencies[0] ; in DN^2/mHz
 d2_power = d2_power/frequencies[0] ; in DN^2/mHz

 if nosignificance eq 0 then begin
     nf = n_elements(frequencies)
     ndata1 = n_elements(data1)
     dt1 = walsa_mode(walsa_diff(time1))
     ndata2 = n_elements(data2)
     dt2 = round(walsa_mode(walsa_diff(time2)))

     coh_perm = fltarr(nf,nperm)
     cross_perm = fltarr(nf,nperm)
     for ip=0L, nperm-1 do begin
         permutation1 = walsa_randperm(ndata1)
         y_perm1 = data1(permutation1)
         permutation2 = walsa_randperm(ndata2)
         y_perm2 = data2(permutation2)
         cospectrumsig = walsa_getcross_spectrum(y_perm1, y_perm2, time, coherence=cohsig, $
                                     fft=fft, lombscargle=lombscargle, hht=hht,$ ; type of analysis
                                     padding=padding, apod=apod, nodetrendapod=nodetrendapod, pxdetrend=pxdetrend, meandetrend=meandetrend,$ ; padding and apodization parameters
                                     polyfit=polyfit, meantemporal=meantemporal, recon=recon,$
                                     stdlimit=stdlimit, nfilter=nfilter, n_segments=n_segments)

         coh_perm[*,ip] = cohsig
         cross_perm[*,ip] = cospectrumsig
         if ip eq 0 then PRINT
         print,string(13b)+' >>> % Running Monte Carlo (significance): ',(ip*100.)/(nperm-1),format='(a,f4.0,$)'
     endfor
     PRINT
     PRINT
     signif_coh = walsa_confidencelevel(coh_perm, siglevel=siglevel, nf=nf)
     signif_cross = walsa_confidencelevel(cross_perm, siglevel=siglevel, nf=nf)
 endif

 print, ''
 print, 'COMPLETED!'
 print,''

end

Cross Correlations: Wavelet power spectra

WaLSA_wavelet_cross_spectrum

As a largely modified/extended variant of the wave_coherency.pro (of Torrence), this code calculates co-spectrum, coherence, and phase relationships between two time series, where the wavelet power spectra are obtained, thus cross-correlation parameters also have two dimensions.

WaLSA_wavelet_cross_spectrum.pro

; -----------------------------------------------------------------------------------------------------
; WaLSAtools: Wave analysis tools
; Copyright (C) 2025 WaLSA Team - Shahin Jafarzadeh et al.
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
; 
; http://www.apache.org/licenses/LICENSE-2.0
; 
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
; 
; Note: If you use WaLSAtools for research, please consider citing:
; Jafarzadeh, S., Jess, D. B., Stangalini, M. et al. 2025, Nature Reviews Methods Primers, in press.
; -----------------------------------------------------------------------------------------------------


;+
; NAME: WaLSA_wavelet_cross_spectrum
;       part of -- WaLSAtools --
;
; PURPOSE:
;   Compute (and optionally plot) wavelet cospectrum (cross-spectrum), coherence, and phase angles
;           between two time series as well as thier statistical significance levels.
;
; CALLING SEQUENCE:
;
;   WaLSA_wavelet_cross_spectrum, data1, time1, data2, time2, /plot
;
; + INPUTS:
;
;   data1:    first (1D) time series (evenly sampled)
;   data2:    second (1D) time series (evenly sampled), co-aligned with data1
;   time:     observing times in seconds (1D array)
;
; + OPTIONAL KEYWORDS:
;
; ----wavelet parameters/options----
;   mother:         wavelet function, providing different localisation/resolution in frequency and in time (also depends on param, m).
;                   currently, 'Morlet','Paul','DOG' (derivative of Gaussian) are available. default: 'Morlet'.
;   param:          optional mother wavelet parameter. 
;                   For 'Morlet' this is k0 (wavenumber), default is 6. 
;                   For 'Paul' this is m (order), default is 4.
;                   For 'DOG' this is m (m-th derivative), default is 2 (i.e., the real-valued Mexican-hat wavelet)
;   dj:             spacing between discrete scales. default: 0.025
;   colornoise:     if set, noise background is based on Auchère et al. 2017, ApJ, 838, 166 / 2016, ApJ, 825, 110
; ----significance-level parameters----
;   siglevel:       significance level (default: 0.05 = 5% significance level = 95% confidence level)
;   nperm:          number of random permutations for the significance test (default=50)
;                   note: the default value is set for quick tests. Choose a large number (e.g., 2000 or larger)
;                         for a better statistical result.
;   nosignificance: if set, the significance levels are calculated.
;                   (thus not overplotted as contours when plot option is set)
; ----padding, detrending, and apodization parameters----
;   padding:        oversampling factor: zero padding (increasing timespan) to increase frequency resolution (NOTE: doesn't add information)
;   apod:           extent of apodization edges (of a Tukey window); default 0.1
;   nodetrendapod:  if set, neither detrending nor apodization is performed!
;   pxdetrend:      subtract linear trend with time per pixel. options: 1=simple, 2=advanced; default: 2
;   polyfit:        the degree of polynomial fit to the data to detrend it.
;                   if set, instead of linear fit this polynomial fit is performed.
;   meantemporal:   if set, only a very simple temporal detrending is performed by subtracting the mean signal from the signal.
;                   i.e., the fitting procedure (linear or higher polynomial degrees) is omitted.
;   meandetrend:    if set, subtract linear trend with time for the image means (i.e., spatial detrending)
; ----plotting----
;   plot:           if set, wavelet power sepctra of the two time series as well as 
;                   their wavelet cospectrum (cross-spectrum) and coherence, along with the 
;                   significance levels as contours, are plotted.
;                   The phase angles between the two time series are also depicted by default. 
;                   Arrows pointing right mark zero phase (meaning in-phase oscillations),
;                   arrows pointing straight up indicate data2 lags behind data1 by 90 degrees.
;   noarrow:        if set, the phase angles are not overplotted as arrows.
;   arrowdensity:   number of arrows (iluustrating phase angles) in x and y directions (default: [30,18])
;   arrowsize:      size of the arrows (default: 1)
;   arrowheadsize:  size of the arrows' head (default: 1)
;   pownormal:      if set, the power is normalised to its maximum value
;   log:            if set, the power spectra and the cospectrum are plotted in log10 scale
;   removespace:    if set, the time-period areas affected by the coi over the entire time range are not plotted.
;   clt:            color table number (idl ctload)
;   koclt:          custom color tables (currently available: 1 and 2)
;
; + OUTPUTS:
;
;   cospectrum:             absolute values of the cross wavelet map
;   coherence:              wavelet coherence map, as a function of time and scale
;   phase_angle:            phase angles in degrees 
;   time:                   time vector, given by the overlap of time1 and time2 
;                           (it is not used: it is assumed the two time series are temporally aligned)
;   frequency:              the frequency vector; in mHz
;   scale:                  scale vector of scale indices, given by the overlap of scale1 and scale2 
;                           computed by WaLSA_wavelet.pro
;   coi:                    vector of the cone-of-influence
;   signif_coh:             significance map for the coherence (same 2D size as the coherence map)
;                           coherence/signif_coh indicates regions above the siglevel
;   signif_cross:           significance map for the cospectrum (same 2D size as the cospectrum map)
;                           cospectrum/signif_coh indicates regions above the siglevel
;   coh_global:             global (or mean) coherence averaged over all times
;   phase_global:           global (or mean) phase averaged over all times
;   cross_global:           global (or mean) cross wavelet averaged over all times
;   coh_oglobal:            global (or mean) coherence averaged over all times, excluding areas affected by COI (oglobal)
;   phase_oglobal:          global (or mean) phase averaged over all times, excluding areas affected by COI (oglobal)
;   cross_oglobal:          global (or mean) cross wavelet averaged over all times, excluding areas affected by COI (oglobal)
;;----------------------------------------------------------------------------
; This routine is originally based on WAVE_COHERENCY.pro
; Copyright (C) 1998-2005, Christopher Torrence
; This software may be used, copied, or redistributed as long as it is not
; sold and this copyright notice is reproduced on each copy made. This
; routine is provided as is without any express or
; implied warranties whatsoever.
;
; Reference: Torrence, C. and P. J. Webster, 1999: Interdecadal changes in the
;            ENSO-monsoon system. <I>J. Climate</I>, 12, 2679-2690.
;;----------------------------------------------------------------------------
; Largely modified/extended by Shahin Jafarzadeh 2016-2021
;-

function walsa_getcorss_spectrum_wavelet, $
         data1,time1,data2,time2, $
         mother=mother, param=param, dj=dj, colornoise=colornoise,$
         coherence=coherence,phase_angle=phase_angle, $
         time_out=time_out,scale_out=scale_out,coi_out=coi_out, $
         cross_wavelet=cross_wavelet,power1=power1,power2=power2, $
         frequency=frequency,period=period,silent=silent, $
         coh_global=coh_global, phase_global=phase_global, cross_global=cross_global, $
         coh_oglobal=coh_oglobal, phase_oglobal=phase_oglobal, cross_oglobal=cross_oglobal, $
         nosignificance=nosignificance,removespace=removespace, $
         nodetrendapod=nodetrendapod,log=log,plot=plot,clt=clt,koclt=koclt,$
         padding=padding,apod=apod,pxdetrend=pxdetrend,meandetrend=meandetrend,$
         polyfit=polyfit,meantemporal=meantemporal

    dt1 = walsa_mode(walsa_diff(time1))
    wave1 = walsa_wavelet(data1,dt1,scale=scale1,nodetrendapod=nodetrendapod,PERIOD=period,COI=coi1,plot=plot,w=4,log=log,silent=silent,$
                          removespace=removespace,koclt=koclt,clt=clt,mother=mother, param=param, dj=dj, colornoise=colornoise)

    frequency = 1000./period ; in mHz
    nf = n_elements(frequency)

    dt2 = walsa_mode(walsa_diff(time2))
    wave2 = walsa_wavelet(data2,dt2,scale=scale2,nodetrendapod=nodetrendapod,plot=plot,w=5,log=log,silent=silent,removespace=removespace,koclt=koclt,$
                          clt=clt,mother=mother, param=param, dj=dj, colornoise=colornoise)

;*** find overlapping times
    time_start = MIN(time1) > MIN(time2)
    time_end = MAX(time1) < MAX(time2)
    time1_start = MIN(WHERE((time1 ge time_start)))
    time1_end = MAX(WHERE((time1 le time_end)))
    time2_start = MIN(WHERE((time2 ge time_start)))
    time2_end = MAX(WHERE((time2 le time_end)))

;*** find overlapping scales
    scale_start = MIN(scale1) > MIN(scale2)
    scale_end = MAX(scale1) < MAX(scale2)
    scale1_start = MIN(WHERE((scale1 ge scale_start)))
    scale1_end = MAX(WHERE((scale1 le scale_end)))
    scale2_start = MIN(WHERE((scale2 ge scale_start)))
    scale2_end = MAX(WHERE((scale2 le scale_end)))

    period = period(scale1_start:scale1_end)

;*** cross wavelet & individual wavelet power
    cross_wavelet = wave1(time1_start:time1_end,scale1_start:scale1_end)*CONJ(wave2(time2_start:time2_end,scale2_start:scale2_end))    
    power1 = ABS(wave1(time1_start:time1_end,scale1_start:scale1_end))^2
    power2 = ABS(wave2(time2_start:time2_end,scale2_start:scale2_end))^2

    dt = dt1
    ntime = time1_end - time1_start + 1
    nj = scale1_end - scale1_start + 1
    if (N_EleMENTS(dj) le 0) then dj = ALOG(scale1(1)/scale1(0))/ALOG(2)
    scale = scale1(scale1_start:scale1_end)
    time_out = time1(time1_start:time1_end)
    scale_out = scale1(scale1_start:scale1_end)
    if (N_EleMENTS(coi1) EQ N_EleMENTS(time1)) then $
        coi_out = coi1(time1_start:time1_end)

    nt = n_elements(time_out)

; calculate global cross-power, coherency, and phase angle
; global wavelet is the time average of the wavelet spectrum
    global1 = TOTAL(power1,1,/nan)/nt
    global2 = TOTAL(power2,1,/nan)/nt
    cross_global = TOTAL(cross_wavelet,1)/nt
    coh_global = ABS(cross_global)^2/(global1*global2)
    phase_global = reform(ATAN(IMAGINARY(cross_global),REAL_PART(cross_global)))*(180./!pi)

    global1 = global1;/frequency[nf-1] ; in DN^2
    global2 = global2;/frequency[nf-1] ; in DN^2
    cross_global = ABS(cross_global);/frequency[nf-1] ; in DN^2

; calculate global cross-power, coherency, and phase angle excluding areas affected by COI (oglobal)
    oglobal_power1 = fltarr(nt,nf)
    oglobal_power2 = fltarr(nt,nf)
    oglobal_cross_wavelet = fltarr(nt,nf)
    for i=0L, nt-1 do begin
        ii = where(reform(period) lt coi_out[i], pnum)
        oglobal_power1[i,ii] = reform(power1[i,ii])
        oglobal_power2[i,ii] = reform(power2[i,ii])
        oglobal_cross_wavelet[i,ii] = reform(cross_wavelet[i,ii])
    endfor
    oglobal_global1 = TOTAL(oglobal_power1,1,/nan)/nt 
    oglobal_global2 = TOTAL(oglobal_power2,1,/nan)/nt 
    cross_oglobal = TOTAL(oglobal_cross_wavelet,1)/nt
    coh_oglobal = ABS(cross_oglobal)^2/(oglobal_global1*oglobal_global2)
    phase_oglobal = reform(ATAN(IMAGINARY(cross_oglobal),REAL_PART(cross_oglobal)))*(180./!pi)

    oglobal1 = oglobal_global1;/frequency[nf-1] ; in DN^2
    oglobal2 = oglobal_global2;/frequency[nf-1] ; in DN^2
    cross_oglobal = ABS(cross_oglobal);/frequency[nf-1] ; in DN^2

    for j=0,nj-1 do begin ;*** time-smoothing
        st1 = SYSTIME(1)
        nt = LONG(4L*scale(j)/dt)/2L*4 + 1L
        time_wavelet = (FINDgeN(nt) - nt/2)*dt/scale(j)
        wave_function = EXP(-time_wavelet^2/2.)   ;*** Morlet
        wave_function = FLOAT(wave_function/TOTAL(wave_function)) ; normalize
        nz = nt/2
        zeros = COMPleX(FltARR(nz),FltARR(nz))
        cross_wave_slice = [zeros,cross_wavelet(*,j),zeros]
        cross_wave_slice = CONVOL(cross_wave_slice,wave_function)
        cross_wavelet(*,j) = cross_wave_slice(nz:ntime+nz-1)
        zeros = FLOAT(zeros)
        power_slice = [zeros,power1(*,j),zeros]
        power_slice = CONVOL(power_slice,wave_function)
        power1(*,j) = power_slice(nz:ntime + nz - 1)
        power_slice = [zeros,power2(*,j),zeros]
        power_slice = CONVOL(power_slice,wave_function)
        power2(*,j) = power_slice(nz:ntime + nz - 1)
    endfor  ;*** time-smoothing

;*** normalize by scale
    scales = REBIN(TRANSPOSE(scale),ntime,nj)
    cross_wavelet = TEMPORARY(cross_wavelet)/scales
    power1 = TEMPORARY(power1)/scales
    power2 = TEMPORARY(power2)/scales

    nweights = FIX(0.6/dj/2 + 0.5)*2 - 1   ; closest (smaller) odd integer
    weights = REPLICATE(1.,nweights)
    weights = weights/TOTAL(weights) ; normalize
    for i=0,ntime-1 do begin ;*** scale-smoothing
        cross_wavelet(i,*) = CONVOL((cross_wavelet(i,*))(*),weights, /EDGE_TRUNCATE)
        power1(i,*) = CONVOL((power1(i,*))(*),weights, /EDGE_TRUNCATE)
        power2(i,*) = CONVOL((power2(i,*))(*),weights, /EDGE_TRUNCATE)
    endfor ;*** scale-smoothing

    wave_phase = reform(ATAN(IMAGINARY(cross_wavelet),REAL_PART(cross_wavelet)))*(180./!pi)
    wave_coher = (ABS(cross_wavelet)^2)/(power1*power2 > 1E-9)

    ; cospectrum = ABS(REAL_PART(cross_wavelet))
    cospectrum = ABS(cross_wavelet);/frequency[nf-1] ; in DN^2
    coherence=reform(wave_coher)
    phase_angle=reform(wave_phase)

 return, cospectrum

end
;==================================================== MAIN ROUTINE ====================================================
pro walsa_wavelet_cross_spectrum, $
    data1,data2,time, $   ;*** required inputs
    mother=mother, param=param, dj=dj, colornoise=colornoise, $
    coherence=coherence,phase_angle=phase_angle, $
    scale=scale,coi=coi, $
    coh_global=coh_global, phase_global=phase_global, cross_global=cross_global, $
    coh_oglobal=coh_oglobal, phase_oglobal=phase_oglobal, cross_oglobal=cross_oglobal, $
    cospectrum=cospectrum, period=period, $
    frequency=frequency,signif_coh=signif_coh,signif_cross=signif_cross, $
    padding=padding, apod=apod, nodetrendapod=nodetrendapod, pxdetrend=pxdetrend, meandetrend=meandetrend,$
    polyfit=polyfit,meantemporal=meantemporal,$
    nosignificance=nosignificance,pownormal=pownormal,siglevel=siglevel, $
    plot=plot,clt=clt,log=log,nperm=nperm,removespace=removespace,koclt=koclt,$
    arrowdensity=arrowdensity,arrowsize=arrowsize,arrowheadsize=arrowheadsize,noarrow=noarrow,silent=silent

    ; assuming the two time series are temporally aligned
    time1 = time
    time2 = time

    if n_elements(plot) eq 0 then plot = 0
    if n_elements(log) eq 0 then log = 0
    if n_elements(dj) eq 0 then dj = 0.025
    if n_elements(nosignificance) eq 0 then nosignificance = 0
    if n_elements(nodetrendapod) eq 0 then nodetrendapod = 0
    if n_elements(nperm) eq 0 then nperm = 50
    if n_elements(siglevel) eq 0 then siglevel=0.05 ; 5% significance level = 95% confidence level
    if n_elements(removespace) eq 0 then removespace=0
    if n_elements(silent) eq 0 then silent=0

    sizecube1 = size(reform(data1))
    sizecube2 = size(reform(data2))

    givewarning = 0
    if sizecube1[0] eq 1 and sizecube1[0] eq 1 then begin
        if sizecube1[1] ne sizecube1[1] then givewarning = 1
        if n_elements(time) ne sizecube1[1] then givewarning = 1
        if n_elements(time) ne sizecube2[1] then givewarning = 1
    endif else givewarning = 1
    if givewarning eq 1 then begin
        print, ' '
        print, ' [!] data1, data2, and time must be one diemnsional and have identical lengths.'
        print, ' '
        stop
    endif

    if silent eq 0 then begin
        cadence=walsa_mode(walsa_diff(time))
        temporal_Nyquist = 1. / (cadence * 2.)
        print,' '
        print,'The input datacubes are of size: ['+ARR2STR(sizecube1[1],/trim)+']'
        print,'Temporally, the important values are:'
        print,'    2-element duration (Nyquist period) = '+ARR2STR((cadence * 2.),/trim)+' seconds'
        print,'    Time series duration = '+ARR2STR(cadence*sizecube1[1],/trim)+' seconds'
        print,'    Nyquist frequency = '+ARR2STR(temporal_Nyquist*1000.,/trim)+' mHz'
        print, ' '
    endif

    cospectrum = walsa_getcorss_spectrum_wavelet(data1,time1,data2,time2,mother=mother,param=param,dj=dj,colornoise=colornoise,coherence=coherence,phase_angle=phase_angle, $
                                         TIME_OUT=time_out,SCAle_OUT=scale_out,COI_OUT=coi_out,CROSS_WAVEleT=cross_wavelet,POWER1=power1,POWER2=power2, $
                                         frequency=frequency,nosignificance=nosignificance,koclt=koclt, $
                                         log=log,period=period,plot=plot,clt=clt,removespace=removespace,$
                                         coh_global=coh_global, phase_global=phase_global, cross_global=cross_global, $
                                         coh_oglobal=coh_oglobal, phase_oglobal=phase_oglobal, cross_oglobal=cross_oglobal, $
                                         padding=padding, apod=apod, nodetrendapod=nodetrendapod, pxdetrend=pxdetrend, meandetrend=meandetrend,$
                                         polyfit=polyfit,meantemporal=meantemporal)

    nxx = n_elements(cospectrum[*,0])
    nyy = n_elements(cospectrum[0,*])

    if nosignificance eq 0 then begin
        ndata1 = n_elements(data1)
        dt1 = walsa_mode(walsa_diff(time1))
        ndata2 = n_elements(data2)
        dt2 = round(walsa_mode(walsa_diff(time2)))

        coh_perm = fltarr(nxx,nyy,nperm)
        cross_perm = fltarr(nxx,nyy,nperm)
        for ip=0L, nperm-1 do begin
            permutation1 = walsa_randperm(ndata1)
            y_perm1 = data1(permutation1)
            permutation2 = walsa_randperm(ndata2)
            y_perm2 = data2(permutation2)
            cospectrumsig = walsa_getcorss_spectrum_wavelet(y_perm1,time1,y_perm2,time2,coherence=cohsig, $
                                                 mother=mother,param=param,dj=dj,colornoise=colornoise, $
                                                 log=0,plot=0,silent=1,padding=padding,apod=apod,nodetrendapod=nodetrendapod, $
                                                 pxdetrend=pxdetrend, meandetrend=meandetrend,polyfit=polyfit,meantemporal=meantemporal)

            coh_perm[*,*,ip] = cohsig
            cross_perm[*,*,ip] = cospectrumsig
            if ip eq 0 then PRINT
            print,string(13b)+' >>> % Running Monte Carlo (significance): ',(ip*100.)/(nperm-1),format='(a,f4.0,$)'
        endfor
        PRINT
        PRINT
        signif_coh = walsa_confidencelevel_wavelet(coh_perm, siglevel=siglevel)
        signif_cross = walsa_confidencelevel_wavelet(cross_perm, siglevel=siglevel)
    endif

    if plot eq 1 then begin
        powplot = cospectrum
        perplot = period
        timplot = time_out
        coiplot = coi_out
        walsa_plot_wavelet_cross_spectrum, powplot, perplot, timplot, coiplot, clt=clt, w=8, phase_angle=phase_angle, log=log, normal=pownormal, $
                                           /crossspectrum, arrowdensity=arrowdensity,arrowsize=arrowsize,arrowheadsize=arrowheadsize, $
                                           noarrow=noarrow, significancelevel=signif_cross, nosignificance=nosignificance,removespace=removespace

        powplot = coherence
        perplot = period
        timplot = time_out
        coiplot = coi_out
        walsa_plot_wavelet_cross_spectrum, powplot, perplot, timplot, coiplot, clt=clt, w=9, phase_angle=phase_angle, log=0, normal=pownormal, $
                                           /coherencespectrum, arrowdensity=arrowdensity,arrowsize=arrowsize,arrowheadsize=arrowheadsize, $
                                           noarrow=noarrow, significancelevel=signif_coh, nosignificance=nosignificance,removespace=removespace
    endif

    frequency = 1000./period ; in mHz
    time = time_out
    coi = coi_out
    scale = scale_out

print,''
print, 'COMPLETED!'
print,''

end

This code also uses the following routine to plot the wavelet co-spectrum and coherence spectrum (along with confidence levels, cone-of-influence regions, and phase lags).

WaLSA_plot_wavelet_cross_spectrum.pro

; -----------------------------------------------------------------------------------------------------
; WaLSAtools: Wave analysis tools
; Copyright (C) 2025 WaLSA Team - Shahin Jafarzadeh et al.
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
; 
; http://www.apache.org/licenses/LICENSE-2.0
; 
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
; 
; Note: If you use WaLSAtools for research, please consider citing:
; Jafarzadeh, S., Jess, D. B., Stangalini, M. et al. 2025, Nature Reviews Methods Primers, in press.
; -----------------------------------------------------------------------------------------------------


function walsa_plot_vector, u, v, x, y, angle=angle, head_len=head_len,$
                 maxmag=maxmag, align=align, length=length,$
                 ref_text=ref_text, COLOR=color, THICK=thick, cstyle=cstyle, myphi=myphi, charsize=charsize

;Procedure to calculate and plot a vector
mylen=300.0 ; length of default arrow in pixels
rev=1.0    
x0=0.0
y0=0.0
x1=u/maxmag*mylen*length
y1=v/maxmag*mylen*length
dx=x1-x0
if (dx LT 0.0) then rev=-1.0
dy=y1-y0
r=SQRT(dx^2+dy^2)
theta=ATAN(dy/dx)     
phi=angle*!dtor
rfrac=head_len 
x2=x1-r*rfrac*rev*COS(theta-phi)
y2=y1-r*rfrac*rev*SIN(theta-phi)
x3=x1-r*rfrac*rev*COS(theta+phi)
y3=y1-r*rfrac*rev*SIN(theta+phi)
x4=x1-rfrac/2*r*rev*COS(theta)
y4=y1-rfrac/2*r*rev*SIN(theta)
;Calculate intersection of vector shaft and head points either
;side of the shaft - see
;http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d
;for more details
ua=((x3-x2)*(y0-y2)-(y3-y2)*(x0-x2))/$
   ((y3-y2)*(x1-x0)-(x3-x2)*(y1-y0))
x5=x0+ua*(x1-x0)
y5=y0+ua*(y1-y0) 

outputval = 1

;Plot vectors in data space - cstyle=0
;Position in device coordinates and then convert to data coordinates
if (cstyle EQ 0) then begin
    pt1=CONVERT_COORD(x,y, /DATA, /TO_DEVICE)
    xpts=[x0,x1,x2,x3,x4,x5]+pt1[0]-align*dx
    ypts=[y0,y1,y2,y3,y4,y5]+pt1[1]-align*dy
    pts=CONVERT_COORD(xpts,ypts, /DEVICE, /TO_DATA)
    xpts=pts[0,*]
    ypts=pts[1,*]
    x0=xpts[0]
    x1=xpts[1]
    x2=xpts[2]
    x3=xpts[3]
    x4=xpts[4]
    x5=xpts[5]
    y0=ypts[0]
    y1=ypts[1]
    y2=ypts[2]
    y3=ypts[3]
    y4=ypts[4]
    y5=ypts[5]   

    ; Plot the vectors omiting any vectors with NaNs
    z=[xpts, ypts]   
    if (TOTAL(FINITE(z)) EQ 12) then begin
        PLOTS, [x0,x5,x3,x1,x2,x5], [y0,y5,y3,y1,y2,y5], COLOR=color, THICK=thick, NOCLIP=0
        POLYFILL, [x2,x1,x3],[y2,y1,y3], COLOR=color, THICK=thick, NOCLIP=0
    endif

endif

;Plot reference vector - cstyle=1
;Position in device coordinates and then convert to data coordinates
if (cstyle EQ 1) then begin
    pt1=CONVERT_COORD(x,y, /NORMAL, /TO_DEVICE)
    xpts=[x0,x1,x2,x3,x4,x5]+pt1[0]
    ypts=[y0,y1,y2,y3,y4,y5]+pt1[1]
    x0=xpts[0]
    x1=xpts[1]
    x2=xpts[2]
    x3=xpts[3]
    x4=xpts[4]
    x5=xpts[5]
    y0=ypts[0]
    y1=ypts[1]
    y2=ypts[2]
    y3=ypts[3]
    y4=ypts[4]
    y5=ypts[5]     

    ; Plot the vectors omiting any vectors with NaNs
    z=[xpts, ypts]   
    if (TOTAL(FINITE(z)) EQ 12) then begin
        PLOTS, [x0,x5,x3,x1,x2,x5], [y0,y5,y3,y1,y2,y5], COLOR=color, THICK=thick, /DEVICE
        POLYFILL, [x2,x1,x3],[y2,y1,y3], COLOR=color, THICK=thick, /DEVICE
    endif

    ;Add the reference vector text
    xoffset = round(abs(x3-x2))/2.
    yoffset = round(abs(y3-y2))
    if myphi eq 0 then CGTEXT, x0, y0+(2.5*yoffset), cgGreek('phi')+'=0'+cgSymbol('deg'), ALIGNMENT=0.0, COLOR=color, /DEVICE, charsize=charsize
    if myphi eq 90 then CGTEXT, x0+(2.*xoffset), y0+(2.*xoffset), cgGreek('phi')+'=90'+cgSymbol('deg'), ALIGNMENT=0.0, COLOR=color, /DEVICE, charsize=charsize
endif

return, outputval

end     
;----------------------------------------------------------------------------
function walsa_vector, u,v,x,y, LENGTH = length,$
        Color=color, XSTRIDE=xstride, YSTRIDE=ystride, ALIGN=align, $
        REF_MAG=ref_mag, ANGLE=angle, HEAD_LEN=head_len,$
        REF_POS=ref_pos, REF_TEXT=ref_text, OVERPLOT=overplot, _EXTRA=extra, THICK=thick, charsize=charsize

a=SIZE(u)
b=SIZE(v)
c=SIZE(x)
d=SIZE(y)

;Initialise parameters if undefined
if (N_ELEMENTS(XSTRIDE) EQ 0) then xstride=0
if (N_ELEMENTS(YSTRIDE) EQ 0) then ystride=0
if N_ELEMENTS(LENGTH) EQ 0 then length=1.0
if N_ELEMENTS(COLOR) EQ 0 then color = !P.COLOR
if (N_ELEMENTS(ANGLE) EQ 0) then angle=22.5
if (N_ELEMENTS(HEAD_LEN) EQ 0) then head_len=0.3
if (N_ELEMENTS(TYPE) EQ 0) then TYPE=0
if (N_ELEMENTS(ALIGN) EQ 0) then align=0.5  
if (N_ELEMENTS(REF_TEXT) EQ 0) then ref_text=' '
if (N_ELEMENTS(REF_MAG) EQ 0) then begin
    maxmag=MAX(ABS(SQRT(u^2.+v^2.)))
endif ELSE begin
    maxmag=ref_mag
endELSE

outputval = 1

;Setup the plot area if undefined
if (NOT KEYWORD_SET(overplot)) then begin
    xs=x[0]-(x(1)-x(0))
    xf=x[N_ELEMENTS(x)-1]+(x(1)-x(0))
    ys=y[0]-(y(1)-y(0))
    yf=y[N_ELEMENTS(y)-1]+(y(1)-y(0))  
    PLOT,[xs,xf],[ys,yf], XSTYLE=1, YSTYLE=1, /NODATA,$
       COLOR=color, _EXTRA = extra
endif

;do stride data reduction if needed 
if (xstride GT 1) then begin
    mypts=FLTARR(a[1], a[2])
    mypts[*,*]=0.0     
    for iy=0,a[2]-1,xstride do begin
        for ix=0,a[1]-1,ystride do begin
            if ( ((ix/xstride) EQ FIX(ix/xstride)) AND $
               ((iy/ystride) EQ FIX(iy/ystride)) ) then mypts[ix,iy]=1.0
        endfor
    endfor
    pts=WHERE(mypts LT 1.0)
    u[pts]=0.0
    v[pts]=0.0
endif 

;Main vector plotting loop
for ix=1, N_ELEMENTS(x)-1 do begin
    for iy=1, N_ELEMENTS(y)-1 do begin
        tempt = walsa_plot_vector(u(ix,iy), v(ix,iy), x(ix), y(iy), $
               angle=angle, head_len=head_len,$
               maxmag=maxmag, align=align, length=length,$
               color=color, cstyle=0, THICK=thick)
    endfor
endfor

;Plot reference arrow(s)
if (N_ELEMENTS(REF_POS) NE 0) then begin
    tempt = walsa_plot_vector(2.*maxmag, 0.0, ref_pos[0], ref_pos[1], $
               angle=angle, ref_text=ref_text, head_len=head_len-0.2,$
               maxmag=maxmag, align=align, length=length,$
               color=color, cstyle=1, thick=thick, myphi=0, charsize=charsize)

    tempt = walsa_plot_vector(0.0, 2.*maxmag, ref_pos[0], ref_pos[1]+0.10, $
               angle=angle, ref_text=ref_text, head_len=head_len-0.2,$
               maxmag=maxmag, align=align, length=length,$
               color=color, cstyle=1, thick=thick, myphi=90, charsize=charsize)
endif

return, outputval

end
;----------------------------------------------------------------------------
pro walsa_plot_wavelet_cross_spectrum, power, period, time, coi, significancelevel=significancelevel, clt=clt, ylog=ylog, $
                                       phase_angle=phase_angle, log=log, crossspectrum=crossspectrum, normal=normal, epsfilename=epsfilename, $ 
                                       coherencespectrum=coherencespectrum, noarrow=noarrow, w=w, nosignificance=nosignificance, maxperiod=maxperiod, $
                                       arrowdensity=arrowdensity,arrowsize=arrowsize,arrowheadsize=arrowheadsize,removespace=removespace, koclt=koclt, arrowthick=arrowthick

if n_elements(crossspectrum) eq 0 then crossspectrum = 0
if n_elements(coherencespectrum) eq 0 then coherencespectrum = 0
if n_elements(log) eq 0 then log = 1
if n_elements(ylog) eq 0 then ylog = 1.
if n_elements(arrowdensity) eq 0 then arrowdensity = [30,18]
if n_elements(arrowsize) eq 0 then arrowsize = 1.
if n_elements(arrowheadsize) eq 0 then arrowheadsize = 1.
if n_elements(arrowthick) eq 0 then arrowthick = 2.
if n_elements(nosignificance) eq 0 then nosignificance = 0
if n_elements(noarrow) eq 0 then noarrow = 0
if n_elements(removespace) eq 0 then removespace=0
if n_elements(normal) eq 0 then normal=0
if n_elements(epsfilename) eq 0 then eps = 0 else eps = 1 

if crossspectrum eq 0 and coherencespectrum eq 0 then begin
    PRINT
    PRINT, ' Please define which one to plot: cross spectrum or coherence'
    PRINT
    stop
endif

if crossspectrum ne 0 and coherencespectrum ne 0 then begin
    PRINT
    PRINT, ' Please define which one to plot: cross spectrum or coherence'
    PRINT
    stop
endif

nt = n_elements(reform(time))
np = n_elements(reform(period))

fundf = 1000./(time[nt-1]) ; fundamental frequency (frequency resolution) in mHz
if n_elements(maxperiod) eq 0 then maxp = 1000./fundf else maxp = maxperiod ; longest period to be plotted
if n_elements(maxperiod) eq 0 then if removespace ne 0 then maxp = max(coi) ; remove areas below the COI
iit = closest_index(maxp,period)

period = period[0:iit]
if nosignificance eq 0 then isig = reform(significancelevel[*,0:iit])
power = reform(power[*,0:iit])
power = reverse(power,2)

np = n_elements(reform(period))

aphaseangle = phase_angle

aphaseangle = reform(aphaseangle[*,0:iit])

if nosignificance eq 0 then isig = reverse(isig,2)
if nosignificance eq 0 then sigi = power/isig

if n_elements(w) eq 0 then w=8

dimensions = GET_SCREEN_SIZE(RESOLUTION=resolution)
xscreensize = dimensions[0]
yscreensize = dimensions[1]
IF (xscreensize le yscreensize) THEN smallest_screensize = xscreensize
IF (yscreensize le xscreensize) THEN smallest_screensize = yscreensize

if EPS eq 1 then begin
    walsa_eps, size=[18,13]
    !p.font=0
    device,set_font='Times-Roman'
    charsize = 1.3
    !x.thick=4.
    !y.thick=4.
    !x.ticklen=-0.033
    !y.ticklen=-0.024
    ; !y.minor = 10
    barthick = 550
    distbar = 550
    coithick = 3.
    arrowsize = 20.*arrowsize
    arrowthick = (3.5/2.)*arrowthick
    c_thick = 3.
    h_thick = 1.4
    ; !y.tickinterval = 100.
    ; arrowheadsize = 10.
endif else begin
    if (xscreensize ge 1000) AND (yscreensize ge 1000) then begin 
        if crossspectrum ne 0 then window, w, xs=900, ys = 650, title=strtrim(w,2)+': Cross Wavelet Spectrum'
        if coherencespectrum ne 0 then window, w, xs=900, ys = 650, title=strtrim(w,2)+': Coherence'
        charsize=2.0
        !x.thick=2.  
        !y.thick=2.  
        !x.ticklen=-0.033
        !y.ticklen=-0.022
        ; !X.MINOR=6
        distbar = 30
        barthick = 30
        coithick = 2
        c_thick = 2.
        h_thick = 1.
    endif
    if (xscreensize lt 1000) OR  (yscreensize lt 1000) then begin 
        if crossspectrum ne 0 then window, w, xs=FIX(smallest_screensize*0.9), ys=FIX(smallest_screensize*0.9), title=strtrim(w,2)+': Cross Wavelet Spectrum'
        if coherencespectrum ne 0 then window, w, xs=FIX(smallest_screensize*0.9), ys=FIX(smallest_screensize*0.9), title=strtrim(w,2)+': Coherence'
        charsize=1.7
        !x.thick=2
        !y.thick=2
        !x.ticklen=-0.033
        !y.ticklen=-0.022
        distbar = 25
        barthick = 25
        coithick = 2
        c_thick = 2.
        h_thick = 1.
    endif
endelse

colset               
device,decomposed=0

xtitle = 'Time (s)'

if crossspectrum ne 0 then begin
    if normal ne 0 then begin
        ztitle = 'Normalised Cross Power!C'
        if log ne 0 then ztitle = 'Log!d10!n(Normalised Cross Power)!C'
    endif else begin
        ztitle = 'Cross Power!C'
        if log ne 0 then ztitle = 'Log!d10!n(Cross Power)!C'
    endelse
endif
if coherencespectrum ne 0 then begin
    ztitle = 'Coherence!C'
    if log ne 0 then ztitle = 'Log!d10!n(Coherence)!C'
endif

ii = where(power lt 0., cii)
if cii gt 0 then power(ii) = 0.
if crossspectrum ne 0 then if normal ne 0 then power = 100.*power/max(power)

xrg = minmax(time)
yrg = [max(period),min(period)]

; userlct,/full,verbose=0,coltab=nnn
if n_elements(clt) eq 0 then clt=20
loadct, clt
if n_elements(koclt) ne 0 then walsa_kopowercolor, koclt

if log ne 0 then power = alog10(power)

walsa_image_plot, power, xrange=xrg, yrange=yrg, $
          nobar=0, zrange=minmax(power,/nan), ylog=ylog, $
          contour=0, /nocolor, charsize=charsize, $
          ztitle=ztitle, xtitle=xtitle,  $
          exact=1, aspect=0, cutaspect=0, ystyle=5, $
          barpos=1, zlen=-0.6, distbar=distbar, $ 
          barthick=barthick, position=[0.14,0.14,0.87,0.87]

cgAxis, YAxis=0, YRange=yrg, ystyle=1, ylog=ylog, charsize=charsize, ytitle='Period (s)'

; Lblv = LOGLEVELS([max(period),min(period)])
; axlabel, Lblv, charsize=charsize, color=cgColor('Black') ,format='(i12)'

cgAxis, YAxis=1, YRange=[1000./yrg[0],1000./yrg[1]], ystyle=1, ylog=ylog, title='Frequency (mHz)', charsize=charsize


; plot phase angles as arrows
angle = aphaseangle
UU = cos(d2r(angle))
VV = sin(d2r(angle))

if noarrow eq 0 then $
tempt = walsa_vector(UU, VV, time, period, /overplot, color=cgColor('Black'), length=0.04*arrowsize, ySTRIDE=round(nt/float(ArrowDensity[0])), $
     xSTRIDE=round(np/float(ArrowDensity[1])), thick=arrowthick, head_len=0.5*arrowheadsize, ref_pos=[0.025,0.815], align=0.5, charsize=charsize)

; plot the Cone-of-Influence
plots, time, coi, noclip=0, linestyle=0, thick=coithick, color=cgColor('Black')
; shade the area above the Cone-of-Influence, with hashed lines:
ncoi = n_elements(coi)
y = fltarr(ncoi)
for j=0, ncoi-1 do y(j) = maxp
walsa_curvefill, time, y, coi, color = cgColor('Black'), thick=h_thick, /LINE_FILL, ORIENTATION=45
walsa_curvefill, time, y, coi, color = cgColor('Black'), thick=h_thick, /LINE_FILL, ORIENTATION=-45

; contours mark significance level
if nosignificance eq 0 then $
cgContour, sigi, /noerase, levels=1., XTICKforMAT="(A1)", YTICKforMAT="(A1)", $
     xthick=1.e-40, ythick=1.e-40, xticklen=1.e-40, yticklen=1.e-40, xticks=1.e-40, yticks=1.e-40, $
     c_colors=[cgColor('Navy')], label=0, $
     c_linestyle=0, c_thick=c_thick

if EPS eq 1 then walsa_endeps, filename=epsfilename, /pdf

end