JuliaLang/julia — pull request #1878 by stevengj (2013-01-03 00:42:18+00:00)
Title: add dct functions and interface to FFTW's r2r transforms

This patch adds access to FFTW's r2r transforms, in order to easily perform multidimensional DCTs, DSTs, and related transforms. The interface is similar to `fft`: `FFTW.r2r(A, kind [, dims])` performs the multidimensional r2r transform of the given `kind` (`FFTW.REDFT00` etcetera as defined in the FFTW manual), optionally along a subset of the dimensions. `kind` may be an array (or tuple) in order to specify different transforms along different dimensions. There is also `r2r!` for in-place transforms, `plan_r2r` for pre-optimized transforms, etcetera.

It also adds a `dct(A [, dims])` function, based on the r2r transforms, which computes a type-II DCT with a unitary normalization, identical to Matlab's `dct` function, along with the inverse `idct` (type-III DCT). Unlike Matlab, it computes a multidimensional DCT by default, with the Matlab behavior obtained via `dct(A, 1)`. And of course there are in-place (`dct!`) and pre-planned (`plan_dct`) variants. This special case is important enough, and tricky enough to compute on top of r2r if you care about normalization, that it is probably worth including in Julia.

Matlab also has a `dst` (discrete sine transform) function, but after some thought I chose not to write an equivalent (although it would be easy). Matlab's DST interface seems poorly thought-out to me. They choose a type-I DST, when for a variety of reasons I would choose a type-II DST if you have to pick one. The main application is for solving PDEs, but they don't provide a multidimensional DST. And for solving PDEs, you probably want to support a variety of boundary conditions (e.g. Dirichlet and Neumann) anyway, so you are probably better off with the lower-level r2r interface anyway.

There are a variety of other functions in Matlab's signal-processing toolbox, like the cepstrum and the fast Hadamard transform, that would be straightforward to provide on top of FFTW. There is also the MDCT, which is notably omitted from Matlab but is easy to implement on top of a type-IV DCT (`REDFT11` in FFTW parlance). I don't know if you would want such things in Julia, however, or as a separate add-on package.