168 lines
7.0 KiB
Plaintext
168 lines
7.0 KiB
Plaintext
|
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"
|
||
|
|
Dec. 2006 : Fix a minor bug in "fftsg.f"
|
||
|
|
|