1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
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"