Initial commit

.gitattributes

@@ -0,0 +1,4 @@
+# git files
+.gitattributes export-ignore
+.gitignore export-ignore
+.gitlab export-ignore

.gitignore

@@ -0,0 +1,39 @@
+###################
+# Compiled source #
+###################
+*.com
+*.class
+*.exe
+*.o
+*.so
+*.mexa64
+
+############
+# Packages #
+############
+*.7z
+*.dmg
+*.gz
+*.iso
+*.jar
+*.rar
+*.tar
+*.zip
+
+######################
+# Logs and databases #
+######################
+*.log
+*.sql
+*.sqlite
+
+######################
+# OS generated files #
+######################
+.DS_Store
+.DS_Store?
+._*
+.Spotlight-V100
+.Trashes
+ehthumbs.db
+Thumbs.db

README.md

@@ -0,0 +1,34 @@
+# MATLAB Discrete Trigonometric Transform Library
+## Author
+
+This library is written by Bradley Treeby, University College London. Contact: [email protected]
+
+## Overview
+
+This repository provides a MATLAB interface to the FFTW implementations of the discrete cosine transforms (DCT) and discrete sine transforms (DST). The functions are implemented in C++ as MATLAB mex functions.
+
+MATLAB natively includes an interface to the complex-to-complex transforms in FFTW via the inbuilt functions `fft`, `fft2`, and `fftn`. However, MATLAB does not include an interface to the real-to-real transforms which correspond to discrete trigonometric transforms (DTTs). This library is intended to fill the gap.
+
+Three functions are included: `dtt1D`, `dtt2D`, and `dtt3D`. These compute DTTs in 1D, 2D, and 3D. The function `dtt1D` can also perform 1D transformations over 2D arrays.
+
+The type of DTT is specified by the input `dtt_type`, where 1 to 4 corresponds to DCTs, and 5 to 8 to DSTs. Currently, only double precisions transforms are supported, thus the input array must be in double precision.
+
+The DTT functions in this library create the FFTW plan using FFTW_ESTIMATE, and then destroy the plan after execution. Consequently, the performance will not match that of the MATLAB inbuilt FFT functions which can re-use wisdom (see `help fftw`).
+
+## Compilation
+
+The mex files can be compiled using `compileDttMex`, which calls `mex`. In many cases, the default paths will need to be changed. Open `compileDttMex` for further details on how to setup and compile. 
+
+Precompiled mex files for Windows and macOS are included in the repository. These have been compiled using MATLAB 2019b. The Windows mex functions were compiled using Windows 10 (1803) and Microsoft Visual C++ 2015. The macOS mex functions were compiled using macOS Catalina (10.15.3) and Xcode Clang++ (11.3.1).
+
+## Examples
+
+An example of using `dtt1D` is included in the function `gradientDTT1D`. This computes a spectral gradient using any of the eight supported DTT symmetries, including an option for grid staggering. Several other example scripts are also included in the examples folder. 
+
+## License
+
+This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along with this program.  If not, see <https://www.gnu.org/licenses/>.

compileDttMex.m

@@ -0,0 +1,64 @@
+%COMPILEDTTMEX Compile mex-functions for dtt1D, dtt2D, and dtt3D.
+%
+% DESCRIPTION:
+%     compileDttMex compiles the mex functions for dtt1D, dtt2D, and dtt3D.
+%
+%     On Windows, the compilation uses a pre-compiled version of FFTW
+%     included in the repository. If re-compiling, .lib files can be
+%     generated from .dll files using a Visual Studio command prompt
+%
+%         lib /machine:x64 /def:libfftw3-3.def
+%     
+%     On Linux or macOS, fftw can be downloaded from http://www.fftw.org/.
+%     Compile using the following options:
+%
+%         sudo ./configure --enable-avx CFLAGS="-m64 -fPIC"
+%         sudo make
+%         sudo make install
+% 
+%     Note, use --enable-sse2 if avx instructions aren't supported on
+%     your processor.
+%
+%     It is assumed that a suitable C++ compiler is installed and selected
+%     by calling mex -setup. See: https://www.mathworks.com/support/...
+%     requirements/supported-compilers.html for more details. On linux,
+%     specify the version of the compiler and FFTW to use in the code
+%     below. On macOS, compiling mex files requires Xcode, which can be
+%     downloaded from the app store.   
+%
+% ABOUT:
+%     author        - Bradley Treeby
+%     date          - 28 September 2017
+%     last update   - 19 March 2020
+%
+% Copyright (C) 2017-2020 Bradley Treeby
+%
+% See also dtt1D, dtt2D, dtt3D
+
+% check for windows, mac, or linux
+if ispc
+
+    % use default compiler and link to pre-compiled FFTW library
+    mex -L"./" -llibfftw3-3 dtt1D.cpp
+    mex -L"./" -llibfftw3-3 dtt2D.cpp
+    mex -L"./" -llibfftw3-3 dtt3D.cpp
+    
+elseif ismac
+    
+    % use default compiler and link to FFTW installed on local machine
+    mex -I/usr/local/include/ -L/usr/local/lib -lfftw3 -lm dtt1D.cpp
+    mex -I/usr/local/include/ -L/usr/local/lib -lfftw3 -lm dtt2D.cpp
+    mex -I/usr/local/include/ -L/usr/local/lib -lfftw3 -lm dtt3D.cpp
+
+else
+    
+    % make user update paths based on examples below
+    error('Please update paths to GCC and FFTW before calling compileDttMex.');
+    
+%     % specify gcc compiler and dynamically link to FFTW (change the path
+%     % locations in the example below)
+%     mex -v GCC='/usr/bin/gcc-7' -I/usr/local/include/ -L/mnt/Apps/software/FFTW/3.3.8-OpenMPI-4.0.1-GCC-7.3.0-2.30/lib -lfftw3 -lm dtt1D.cpp
+%     mex -v GCC='/usr/bin/gcc-7' -I/usr/local/include/ -L/mnt/Apps/software/FFTW/3.3.8-OpenMPI-4.0.1-GCC-7.3.0-2.30/lib -lfftw3 -lm dtt2D.cpp
+%     mex -v GCC='/usr/bin/gcc-7' -I/usr/local/include/ -L/mnt/Apps/software/FFTW/3.3.8-OpenMPI-4.0.1-GCC-7.3.0-2.30/lib -lfftw3 -lm dtt3D.cpp
+
+end

dtt1D.cpp

@@ -0,0 +1,216 @@
+/**************************************************************************
+ * MEX file to compute 1D discrete trigonometric transforms in double 
+ * precision using FFTW. See dtt1D.m for usage notes.
+ *
+ * author: Bradley Treeby
+ * date: 31 May 2012
+ * last update: 19 March 2020
+ *
+ * Copyright (C) 2012-2020 Bradley Treeby
+ *
+ * This program is free software: you can redistribute it and/or modify it 
+ * under the terms of the GNU General Public License as published by the 
+ * Free Software Foundation, either version 3 of the License, or (at your 
+ * option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but 
+ * WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License along 
+ * with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *************************************************************************/
+
+#include <matrix.h>
+#include <mex.h>
+#include "fftw3.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+
+    //--------------------------------------------
+    // DECLARE VARIABLES
+    //--------------------------------------------
+
+    mxArray *output_mat;
+    const mwSize *dims;
+    mwSize check_el_num;
+    mwSize numelements;
+    int NX, NY, N, numdims;
+    int i;
+    double *input_ptr, *output_ptr;
+    fftw_plan plan;
+    int DIM = 1;    // set default DIM to 1 if not given by user
+    int rank = 1;   // DTTs are all 1D or rank 1
+    int n[1];       // DTT size is thus also rank 1
+    int howmany;
+    int idist, odist;
+    int istride, ostride;
+
+    //--------------------------------------------
+    // CHECK NUMBER OF INPUTS AND OUTPUTS
+    //--------------------------------------------
+
+    //check for proper number of input and output arguments
+    if( !(nrhs == 2 || nrhs == 3) ) {
+        mexErrMsgTxt("Two or three inputs are required.");
+	} else if(nlhs!=1) {
+        mexErrMsgTxt("One output is required.");
+	}
+
+    //--------------------------------------------
+    // CHECK AND ALLOCATE DTT TYPE INPUT
+    //--------------------------------------------
+
+    //check DTT_type input is real, scalar, and double precision
+    check_el_num = mxGetNumberOfElements(prhs[1]);
+    if( !(mxIsDouble(prhs[1]) && !mxIsComplex(prhs[1]) && check_el_num == 1) ){
+        mexErrMsgTxt("Input for DTT_TYPE must be real, scalar, and double precision.");
+    }
+
+    //get the DTT type (cast double input to 32-bit integer)
+    int DTT_type = (* (double *) mxGetPr(prhs[1]));
+
+    //check the value
+    if ( !((DTT_type == 1) || (DTT_type == 2) || (DTT_type == 3) || (DTT_type == 4) || (DTT_type == 5) || (DTT_type == 6) || (DTT_type == 7) || (DTT_type == 8)) ){
+        mexErrMsgTxt("Input for DTT_TYPE must be an integer between 1 and 8.");
+    }
+
+    //--------------------------------------------
+    // CHECK AND ALLOCATE DIM INPUT
+    //--------------------------------------------
+
+    //if DIM input is given, check it
+    if (nrhs == 3){
+
+        //get the number of elements
+        check_el_num = mxGetNumberOfElements(prhs[2]);
+
+        //check DIM input is real, scalar, and double precision
+        if( !(mxIsDouble(prhs[1]) && !mxIsComplex(prhs[1]) && check_el_num == 1) ){
+            mexErrMsgTxt("Input for DIM must be real, scalar, and double precision.");
+        }
+
+        //get the value of the DIM input (cast double input to 32-bit integer)
+        // (* something) means take the first element == something[0]
+        DIM = (int)(* (double *) mxGetPr(prhs[2]));
+
+        //check the value is 1 or 2
+        if ( !( (DIM == 1) || (DIM == 2) ) ){
+            mexErrMsgTxt("Input for DIM must be 1 or 2.");
+        }
+
+    }
+
+    //--------------------------------------------
+    // CHECK AND ALLOCATE INPUT AND OUTPUT ARRAYS
+    //--------------------------------------------
+
+    //check the input matrix is real and double precision
+    if( !(mxIsDouble(prhs[0]) && !mxIsComplex(prhs[0])) ) {
+        mexErrMsgTxt("Input array must be double precision and real.");
+    }  
+
+    //check that the input is either 1D or 2D
+    numdims = mxGetNumberOfDimensions(prhs[0]);     //number of dimensions is always >= 2
+    if ( numdims > 2){
+        mexErrMsgTxt("Input array must be 1D or 2D.");
+    }
+
+    //get the dimensions of the input array
+    dims = mxGetDimensions(prhs[0]);
+    numelements = mxGetNumberOfElements(prhs[0]);
+    NX = (int)dims[0];
+    NY = (int)dims[1];
+    N = NX*NY;
+
+    //check if 1D and force the correct DIM (input DIM is not used)
+    if (NX == 1) {
+        DIM = 2;
+    }
+    else if (NY == 1) {
+        DIM = 1;
+    }
+
+    //print dimensions of the input array
+    //mexPrintf("DTT Type %d, Array Dimensions %d by %d, DTT on Dimension %d\n", DTT_type, NX, NY, DIM);
+
+    //create MATLAB output
+    output_mat = plhs[0] = mxCreateNumericArray(numdims, dims, mxDOUBLE_CLASS, mxREAL);
+
+    //get pointer to input and output arrays
+    input_ptr = (double *) mxGetPr(prhs[0]);
+    output_ptr = (double *) mxGetPr(output_mat);
+
+    //--------------------------------------------
+    // DEFINE PLAN VARIABLES
+    //--------------------------------------------
+
+    //define plan_many variables based on DIM setting (these are switched 
+    //to c++ as MATLAB inputs are column major)
+    switch ( DIM ) {
+        case 1 :
+            //perform DTT over columns of input array
+            n[0] = NX;
+            howmany = NY;
+            idist = odist = NX;
+            istride = ostride = 1;            
+            break;
+        case 2 :      
+            //perform DTT over rows of input array
+            n[0] = NY;
+            howmany = NX;
+            idist = odist = 1;
+            istride = ostride = NX;
+            break;
+        default :
+            mexErrMsgTxt("Input for DIM must be 1 or 2.");
+    }
+
+    //plan_many requires a pointer to the input type. However, FFTW_REDFT00 etc
+    //are defined using #define so we cannot directly get a pointer to them. 
+    //Instead we allocate some memory, copy the value of the input type, and use the 
+    //pointer to this memory.
+    fftw_r2r_kind * fKind = (fftw_r2r_kind *) fftw_malloc(sizeof(fftw_r2r_kind));
+
+    //set DTT type
+    switch ( DTT_type ) {
+        case 1: fKind[0] = FFTW_REDFT00; break;
+        case 2: fKind[0] = FFTW_REDFT10; break;
+        case 3: fKind[0] = FFTW_REDFT01; break;
+        case 4: fKind[0] = FFTW_REDFT11; break;
+        case 5: fKind[0] = FFTW_RODFT00; break;
+        case 6: fKind[0] = FFTW_RODFT10; break;
+        case 7: fKind[0] = FFTW_RODFT01; break;
+        case 8: fKind[0] = FFTW_RODFT11; break;  		
+        default: mexErrMsgTxt("Input for DTT_TYPE must be an integer between 1 and 8.");
+    }
+
+    //--------------------------------------------
+    // CREATE AND EXECUTE FFTW PLAN
+    //--------------------------------------------
+
+    /*fftw_plan fftw_plan_many_r2r(int rank, const int *n, int howmany,
+                                  double *in, const int *inembed,
+                                  int istride, int idist,
+                                  double *out, const int *onembed,
+                                  int ostride, int odist,
+                                  const fftw_r2r_kind *kind, unsigned flags);*/
+
+    //create plan using the input and output pointers directly (out of place transform)
+    plan = fftw_plan_many_r2r(rank, (int *) &n, howmany, (double *) input_ptr, NULL, istride, idist, (double *) output_ptr, NULL, ostride, odist, (fftw_r2r_kind*) fKind, FFTW_ESTIMATE);    
+
+    //execute plan
+    fftw_execute(plan);    
+
+    //--------------------------------------------
+    // CLEANUP ALLOCATED VARIABLES
+    //--------------------------------------------
+
+    //cleanup
+    fftw_destroy_plan(plan);
+    fftw_free((fftw_r2r_kind*)fKind);
+
+    return;
+}