| 1 | ==============================================================================
|
|---|
| 2 |
|
|---|
| 3 | FFTReal
|
|---|
| 4 | Version 2.00, 2005/10/18
|
|---|
| 5 |
|
|---|
| 6 | Fourier transformation (FFT, IFFT) library specialised for real data
|
|---|
| 7 | Portable ISO C++
|
|---|
| 8 |
|
|---|
| 9 | (c) Laurent de Soras <[email protected]>
|
|---|
| 10 | Object Pascal port (c) Frederic Vanmol <[email protected]>
|
|---|
| 11 |
|
|---|
| 12 | ==============================================================================
|
|---|
| 13 |
|
|---|
| 14 |
|
|---|
| 15 |
|
|---|
| 16 | 1. Legal
|
|---|
| 17 | --------
|
|---|
| 18 |
|
|---|
| 19 | This library is free software; you can redistribute it and/or
|
|---|
| 20 | modify it under the terms of the GNU Lesser General Public
|
|---|
| 21 | License as published by the Free Software Foundation; either
|
|---|
| 22 | version 2.1 of the License, or (at your option) any later version.
|
|---|
| 23 |
|
|---|
| 24 | This library is distributed in the hope that it will be useful,
|
|---|
| 25 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 26 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|---|
| 27 | Lesser General Public License for more details.
|
|---|
| 28 |
|
|---|
| 29 | You should have received a copy of the GNU Lesser General Public
|
|---|
| 30 | License along with this library; if not, write to the Free Software
|
|---|
| 31 | Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
|---|
| 32 |
|
|---|
| 33 | Check the file license.txt to get full information about the license.
|
|---|
| 34 |
|
|---|
| 35 |
|
|---|
| 36 |
|
|---|
| 37 | 2. Content
|
|---|
| 38 | ----------
|
|---|
| 39 |
|
|---|
| 40 | FFTReal is a library to compute Discrete Fourier Transforms (DFT) with the
|
|---|
| 41 | FFT algorithm (Fast Fourier Transform) on arrays of real numbers. It can
|
|---|
| 42 | also compute the inverse transform.
|
|---|
| 43 |
|
|---|
| 44 | You should find in this package a lot of files ; some of them are of interest:
|
|---|
| 45 | - readme.txt: you are reading it
|
|---|
| 46 | - FFTReal.h: FFT, length fixed at run-time
|
|---|
| 47 | - FFTRealFixLen.h: FFT, length fixed at compile-time
|
|---|
| 48 | - FFTReal.pas: Pascal implementation (working but not up-to-date)
|
|---|
| 49 | - stopwatch directory
|
|---|
| 50 |
|
|---|
| 51 |
|
|---|
| 52 |
|
|---|
| 53 | 3. Using FFTReal
|
|---|
| 54 | ----------------
|
|---|
| 55 |
|
|---|
| 56 | Important - if you were using older versions of FFTReal (up to 1.03), some
|
|---|
| 57 | things have changed. FFTReal is now a template. Therefore use FFTReal<float>
|
|---|
| 58 | or FFTReal<double> in your code depending on the application datatype. The
|
|---|
| 59 | flt_t typedef has been removed.
|
|---|
| 60 |
|
|---|
| 61 | You have two ways to use FFTReal. In the first way, the FFT has its length
|
|---|
| 62 | fixed at run-time, when the object is instanciated. It means that you have
|
|---|
| 63 | not to know the length when you write the code. This is the usual way of
|
|---|
| 64 | proceeding.
|
|---|
| 65 |
|
|---|
| 66 |
|
|---|
| 67 | 3.1 FFTReal - Length fixed at run-time
|
|---|
| 68 | --------------------------------------
|
|---|
| 69 |
|
|---|
| 70 | Just instanciate one time a FFTReal object. Specify the data type you want
|
|---|
| 71 | as template parameter (only floating point: float, double, long double or
|
|---|
| 72 | custom type). The constructor precompute a lot of things, so it may be a bit
|
|---|
| 73 | long. The parameter is the number of points used for the next FFTs. It must
|
|---|
| 74 | be a power of 2:
|
|---|
| 75 |
|
|---|
| 76 | #include "FFTReal.h"
|
|---|
| 77 | ...
|
|---|
| 78 | long len = 1024;
|
|---|
| 79 | ...
|
|---|
| 80 | FFTReal <float> fft_object (len); // 1024-point FFT object constructed.
|
|---|
| 81 |
|
|---|
| 82 | Then you can use this object to compute as many FFTs and IFFTs as you want.
|
|---|
| 83 | They will be computed very quickly because a lot of work has been done in the
|
|---|
| 84 | object construction.
|
|---|
| 85 |
|
|---|
| 86 | float x [1024];
|
|---|
| 87 | float f [1024];
|
|---|
| 88 |
|
|---|
| 89 | ...
|
|---|
| 90 | fft_object.do_fft (f, x); // x (real) --FFT---> f (complex)
|
|---|
| 91 | ...
|
|---|
| 92 | fft_object.do_ifft (f, x); // f (complex) --IFFT--> x (real)
|
|---|
| 93 | fft_object.rescale (x); // Post-scaling should be done after FFT+IFFT
|
|---|
| 94 | ...
|
|---|
| 95 |
|
|---|
| 96 | x [] and f [] are floating point number arrays. x [] is the real number
|
|---|
| 97 | sequence which we want to compute the FFT. f [] is the result, in the
|
|---|
| 98 | "frequency" domain. f has the same number of elements as x [], but f []
|
|---|
| 99 | elements are complex numbers. The routine uses some FFT properties to
|
|---|
| 100 | optimize memory and to reduce calculations: the transformaton of a real
|
|---|
| 101 | number sequence is a conjugate complex number sequence: F [k] = F [-k]*.
|
|---|
| 102 |
|
|---|
| 103 |
|
|---|
| 104 | 3.2 FFTRealFixLen - Length fixed at compile-time
|
|---|
| 105 | ------------------------------------------------
|
|---|
| 106 |
|
|---|
| 107 | This class is significantly faster than the previous one, giving a speed
|
|---|
| 108 | gain between 50 and 100 %. The template parameter is the base-2 logarithm of
|
|---|
| 109 | the FFT length. The datatype is float; it can be changed by modifying the
|
|---|
| 110 | DataType typedef in FFTRealFixLenParam.h. As FFTReal class, it supports
|
|---|
| 111 | only floating-point types or equivalent.
|
|---|
| 112 |
|
|---|
| 113 | To instanciate the object, just proceed as below:
|
|---|
| 114 |
|
|---|
| 115 | #include "FFTRealFixLen.h"
|
|---|
| 116 | ...
|
|---|
| 117 | FFTRealFixLen <10> fft_object; // 1024-point (2^10) FFT object constructed.
|
|---|
| 118 |
|
|---|
| 119 | Use is similar as the one of FFTReal.
|
|---|
| 120 |
|
|---|
| 121 |
|
|---|
| 122 | 3.3 Data organisation
|
|---|
| 123 | ---------------------
|
|---|
| 124 |
|
|---|
| 125 | Mathematically speaking, the formulas below show what does FFTReal:
|
|---|
| 126 |
|
|---|
| 127 | do_fft() : f(k) = sum (p = 0, N-1, x(p) * exp (+j*2*pi*k*p/N))
|
|---|
| 128 | do_ifft(): x(k) = sum (p = 0, N-1, f(p) * exp (-j*2*pi*k*p/N))
|
|---|
| 129 |
|
|---|
| 130 | Where j is the square root of -1. The formulas differ only by the sign of
|
|---|
| 131 | the exponential. When the sign is positive, the transform is called positive.
|
|---|
| 132 | Common formulas for Fourier transform are negative for the direct tranform and
|
|---|
| 133 | positive for the inverse one.
|
|---|
| 134 |
|
|---|
| 135 | However in these formulas, f is an array of complex numbers and doesn't
|
|---|
| 136 | correspound exactly to the f[] array taken as function parameter. The
|
|---|
| 137 | following table shows how the f[] sequence is mapped onto the usable FFT
|
|---|
| 138 | coefficients (called bins):
|
|---|
| 139 |
|
|---|
| 140 | FFTReal output | Positive FFT equiv. | Negative FFT equiv.
|
|---|
| 141 | ---------------+-----------------------+-----------------------
|
|---|
| 142 | f [0] | Real (bin 0) | Real (bin 0)
|
|---|
| 143 | f [...] | Real (bin ...) | Real (bin ...)
|
|---|
| 144 | f [length/2] | Real (bin length/2) | Real (bin length/2)
|
|---|
| 145 | f [length/2+1] | Imag (bin 1) | -Imag (bin 1)
|
|---|
| 146 | f [...] | Imag (bin ...) | -Imag (bin ...)
|
|---|
| 147 | f [length-1] | Imag (bin length/2-1) | -Imag (bin length/2-1)
|
|---|
| 148 |
|
|---|
| 149 | And FFT bins are distributed in f [] as above:
|
|---|
| 150 |
|
|---|
| 151 | | | Positive FFT | Negative FFT
|
|---|
| 152 | Bin | Real part | imaginary part | imaginary part
|
|---|
| 153 | ------------+----------------+-----------------+---------------
|
|---|
| 154 | 0 | f [0] | 0 | 0
|
|---|
| 155 | 1 | f [1] | f [length/2+1] | -f [length/2+1]
|
|---|
| 156 | ... | f [...], | f [...] | -f [...]
|
|---|
| 157 | length/2-1 | f [length/2-1] | f [length-1] | -f [length-1]
|
|---|
| 158 | length/2 | f [length/2] | 0 | 0
|
|---|
| 159 | length/2+1 | f [length/2-1] | -f [length-1] | f [length-1]
|
|---|
| 160 | ... | f [...] | -f [...] | f [...]
|
|---|
| 161 | length-1 | f [1] | -f [length/2+1] | f [length/2+1]
|
|---|
| 162 |
|
|---|
| 163 | f [] coefficients have the same layout for FFT and IFFT functions. You may
|
|---|
| 164 | notice that scaling must be done if you want to retrieve x after FFT and IFFT.
|
|---|
| 165 | Actually, IFFT (FFT (x)) = x * length(x). This is a not a problem because
|
|---|
| 166 | most of the applications don't care about absolute values. Thus, the operation
|
|---|
| 167 | requires less calculation. If you want to use the FFT and IFFT to transform a
|
|---|
| 168 | signal, you have to apply post- (or pre-) processing yourself. Multiplying
|
|---|
| 169 | or dividing floating point numbers by a power of 2 doesn't generate extra
|
|---|
| 170 | computation noise.
|
|---|
| 171 |
|
|---|
| 172 |
|
|---|
| 173 |
|
|---|
| 174 | 4. Compilation and testing
|
|---|
| 175 | --------------------------
|
|---|
| 176 |
|
|---|
| 177 | Drop the following files into your project or makefile:
|
|---|
| 178 |
|
|---|
| 179 | Array.*
|
|---|
| 180 | def.h
|
|---|
| 181 | DynArray.*
|
|---|
| 182 | FFTReal*.cpp
|
|---|
| 183 | FFTReal*.h*
|
|---|
| 184 | OscSinCos.*
|
|---|
| 185 |
|
|---|
| 186 | Other files are for testing purpose only, do not include them if you just need
|
|---|
| 187 | to use the library ; they are not needed to use FFTReal in your own programs.
|
|---|
| 188 |
|
|---|
| 189 | FFTReal may be compiled in two versions: release and debug. Debug version
|
|---|
| 190 | has checks that could slow down the code. Define NDEBUG to set the Release
|
|---|
| 191 | mode. For example, the command line to compile the test bench on GCC would
|
|---|
| 192 | look like:
|
|---|
| 193 |
|
|---|
| 194 | Debug mode:
|
|---|
| 195 | g++ -Wall -o fftreal_debug.exe *.cpp stopwatch/*.cpp
|
|---|
| 196 |
|
|---|
| 197 | Release mode:
|
|---|
| 198 | g++ -Wall -o fftreal_release.exe -DNDEBUG -O3 *.cpp stopwatch/*.cpp
|
|---|
| 199 |
|
|---|
| 200 | It may be tricky to compile the test bench because the speed tests use the
|
|---|
| 201 | stopwatch sub-library, which is not that cross-platform. If you encounter
|
|---|
| 202 | any problem that you cannot easily fix while compiling it, edit the file
|
|---|
| 203 | test_settings.h and un-define the speed test macro. Remove the stopwatch
|
|---|
| 204 | directory from your source file list, too.
|
|---|
| 205 |
|
|---|
| 206 | If it's not done by default, you should activate the exception handling
|
|---|
| 207 | of your compiler to get the class memory-leak-safe. Thus, when a memory
|
|---|
| 208 | allocation fails (in the constructor), an exception is thrown and the entire
|
|---|
| 209 | object is safely destructed. It reduces the permanent error checking overhead
|
|---|
| 210 | in the client code. Also, the test bench requires Run-Time Type Information
|
|---|
| 211 | (RTTI) to be enabled in order to display the names of the tested classes -
|
|---|
| 212 | sometimes mangled, depending on the compiler.
|
|---|
| 213 |
|
|---|
| 214 | The test bench may take a long time to compile, especially in Release mode,
|
|---|
| 215 | because a lot of recursive templates are instanciated.
|
|---|
| 216 |
|
|---|
| 217 |
|
|---|
| 218 |
|
|---|
| 219 | 5. History
|
|---|
| 220 | ----------
|
|---|
| 221 |
|
|---|
| 222 | v2.00 (2005.10.18)
|
|---|
| 223 | - Turned FFTReal class into template (data type as parameter)
|
|---|
| 224 | - Added FFTRealFixLen
|
|---|
| 225 | - Trigonometric tables are size-limited in order to preserve cache memory;
|
|---|
| 226 | over a given size, sin/cos functions are computed on the fly.
|
|---|
| 227 | - Better test bench for accuracy and speed
|
|---|
| 228 |
|
|---|
| 229 | v1.03 (2001.06.15)
|
|---|
| 230 | - Thanks to Frederic Vanmol for the Pascal port (works with Delphi).
|
|---|
| 231 | - Documentation improvement
|
|---|
| 232 |
|
|---|
| 233 | v1.02 (2001.03.25)
|
|---|
| 234 | - sqrt() is now precomputed when the object FFTReal is constructed, resulting
|
|---|
| 235 | in speed impovement for small size FFT.
|
|---|
| 236 |
|
|---|
| 237 | v1.01 (2000)
|
|---|
| 238 | - Small modifications, I don't remember what.
|
|---|
| 239 |
|
|---|
| 240 | v1.00 (1999.08.14)
|
|---|
| 241 | - First version released
|
|---|
| 242 |
|
|---|