Cinder  0.8.6
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
fftsg.h
Go to the documentation of this file.
1 /*
2 
3 This is the OOURA implementation of various FFT algorithms, which was found at:
4 http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html
5 
6 The original package contains multiple implementations, we are using fftsg.c, also
7 called 'Fast Version III (Split-Radix)'. Documentation of the algorithms and their
8 parameters remain in the corresponding implementation file.
9 
10 Modifications (August 2013):
11 - The functions have been templated on the scalar type, which was originally double.
12 Definitions for T = float and double are declared at the bottom of fftsg.cpp
13 - Removed pthreads support.
14 
15 ========================================================================
16 
17 Original docs from OOURA's readme.txt:
18 
19 General Purpose FFT (Fast Fourier/Cosine/Sine Transform) Package
20 
21 Description:
22 A package to calculate Discrete Fourier/Cosine/Sine Transforms of
23 1-dimensional sequences of length 2^N.
24 
25 Files:
26 fft4g.c : FFT Package in C - Fast Version I (radix 4,2)
27 fft4g.f : FFT Package in Fortran - Fast Version I (radix 4,2)
28 fft4g_h.c : FFT Package in C - Simple Version I (radix 4,2)
29 fft8g.c : FFT Package in C - Fast Version II (radix 8,4,2)
30 fft8g.f : FFT Package in Fortran - Fast Version II (radix 8,4,2)
31 fft8g_h.c : FFT Package in C - Simple Version II (radix 8,4,2)
32 fftsg.c : FFT Package in C - Fast Version III (Split-Radix)
33 fftsg.f : FFT Package in Fortran - Fast Version III (Split-Radix)
34 fftsg_h.c : FFT Package in C - Simple Version III (Split-Radix)
35 readme.txt : Readme File
36 sample1/ : Test Directory
37 Makefile : for gcc, cc
38 Makefile.f77: for Fortran
39 testxg.c : Test Program for "fft*g.c"
40 testxg.f : Test Program for "fft*g.f"
41 testxg_h.c : Test Program for "fft*g_h.c"
42 sample2/ : Benchmark Directory
43 Makefile : for gcc, cc
44 Makefile.pth: POSIX Thread version
45 pi_fft.c : PI(= 3.1415926535897932384626...) Calculation Program
46 for a Benchmark Test for "fft*g.c"
47 
48 Difference of the Files:
49 C and Fortran versions are equal and
50 the same routines are in each version.
51 "fft4g*.*" are optimized for most machines.
52 "fft8g*.*" are fast on the UltraSPARC.
53 "fftsg*.*" are optimized for the machines that
54 have the multi-level (L1,L2,etc) cache.
55 The simple versions "fft*g_h.c" use no work area, but
56 the fast versions "fft*g.*" use work areas.
57 The fast versions "fft*g.*" have the same specification.
58 
59 Routines in the Package:
60 cdft: Complex Discrete Fourier Transform
61 rdft: Real Discrete Fourier Transform
62 ddct: Discrete Cosine Transform
63 ddst: Discrete Sine Transform
64 dfct: Cosine Transform of RDFT (Real Symmetric DFT)
65 dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
66 
67 Usage:
68 Please refer to the comments in the "fft**.*" file which
69 you want to use. Brief explanations are in the block
70 comments of each package. The examples are also given in
71 the test programs.
72 
73 Method:
74 -------- cdft --------
75 fft4g*.*, fft8g*.*:
76 A method of in-place, radix 2^M, Sande-Tukey (decimation in
77 frequency). Index of the butterfly loop is in bit
78 reverse order to keep continuous memory access.
79 fftsg*.*:
80 A method of in-place, Split-Radix, recursive fast
81 algorithm.
82 -------- rdft --------
83 A method with a following butterfly operation appended to "cdft".
84 In forward transform :
85 A[k] = sum_j=0^n-1 a[j]*W(n)^(j*k), 0<=k<=n/2,
86 W(n) = exp(2*pi*i/n),
87 this routine makes an array x[] :
88 x[j] = a[2*j] + i*a[2*j+1], 0<=j<n/2
89 and calls "cdft" of length n/2 :
90 X[k] = sum_j=0^n/2-1 x[j] * W(n/2)^(j*k), 0<=k<n.
91 The result A[k] are :
92 A[k] = X[k] - (1+i*W(n)^k)/2 * (X[k]-conjg(X[n/2-k])),
93 A[n/2-k] = X[n/2-k] +
94 conjg((1+i*W(n)^k)/2 * (X[k]-conjg(X[n/2-k]))),
95 0<=k<=n/2
96 (notes: conjg() is a complex conjugate, X[n/2]=X[0]).
97 -------- ddct --------
98 A method with a following butterfly operation appended to "rdft".
99 In backward transform :
100 C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n,
101 this routine makes an array r[] :
102 r[0] = a[0],
103 r[j] = Re((a[j] - i*a[n-j]) * W(4*n)^j*(1+i)/2),
104 r[n-j] = Im((a[j] - i*a[n-j]) * W(4*n)^j*(1+i)/2),
105 0<j<=n/2
106 and calls "rdft" of length n :
107 A[k] = sum_j=0^n-1 r[j]*W(n)^(j*k), 0<=k<=n/2,
108 W(n) = exp(2*pi*i/n).
109 The result C[k] are :
110 C[2*k] = Re(A[k] * (1-i)),
111 C[2*k-1] = -Im(A[k] * (1-i)).
112 -------- ddst --------
113 A method with a following butterfly operation appended to "rdft".
114 In backward transform :
115 S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n,
116 this routine makes an array r[] :
117 r[0] = a[0],
118 r[j] = Im((a[n-j] - i*a[j]) * W(4*n)^j*(1+i)/2),
119 r[n-j] = Re((a[n-j] - i*a[j]) * W(4*n)^j*(1+i)/2),
120 0<j<=n/2
121 and calls "rdft" of length n :
122 A[k] = sum_j=0^n-1 r[j]*W(n)^(j*k), 0<=k<=n/2,
123 W(n) = exp(2*pi*i/n).
124 The result S[k] are :
125 S[2*k] = Re(A[k] * (1+i)),
126 S[2*k-1] = -Im(A[k] * (1+i)).
127 -------- dfct --------
128 A method to split into "dfct" and "ddct" of half length.
129 The transform :
130 C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
131 is divided into :
132 C[2*k] = sum'_j=0^n/2 (a[j]+a[n-j])*cos(pi*j*k/(n/2)),
133 C[2*k+1] = sum_j=0^n/2-1 (a[j]-a[n-j])*cos(pi*j*(k+1/2)/(n/2))
134 (sum' is a summation whose last term multiplies 1/2).
135 This routine uses "ddct" recursively.
136 To keep the in-place operation, the data in fft*g_h.*
137 are sorted in bit reversal order.
138 -------- dfst --------
139 A method to split into "dfst" and "ddst" of half length.
140 The transform :
141 S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
142 is divided into :
143 S[2*k] = sum_j=1^n/2-1 (a[j]-a[n-j])*sin(pi*j*k/(n/2)),
144 S[2*k+1] = sum'_j=1^n/2 (a[j]+a[n-j])*sin(pi*j*(k+1/2)/(n/2))
145 (sum' is a summation whose last term multiplies 1/2).
146 This routine uses "ddst" recursively.
147 To keep the in-place operation, the data in fft*g_h.*
148 are sorted in bit reversal order.
149 
150 Reference:
151 * Masatake MORI, Makoto NATORI, Tatuo TORII: Suchikeisan,
152 Iwanamikouzajyouhoukagaku18, Iwanami, 1982 (Japanese)
153 * Henri J. Nussbaumer: Fast Fourier Transform and Convolution
154 Algorithms, Springer Verlag, 1982
155 * C. S. Burrus, Notes on the FFT (with large FFT paper list)
156 http://www-dsp.rice.edu/research/fft/fftnote.asc
157 
158 Copyright:
159 Copyright(C) 1996-2001 Takuya OOURA
160 email: ooura@mmm.t.u-tokyo.ac.jp
161 download: http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
162 You may use, copy, modify this code for any purpose and
163 without fee. You may distribute this ORIGINAL package.
164 
165 History:
166 ...
167 Dec. 1995 : Edit the General Purpose FFT
168 Mar. 1996 : Change the specification
169 Jun. 1996 : Change the method of trigonometric function table
170 Sep. 1996 : Modify the documents
171 Feb. 1997 : Change the butterfly loops
172 Dec. 1997 : Modify the documents
173 Dec. 1997 : Add "fft4g.*"
174 Jul. 1998 : Fix some bugs in the documents
175 Jul. 1998 : Add "fft8g.*" and delete "fft4f.*"
176 Jul. 1998 : Add a benchmark program "pi_fft.c"
177 Jul. 1999 : Add a simple version "fft*g_h.c"
178 Jul. 1999 : Add a Split-Radix FFT package "fftsg*.c"
179 Sep. 1999 : Reduce the memory operation (minor optimization)
180 Oct. 1999 : Change the butterfly structure of "fftsg*.c"
181 Oct. 1999 : Save the code size
182 Sep. 2001 : Add "fftsg.f"
183 Sep. 2001 : Add Pthread & Win32thread routines to "fftsg*.c"
184 Dec. 2006 : Fix a minor bug in "fftsg.f"
185 
186 
187 */
188 
189 #pragma once
190 
191 namespace cinder { namespace audio { namespace dsp { namespace ooura {
192 
193 template <typename T>
194 void cdft(int n, int isgn, T *a, int *ip, T *w);
195 
196 template <typename T>
197 void rdft(int n, int isgn, T *a, int *ip, T *w);
198 
199 template <typename T>
200 void ddct(int n, int isgn, T *a, int *ip, T *w);
201 
202 template <typename T>
203 void ddst(int n, int isgn, T *a, int *ip, T *w);
204 
205 template <typename T>
206 void dfct(int n, T *a, T *t, int *ip, T *w);
207 
208 template <typename T>
209 void dfst(int n, T *a, T *t, int *ip, T *w);
210 
211 } } } } // namespace cinder::audio::dsp::ooura
void dfct(int n, T *a, T *t, int *ip, T *w)
Definition: fftsg.cpp:627
void cdft(int n, int isgn, T *a, int *ip, T *w)
Definition: fftsg.cpp:435
void ddst(int n, int isgn, T *a, int *ip, T *w)
Definition: fftsg.cpp:566
GLenum GLsizei n
Definition: GLee.h:5780
void ddct(int n, int isgn, T *a, int *ip, T *w)
Definition: fftsg.cpp:505
void rdft(int n, int isgn, T *a, int *ip, T *w)
Definition: fftsg.cpp:456
void dfst(int n, T *a, T *t, int *ip, T *w)
Definition: fftsg.cpp:731
GLboolean GLboolean GLboolean GLboolean a
Definition: GLee.h:2964
GLubyte GLubyte GLubyte GLubyte w
Definition: GLee.h:2685
GLdouble GLdouble t
Definition: GLee.h:1426