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
Copyright:
Copyright(C) 1996-2001 Takuya OOURA
email: ooura@mmm.t.u-tokyo.ac.jp
download: http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
You may use, copy, modify this code for any purpose and
without fee. You may distribute this ORIGINAL package.
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"
评论3