fftsg.h

Raw

                  17
                 
                  135

This is the OOURA implementation of various FFT algorithms, which was found at:

http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html

The original package contains multiple implementations, we are using fftsg.c, also

called 'Fast Version III (Split-Radix)'. Documentation of the algorithms and their

parameters remain in the corresponding implementation file.

Modifications (August 2013):

- The functions have been templated on the scalar type, which was originally double.

Definitions for T = float and double are declared at the bottom of fftsg.cpp

- Removed pthreads support.

========================================================================

Original docs from OOURA's readme.txt:

General Purpose FFT (Fast Fourier/Cosine/Sine Transform) Package

Description:

A package to calculate Discrete Fourier/Cosine/Sine Transforms of

1-dimensional sequences of length 2^N.

Files:

fft4g.c : FFT Package in C - Fast Version I (radix 4,2)

fft4g.f : FFT Package in Fortran - Fast Version I (radix 4,2)

fft4g_h.c : FFT Package in C - Simple Version I (radix 4,2)

fft8g.c : FFT Package in C - Fast Version II (radix 8,4,2)

fft8g.f : FFT Package in Fortran - Fast Version II (radix 8,4,2)

fft8g_h.c : FFT Package in C - Simple Version II (radix 8,4,2)

fftsg.c : FFT Package in C - Fast Version III (Split-Radix)

fftsg.f : FFT Package in Fortran - Fast Version III (Split-Radix)

fftsg_h.c : FFT Package in C - Simple Version III (Split-Radix)

readme.txt : Readme File

sample1/ : Test Directory

Makefile : for gcc, cc

Makefile.f77: for Fortran

testxg.c : Test Program for "fft*g.c"

testxg.f : Test Program for "fft*g.f"

testxg_h.c : Test Program for "fft*g_h.c"

sample2/ : Benchmark Directory

Makefile : for gcc, cc

Makefile.pth: POSIX Thread version

pi_fft.c : PI(= 3.1415926535897932384626...) Calculation Program

for a Benchmark Test for "fft*g.c"

Difference of the Files:

C and Fortran versions are equal and

the same routines are in each version.

"fft4g*.*" are optimized for most machines.

"fft8g*.*" are fast on the UltraSPARC.

"fftsg*.*" are optimized for the machines that

have the multi-level (L1,L2,etc) cache.

The simple versions "fft*g_h.c" use no work area, but

the fast versions "fft*g.*" use work areas.

The fast versions "fft*g.*" have the same specification.

Routines in the Package:

cdft: Complex Discrete Fourier Transform

rdft: Real Discrete Fourier Transform

ddct: Discrete Cosine Transform

ddst: Discrete Sine Transform

dfct: Cosine Transform of RDFT (Real Symmetric DFT)

dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)

Usage:

Please refer to the comments in the "fft**.*" file which

you want to use. Brief explanations are in the block

comments of each package. The examples are also given in

the test programs.

Method:

-------- cdft --------

fft4g*.*, fft8g*.*:

A method of in-place, radix 2^M, Sande-Tukey (decimation in

frequency). Index of the butterfly loop is in bit

reverse order to keep continuous memory access.

fftsg*.*:

A method of in-place, Split-Radix, recursive fast

algorithm.

-------- rdft --------

A method with a following butterfly operation appended to "cdft".

In forward transform :

A[k] = sum_j=0^n-1 a[j]*W(n)^(j*k), 0<=k<=n/2,

W(n) = exp(2*pi*i/n),

this routine makes an array x[] :

x[j] = a[2*j] + i*a[2*j+1], 0<=j<n/2

and calls "cdft" of length n/2 :

X[k] = sum_j=0^n/2-1 x[j] * W(n/2)^(j*k), 0<=k<n.

The result A[k] are :

A[k] = X[k] - (1+i*W(n)^k)/2 * (X[k]-conjg(X[n/2-k])),

A[n/2-k] = X[n/2-k] +

conjg((1+i*W(n)^k)/2 * (X[k]-conjg(X[n/2-k]))),

0<=k<=n/2

(notes: conjg() is a complex conjugate, X[n/2]=X[0]).

-------- ddct --------

A method with a following butterfly operation appended to "rdft".

In backward transform :

C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n,

this routine makes an array r[] :

r[0] = a[0],

r[j] = Re((a[j] - i*a[n-j]) * W(4*n)^j*(1+i)/2),

r[n-j] = Im((a[j] - i*a[n-j]) * W(4*n)^j*(1+i)/2),

0<j<=n/2

and calls "rdft" of length n :

A[k] = sum_j=0^n-1 r[j]*W(n)^(j*k), 0<=k<=n/2,

W(n) = exp(2*pi*i/n).

The result C[k] are :

C[2*k] = Re(A[k] * (1-i)),

C[2*k-1] = -Im(A[k] * (1-i)).

-------- ddst --------

A method with a following butterfly operation appended to "rdft".

In backward transform :

S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n,

this routine makes an array r[] :

r[0] = a[0],

r[j] = Im((a[n-j] - i*a[j]) * W(4*n)^j*(1+i)/2),

r[n-j] = Re((a[n-j] - i*a[j]) * W(4*n)^j*(1+i)/2),

0<j<=n/2

and calls "rdft" of length n :

A[k] = sum_j=0^n-1 r[j]*W(n)^(j*k), 0<=k<=n/2,

W(n) = exp(2*pi*i/n).

The result S[k] are :

S[2*k] = Re(A[k] * (1+i)),

S[2*k-1] = -Im(A[k] * (1+i)).

-------- dfct --------

A method to split into "dfct" and "ddct" of half length.

The transform :

C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n

is divided into :

C[2*k] = sum'_j=0^n/2 (a[j]+a[n-j])*cos(pi*j*k/(n/2)),

C[2*k+1] = sum_j=0^n/2-1 (a[j]-a[n-j])*cos(pi*j*(k+1/2)/(n/2))

(sum' is a summation whose last term multiplies 1/2).

This routine uses "ddct" recursively.

To keep the in-place operation, the data in fft*g_h.*

are sorted in bit reversal order.

-------- dfst --------

A method to split into "dfst" and "ddst" of half length.

The transform :

S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n

is divided into :

S[2*k] = sum_j=1^n/2-1 (a[j]-a[n-j])*sin(pi*j*k/(n/2)),

S[2*k+1] = sum'_j=1^n/2 (a[j]+a[n-j])*sin(pi*j*(k+1/2)/(n/2))

(sum' is a summation whose last term multiplies 1/2).

This routine uses "ddst" recursively.

To keep the in-place operation, the data in fft*g_h.*

are sorted in bit reversal order.

Reference:

* Masatake MORI, Makoto NATORI, Tatuo TORII: Suchikeisan,

Iwanamikouzajyouhoukagaku18, Iwanami, 1982 (Japanese)

* Henri J. Nussbaumer: Fast Fourier Transform and Convolution

Algorithms, Springer Verlag, 1982

* C. S. Burrus, Notes on the FFT (with large FFT paper list)

http://www-dsp.rice.edu/research/fft/fftnote.asc

History:

...

Dec. 1995 : Edit the General Purpose FFT

Mar. 1996 : Change the specification

Jun. 1996 : Change the method of trigonometric function table

Sep. 1996 : Modify the documents

Feb. 1997 : Change the butterfly loops

Dec. 1997 : Modify the documents

Dec. 1997 : Add "fft4g.*"

Jul. 1998 : Fix some bugs in the documents

Jul. 1998 : Add "fft8g.*" and delete "fft4f.*"

Jul. 1998 : Add a benchmark program "pi_fft.c"

Jul. 1999 : Add a simple version "fft*g_h.c"

Jul. 1999 : Add a Split-Radix FFT package "fftsg*.c"

Sep. 1999 : Reduce the memory operation (minor optimization)

Oct. 1999 : Change the butterfly structure of "fftsg*.c"

Oct. 1999 : Save the code size

Sep. 2001 : Add "fftsg.f"

Sep. 2001 : Add Pthread & Win32thread routines to "fftsg*.c"

Dec. 2006 : Fix a minor bug in "fftsg.f"

# pragma once

namespace cinder { namespace audio { namespace dsp { namespace ooura {

template < typename T>

void cdft ( int n, int isgn, T *a, int *ip, T *w);

template < typename T>

void rdft ( int n, int isgn, T *a, int *ip, T *w);

template < typename T>

void ddct ( int n, int isgn, T *a, int *ip, T *w);

template < typename T>

void ddst ( int n, int isgn, T *a, int *ip, T *w);

template < typename T>

void dfct ( int n, T *a, T *t, int *ip, T *w);

template < typename T>

void dfst ( int n, T *a, T *t, int *ip, T *w);

} } } } // namespace cinder::audio::dsp::ooura

fftsg.h

Latest commit

History

fftsg.h

File metadata and controls