Code Monkey home page Code Monkey logo

minfft's Introduction

minfft

A minimalistic Fast Fourier Transform library.

It achieves high performance by simple means.

Overview

The library routines compute:

  • Forward and inverse complex DFT,
  • Forward and inverse DFT of real data,
  • Cosine and sine transforms of types 2, 3, 4

of any dimensionality and power-of-two lengths.

The library provides C and Fortran interfaces.

Contents

Interface

All transform routines take three arguments:

  • input data pointer x,
  • output data pointer y,
  • auxiliary data pointer a.

The transform routines are capable of both in-place and out-of-place operation. In the latter case the input is left intact.

Auxiliary data contain chains of precomputed constants and temporary memory buffers, required for a transform routine to do its job. Once prepared, the auxiliary data can be reused as many times as needed. Also, the same auxiliary data fit for both forward and inverse transforms of the same kind.

These examples in C and Fortran show how the library functions are used:

	#include "minfft.h"
	minfft_cmpl x[N],y[N]; // input and output buffers
	minfft_aux *a; // aux data
	// prepare aux data
	a=minfft_mkaux_dft_1d(N);
	// do transforms
	minfft_dft(x,y,a);
	minfft_invdft(y,x,a);
	// free aux data
	minfft_free_aux(a);
	use minfft
	complex(minfft_cmpl),dimension(n) :: x,y ! input and output buffers
	type(minfft_aux) :: a ! aux data
	! prepare aux data
	a=minfft_mkaux_dft_1d(n)
	! do transforms
	call minfft_dft(x,y,a)
	call minfft_invdft(y,x,a)
	! free aux data
	call minfft_free_aux(a)

Data types

The library defines its own types for real and complex numbers, and for auxiliary data:

C Fortran
minfft_real real(minfft_real)
minfft_cmpl complex(minfft_cmpl)
minfft_aux type(minfft_aux)

By default, minfft_real is double. If C99 native complex is available, then minfft_cmpl is double complex. Otherwise, minfft_cmpl is an array of two minfft_reals. It is binary-compatible with other known complex number types, such as MSVC _Dcomplex or C++ complex<double>.

To build a single precision version, define the MINFFT_SINGLE macro.

Transforms

Below is a list of transforms with their definitions, auxiliary data makers, and transform routines.

For convenience, we provide aux data makers for one-, two- and three-dimensional transforms, along with a generic any-dimensional one. The dimensions are passed to aux maker routines in the C order (most rapidly varying index is the last). So, when calling from Fortran, array dimensions must be passed in the reverse order:

	complex(minfft_cmpl),dimension(n1,n2,n3) :: z
	a=minfft_mkaux_dft_3d(n3,n2,n1)

Auxiliary data makers return NULL if an error occured.

Our definitions of transforms, and formats of input and output data, are fully compatible with FFTW.

Complex DFT

$$ y_n=\sum_{k=0}^{N-1}x_k\exp(-\frac{2\pi i}{N}kn), ~ n=0,1,\dots,N{-}1 $$

minfft_aux* minfft_mkaux_dft_1d (int N);
minfft_aux* minfft_mkaux_dft_2d (int N1, int N2);
minfft_aux* minfft_mkaux_dft_3d (int N1, int N2, int N3);
minfft_aux* minfft_mkaux_dft (int d, int *Ns);
void minfft_dft (minfft_cmpl *x, minfft_cmpl *y, const minfft_aux *a);

Inverse complex DFT

$$ y_n=\sum_{k=0}^{N-1}x_k\exp(\frac{2\pi i}{N}kn), ~ n=0,1,\dots,N{-}1 $$

minfft_aux* minfft_mkaux_dft_1d (int N);
minfft_aux* minfft_mkaux_dft_2d (int N1, int N2);
minfft_aux* minfft_mkaux_dft_3d (int N1, int N2, int N3);
minfft_aux* minfft_mkaux_dft (int d, int *Ns);
void minfft_invdft (minfft_cmpl *x, minfft_cmpl *y, const minfft_aux *a);

Real DFT

This transform returns mostly non-redundant part of the complex DFT of real data.

For a real array $N_1 \times \dots \times N_{d-1} \times N_d$ it produces a complex array $N_1 \times \dots \times N_{d-1} \times (\frac{N_d}{2}+1)$.

Note that output takes a little more space than input. For in-place operation, make sure the data buffer is big enough to contain output.

minfft_aux* minfft_mkaux_realdft_1d (int N);
minfft_aux* minfft_mkaux_realdft_2d (int N1, int N2);
minfft_aux* minfft_mkaux_realdft_3d (int N1, int N2, int N3);
minfft_aux* minfft_mkaux_realdft (int d, int *Ns);
void minfft_realdft (minfft_real *x, minfft_cmpl *z, const minfft_aux *a);

Inverse real DFT

This is the inversion of the real DFT.

It takes a complex array $N_1 \times \dots \times N_{d-1} \times (\frac{N_d}{2}+1)$ and returns a real array $N_1 \times \dots \times N_{d-1} \times N_d$.

NB: Multidimensional inverse real DFT does not preserve input.

minfft_aux* minfft_mkaux_realdft_1d (int N);
minfft_aux* minfft_mkaux_realdft_2d (int N1, int N2);
minfft_aux* minfft_mkaux_realdft_3d (int N1, int N2, int N3);
minfft_aux* minfft_mkaux_realdft (int d, int *Ns);
void minfft_invrealdft (minfft_cmpl *z, minfft_real *y, const minfft_aux *a);

DCT-2

$$ y_n=2\sum_{k=0}^{N-1}x_k\cos(\frac{\pi}{N}(k+\frac{1}{2})n), ~ n=0,1,\dots,N{-}1 $$

minfft_aux* minfft_mkaux_t2t3_1d (int N);
minfft_aux* minfft_mkaux_t2t3_2d (int N1, int N2);
minfft_aux* minfft_mkaux_t2t3_3d (int N1, int N2, int N3);
minfft_aux* minfft_mkaux_t2t3 (int d, int *Ns);
void minfft_dct2 (minfft_real *x, minfft_real *y, const minfft_aux *a);

DST-2

$$ y_n=2\sum_{k=0}^{N-1}x_k\sin(\frac{\pi}{N}(k+\frac{1}{2})(n+1)), ~ n=0,1,\dots,N{-}1 $$

minfft_aux* minfft_mkaux_t2t3_1d (int N);
minfft_aux* minfft_mkaux_t2t3_2d (int N1, int N2);
minfft_aux* minfft_mkaux_t2t3_3d (int N1, int N2, int N3);
minfft_aux* minfft_mkaux_t2t3 (int d, int *Ns);
void minfft_dst2 (minfft_real *x, minfft_real *y, const minfft_aux *a);

DCT-3

$$ y_n=x_0+2\sum_{k=1}^{N-1}x_k\cos(\frac{\pi}{N}k(n+\frac{1}{2})), ~ n=0,1,\dots,N{-}1 $$

minfft_aux* minfft_mkaux_t2t3_1d (int N);
minfft_aux* minfft_mkaux_t2t3_2d (int N1, int N2);
minfft_aux* minfft_mkaux_t2t3_3d (int N1, int N2, int N3);
minfft_aux* minfft_mkaux_t2t3 (int d, int *Ns);
void minfft_dct3 (minfft_real *x, minfft_real *y, const minfft_aux *a);

DST-3

$$ y_n=2\sum_{k=0}^{N-2}x_k\sin(\frac{\pi}{N}(k+1)(n+\frac{1}{2}))+(-1)^{n}x_{N-1}, ~ n=0,1,\dots,N{-}1 $$

minfft_aux* minfft_mkaux_t2t3_1d (int N);
minfft_aux* minfft_mkaux_t2t3_2d (int N1, int N2);
minfft_aux* minfft_mkaux_t2t3_3d (int N1, int N2, int N3);
minfft_aux* minfft_mkaux_t2t3 (int d, int *Ns);
void minfft_dst3 (minfft_real *x, minfft_real *y, const minfft_aux *a);

DCT-4

$$ y_n=2\sum_{k=0}^{N-1}x_k\cos(\frac{\pi}{N}(k+\frac{1}{2})(n+\frac{1}{2})), ~ n=0,1,\dots,N{-}1 $$

minfft_aux* minfft_mkaux_t4_1d (int N);
minfft_aux* minfft_mkaux_t4_2d (int N1, int N2);
minfft_aux* minfft_mkaux_t4_3d (int N1, int N2, int N3);
minfft_aux* minfft_mkaux_t4 (int d, int *Ns);
void minfft_dct4 (minfft_real *x, minfft_real *y, const minfft_aux *a);

DST-4

$$ y_n=2\sum_{k=0}^{N-1}x_k\sin(\frac{\pi}{N}(k+\frac{1}{2})(n+\frac{1}{2})), ~ n=0,1,\dots,N{-}1 $$

minfft_aux* minfft_mkaux_t4_1d (int N);
minfft_aux* minfft_mkaux_t4_2d (int N1, int N2);
minfft_aux* minfft_mkaux_t4_3d (int N1, int N2, int N3);
minfft_aux* minfft_mkaux_t4 (int d, int *Ns);
void minfft_dst4 (minfft_real *x, minfft_real *y, const minfft_aux *a);

Freeing auxiliary data

If not needed anymore, the memory consumed by the auxiliary data can be freed by the minfft_free_aux() routine:

void minfft_free_aux (minfft_aux *a);

Memory requirements

Our library does not try to save memory, and allocates temporary buffers wherever it benefits performance.

The amounts of memory, allocated for the auxiliary data of the one-dimensional transforms, are given below:

Transform Auxiliary data size
Complex DFT of length N 2N complex numbers
Real DFT of length N 3.5N real numbers
Type-2 or Type-3 transform of length N 5.5N real numbers
Type-4 transform of length N 6N real numbers

Multi-dimensional transforms use a temporary buffer of the same size as the input data. This value is the dominant term in their auxiliary data size.

Implementation details

The complex DFT is computed by a split-radix (2/4), decimation in frequency, explicitly recursive fast Fourier transform. This method achieves a remarkable balance between performance and simplicity, and it behaves particularly cache-friendly, since it refers mostly to adjacent memory locations.

The real transforms are reduced eventually to a half-length complex transform.

For each transform, we first implement its one-dimensional, out-of-place, input-preserving, sequential input, strided output routine. This allows us to compute a multi-dimensional transform by repeated application of its one-dimensional routine along each dimension.

Performance

Below are the plots of speed and accuracy of our library, compared with the libraries of similar design โ€” KissFFT and PocketFFT. Performance of a highly optimized machine-specific version of the FFTW library is also shown for reference.

Standards

C89, Fortran 2003.

License

MIT.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.