extension of class NumRu::GPhys -- Fast Fourier transformation and its applications

This manual documents the methods of NumRu::GPhys defined in gphys_fft.rb

class methods

GPhys::fft_ignore_missing( ignore=true, replace_val=nil )

Set a flag (class variable) to ignore missing values. This is for data that do not have missing but is treated as potentially having missing (often by having the valid_* attributes of NetCDF.) If replace_val is specified, data missing with replaced with that value.

methods

fft(backward=false, *dims)

Fast Fourier Transformation (FFT) by using (FFTW) ver 3 or ver 2. A FFTW ver.2 interface is included in NArray, while to use FFTW ver.3, you have to install separately. Dimension specification by the argument dims is available only with ver.3. By default, FFT is applied to all dimensions.

The transformation is complex. If the input data is not complex, it will be coerced to complex before transformation.

When the FT is forward, the result is normalized (i.e., divided by the data number), unlike the default behavior of FFTW.

Each coordinate is assumed to be equally spaced without checking. The new coordinate variables will be set equal to wavenumbers, derived as 2*PI/(length of the axis)*[0,1,2,..], where the length of the axis is derived as (coord.val.max - coord.val.min)*(n+1)/n.

REMARK

• If the units of the original coordinate is degree (or its equivalent ones such as degrees_east), the wavenumber was made in integers by converting the coordinate based on radian.

ARGUMENTS

• backward (true of false) : when true, backward FT is done; otherwise forward FT is done.
• dims (integers) : dimensions to apply FFT

RETURN VALUE

• a GPhys

EXAMPLE

gphy.fft           # forward, for all dimensions
gphy.fft(true)     # backward, for all dimensions
gphy.fft(nil, 0,1) # forward, for the first and second dimensions.
gphy.fft(true, -1) # backward, for the last dimension.

detrend(dim1[,dim2[,...]])

Remove means and linear trends along dimension(s) specified. Algorithm: 1st order polynomial fitting.

ARGUMENTS

• dim? (Integer of String): the dimension along which you want to remove trends.

RETURN VALUE

• a GPhys

EXAMPLE

cos_taper(dim1[,dim2[,...]])

Cosine tapering along dimension(s) specified.

Algorithm: to multiply with the half cosine curves at the both 1/10 ends of the data.

cos taper shape:
_____________
_/             \_
->   <-       ->   <-
T/10           T/10
half-cosine     half-cosine
shaped         shaped

The spectra of tapered data should be multiplied by 1/0.875, which is stored as GPhys::COS_TAPER_SP_FACTOR (==1/0.875).

ARGUMENTS

• dim? (Integer of String): the dimension along which you want to remove trends.

RETURN VALUE

• a GPhys

EXAMPLE

dim = 0    # for the 1st dimension
fc = gphys.detrend(dim).cos_taper(dim).fft(nil,dim)
sp = fc.abs**2 * GPhys::COS_TAPER_SP_FACTOR

spect_zero_centering(dim)

Shifts the wavenumber axis to cover from -K/2 to K/2 instead of from 0 to K-1, where the wavenumber is symbolically treated as integer, which is actually not the case, though. Since the first (-K/2) and the last (K/2) elements are duplicated, both are divided by 2. Therefore, this method is to be used for spectra (squared quantity) rather than the raw Fourier coefficients. (That is why the method name is prefixed by "spect_").

The method is applied for a single dimension (specified by the argument dim). If you need to apply for multiple dimensions, use it for multiple times.

ARGUMENTS

• dim (integer): the dimension you want to shift spectra elements. Count starts from zero.

RETURN VALUE

• a GPhys

EXAMPLE

• To get a spectra of a variable var along the 1st and 2nd dimensions:

fc = var.fft(nil, 0,1)    # --> Fourier coef
sp = ( fc.abs**2 ).spect_zero_centering(0).spect_zero_centering(1)

Note that spect_zero_centering is applied after taking |fc|^2.

• Same but if you want to have the 2nd dimension one-sided:

fc = var.fft(nil, 0,1)
sp = ( fc.abs**2 ).spect_zero_centering(0).spect_one_sided(1)
• Similar to the first example but for cross spectra:

fc1 = var1.fft(nil, 0,1)
fc2 = var2.fft(nil, 0,1)
xsp = (fc1 * fc2.conj).spect_zero_centering(0).spect_zero_centering(1)

spect_one_sided(dim)

Similar to spect_zero_centering but to make one-sided spectra. Namely, to convert from 0..K-1 to 0..K/2. To be applied for spectra; wavenumber 2..K/2-1 are multiplied by 2.

ARGUMENTS

• dim (integer): the dimension you want to shift spectra elements. Count starts from zero.

RETURN VALUE

• a GPhys

EXAMPLE

rawspect2powerspect(*dims)

Converts raw spectra obtained by gphys.fft.abs**2 into power spectra by dividing by wavenumber increments along the dimensions specified by dims.

ARGUMENTS

• dims (integers): the dimensions corresponding to wavenumbers.

RETURN VALUE

• a GPhys

EXAMPLE

• Suppose a 2 (or more) dimensional data gphys.

fc = gphys.fft(nil, 0, 1)
sp = fc.abs**2
ps = sp.rawspect2powerspect(0,1)

Here, sp is the raw spectrum of gphys, and ps is the power spectrum. The Parseval relation for them are as follows:

(gphys**2).mean == sp.sum
== pw.sum*dk*dl (== \int pw dk dl, mathematically),

where, dk = (pw.coord(0)[1] - pw.coord(0)[0]), and dl = (pw.coord(1)[1] - pw.coord(1)[0]).

phase_velocity_filter(xdim, tdim, cmin=nil, cmax=nil, xconv=nil, tconv=nil, remove_xtmean=false)

Filtering by phase velocity (between cmin and cmax)

REMARKS

• If the number of the grid points along x or t is an even number, the maximum wavenumber or frequency is treated as positive and negative, respectively, which results in an asymmetry of the treatment of positive and negative phase speeds. (That should be ok. -- In case its effect is significant, to do the filtering itself is not meaningful.)

ARGUMENTS

• xdim (Integer or String): spacial dimension
• tdim (Integer or String): time dimension
• cmin (Float or nil): minimum phase velocity. nil means no specification. (at least cmin or cmax must be given by Float)
• cmax (Float or nil): maximum phase velocity. nil means no specification. (at least cmin or cmax must be given by Float)
• xconv (nil or UNumeric) : (optional) if given, xconv is multiplied with the x axis before computing the phase velocity (kconv=1/xconv is used to scale wavenumbers)
• tconv (nil or UNumeric) : (optional) if given, tconv is multiplied with the t axis before computing the phase velocity (fconv=1/tconv is used to scale frequency)
• remove_xtmean (false or true) : if false (default), components with k=0 and f=0 are counted as c=0 (stationary), (unlike phase_velocity_binning), so they are included if cmin*cmax <= 0; if true, k=0 & f=0 components are always removed.

RETURN VALUE

• a GPhys

EXAMPLE

• For a 4D data with [x,y,z,t] dimensions, filtering by the phase velocity in the y dimension greater than 10.0 (in the unit of y/t) can be made by

cmin = 10.0; cmax = nil
gpfilt = gp.phase_velocity_filter(1, 3, cmin, cmax)
• For a global data (on the Earth's surface) with [lon, lat, z, time] axes, where the units of lon is "degrees" (or "degrees_east" or "radian") and the units of time is "hours", to filter disturbances whose zonal phase speed MEASURED AT THE EQUATOR is less or equal to 30 m/s can be made by

cmin = -30.0; cmax = 30.0
# This is a special case since "radian" is exceptionally omitted.
# See the private method __predefined_coord_units_conversion.
tconv = UNumeric[3.6e3, "s/hours"]
gpfilt = gp.phase_velocity_filter(1, 3, cmin, cmax, xconv, tconv)

phase_velocity_binning_iso_norml(kdim, fdim, cmin, cmax, cint, kconv=nil, fconv=nil)

Same as phase_velocity_binning but exclusively for equal phase velocity spacing. Also, a normalization is additionally made, to scale spectra in terms of integration along phase velocity axis --- The result of phase_velocity_binning called inside this method is divided by cint along with corresponding units conversion. Therefore, if this method is applied to spectra, a normalization is made such that an integration (not summation) along the phase velocity gives the variance (or covariance etc.) -- This normalization is suitable to quadratic quantities (such as spectra) but is not suitable to raw Fourier coefficients.

ARGUMENTS

RETURN VALUE

• a GPhys

phase_velocity_binning(kdim, fdim, cbins, kconv=nil, fconv=nil)

Bin a 2D spectrum in space and time based on phase velocity. The operand (self) must be Fourier coefficients or spectra, whose grid has not been altered since the call of the method fft (i.e., those that have not applied with zero centering etc, since it is done in this method).

Binning by this method is based on summation, leaving the units unchanged.

REMARKS

• Components whose phase velocities are exactly equal to one of the boundaries are divided into the two bins half by half
• components with k=0 and f=0 are excluded -- the spatio-temporal mean do not reflect in the result

ARGUMENTS

• kdim (Integer or String): wavenumber dimension (from spacial dimension)
• fdim (Integer or String): frequency dimension (from time dimension)
• cbins : an Array of bin bounds or a Hash of max, min, int e.g., [-10,-1,-0.1,0.1,11,10], {"min"=>-30,"max"=>30,"int"=>5}
• kconv (nil or UNumeric) : (optional) if given, kconv is multiplied with the wavenumber axis before computing the phase velocity
• fconv (nil or UNumeric) : (optional) if given, fconv is multiplied with the frequency axis before computing the phase velocity

RETURN VALUE

• a GPhys

EXAMPLES

• Example A

fu = u.fft(nil, 0, 2)
cfu = fu.phase_velocity_binning(0, 2, {"min"=>-1,"max"=>1,"int"=>0.1})
• Example B

fu = u.fft(nil, 0, 2)
pw = fu.abs**2rawspect2powerspect(0,2)          # power spectrum
cbins = [-100.0, -10.0, -1.0, 1.0, 10.0, 100.0] # logarithmic spacing
cpw = pw.phase_velocity_binning(0, 2, cbins)
• Example C

fu = u.fft(nil, 0, 3)
fv = v.fft(nil, 0, 3)
kconv = UNumeric[1/6.37e6, "m-1"]
fconv = UNumeric[1/3.6e3, "hours/s"]
fuv = (fu * fv.conj)   # cross spectra
cfuv = fuv.phase_velocity_binning(0, 3, {"min"=>-50,"max"=>50,"int"=>5},
kconv, fconv)